Orthorhombic distortion and rectangular skyrmion crystal in a centrosymmetric tetragonal host

We theoretically investigate the stability of a rectangular skyrmion crystal without fourfold rotational symmetry under an orthorhombic distortion in centrosymmetric tetragonal magnets. The results are obtained by numerically simulated annealing for an effective spin model with competing interactions in momentum space and magnetic anisotropy. By constructing the low-temperature phase diagram while changing the interaction ratio arising from the orthorhombic distortion, we find that the rectangular skyrmion crystal remains stable in an external magnetic field against distortion. We show that the degree of fourfold rotational symmetry tends to recover when the magnetic field is increased. The relevance to the skyrmion-hosting material EuAl4 is also discussed.

The spin textures in the SkX are described by a superposition of multiple spin density waves, i.e. a multiple-Q state [38][39][40][41]. When a superposition of symmetry-related wave vectors dictated by the lattice structures lowers the energy compared to that of the single-Q (1Q) spiral state, the triple-Q triangular SkX appears in the threefold-symmetric lattice structures, such as cubic, hexagonal, and trigonal systems, and the double-Q (2Q) square SkX (S-SkX) appears in the fourfold-symmetric tetragonal lattice structures. In these situations, the SkX is invariant under threefold or fourfold rotation. Meanwhile, it was revealed that SkXs remain stable in the low-symmetric lattice structures [42][43][44][45][46][47] and nonsymmorphic lattice structures with the screw axis [48]. In addition, a recent experiment has shown that a rhombic SkX was observed in the centrosymmetric tetragonal magnet EuAl 4 [49][50][51][52][53], where the importance of a spin-lattice coupling through a tetragonal-to-orthorhombic structural distortion was indicated in realizing such a low-symmetric SkX [54,55]. Based on this observation, it is important to clarify the stability of the SkX in the low-symmetric lattice structures without a threefold or fourfold rotational axis.
In the present study, we numerically investigate the stability of the SkX against orthorhombic distortion in centrosymmetric tetragonal magnets. Starting from an effective spin model on a square lattice, which was shown to exhibit the instabilities toward a fourfold-symmetric S-SkX and a fourfold-asymmetric rectangular SkX (R-SkX) depending on an external magnetic field, we examine the effect of the uniaxial distortion on these SkXs. Through the simulated annealing, we find that the R-SkX is robustly stabilized against the orthorhombic distortion when the competing interactions in momentum space are present. Although the S-SkX is replaced by the R-SkX once the orthorhombic distortion is turned on, the fourfold-symmetric nature tends to recover in the high-field region. We also discuss the relevance to EuAl 4 . Our result indicates that the low-symmetric lattice structures without more than a twofold rotational axis can be platforms to host the SkX in an external magnetic field.

Model and method
Let us start by considering the following effective spin model on a fourfold-symmetric square lattice in two dimensions, whose Hamiltonian is given by where The spin length is fixed at |S i | = 1, where we use reduced units for S = 7/2 of the Eu 2+ ion, and the lattice constant of the square lattice is taken as unity. The first term represents the bilinear exchange interaction in the Q ν channel. Although Q ν generally represents all the wave vectors in the Brillouin zone, we only consider a few wave vectors that give dominantly important contributions to the interactions in momentum space. We assume the appearance of magnetic structures characterized by Q 1 = (Q, 0), Q 2 = (0, Q), Q 3 = (Q, Q)/2, and 4 and Q 2 = Q 3 + Q 4 ; Q 1 and Q 2 are higher-harmonic wave vectors of Q 3 and Q 4 , as shown in figure 1(a). We also assume that the instabilities at Q 1 and Q 2 are larger than those at Q 3 and Q 4 at zero field. We neglect the contributions from the other wave vectors including other higher harmonics of Q 1 -Q 4 . Such a simplification can be justified when the interactions at Q 1 and Q 2 are larger than those at other wave vectors. The situation is achieved by considering short-ranged competing exchange interactions and long-ranged exchange interactions, the latter of which often arises from the Fermi surface instability in itinerant electron systems [56,57]. In the end, a low-temperature phase diagram is effectively determined by a few dominant Q ν -channel interactions in equation (1). Hereafter, we set J Q 1 = 1 as the energy unit of the model, whose energy scale roughly corresponds to the transition temperature. When the tetragonal symmetry is supposed, the relation of In addition, we consider the anisotropic form factor of Γ αβ Q ν in spin space, which microscopically arises from the interplay among relativistic spin-orbit coupling, crystalline electric field, and electronic band structure [58]. Based on tetragonal symmetry, we set Γ yy 6 ; the other components are zero. It is noted that the anisotropic form factor is characterized by the symmetric tensor satisfying which corresponds to the Dzyaloshinskii-Moriya interaction, appears due to the presence of the spatial inversion symmetry. We choose the model parameters as (γ 1 , γ 2 , γ 3 , γ 4 , γ 5 , γ 6 ) = (0.85, 0.95γ 1 , 1, 0.9, 0.025, 1) and J Q 3 /J Q 1 = J Q 4 /J Q 1 = 0.95 to reproduce the low-temperature phase diagram in EuAl 4 when the tetragonal symmetry is supposed, i.e. J Q 1 = J Q 2 [52]. Especially, the latter parameter J Q 3 /J Q 1 = J Q 4 /J Q 1 ∼ 1 indicates that the effect of magnetic frustration in momentum space is large, which becomes the origin of two types of the SkXs (S-SkX and R-SkX) [59,60]. Indeed, for this set of model parameters, it was revealed that the above two SkXs appear in the presence of an external magnetic field H in the second term in equation (1) [52].
In the following, we further investigate the stability of the SkXs against the uniaxial distortion under the [100] direction, where the tetragonal crystal symmetry is lowered to the orthorhombic one. Reflecting the inequivalence between the [100] and [010] directions, We construct the phase diagram while changing J Q 2 /J Q 1 . For simplicity, we fix the position of the ordering vectors and the other parameters against the distortion.
The low-temperature phase diagram of the model with the system size of N = 100 2 spins is investigated by simulated annealing based on the standard Metropolis local updates. Starting from a random spin configuration at high temperatures T/J Q 1 > 1, we gradually reduce the temperature following the relation T n+1 = α A T n until the temperature is reached to the final temperature T = 0.01 so that the effect of thermal fluctuations are almost negligible; T n is the nth-step temperature and α A = 0.999 95-0.999 99. At the final temperature, we perform a total of 10 5 -10 6 Monte Carlo sweeps to calculate physical quantities.

Results
is the unit vector in the x (y) direction]. Among the obtained eight phases, only the S-SkX and the R-SkX have nonzero χ 2 sc ; both states show an integer skyrmion number of −1, whose sign is determined by the parameters γ 1 , γ 2 , and γ 5 [61,62]. To clearly represent the similarity and difference of the magnetic textures in each phase, we show the contours of magnetic moments for the Q 1 component in figure 2(a), the Q 2 component in figure 2(b), and the Q 3 component in figure 2(c); the q component of the magnetic moments is given by (m α q ) 2 = (1/N 2 ) i,j S α i S α j e iq·(r i −r j ) for α = x, y, z, where r i is the position vector at spin i. We also define (m q ) 2 = (m x q ) 2 + (m y q ) 2 + (m z q ) 2 . It is noted that (m Q 3 ) 2 = (m Q 4 ) 2 in the obtained phases.
For J Q 2 /J Q 1 = 1, five phases appear against H: the 1Q proper-screw spiral, 1Q ′ , R-SkX, S-SkX, and 2Q states. In the low-field region, the 1Q proper-screw spiral state is characterized by m y Q 1 and m z Q 1 or m x Q 2 and m z Q 2 , while the 1Q ′ state has additional small intensities at Q 3 and Q 4 , as shown in figures 2(a)-(c). As shown In the R-SkX, the intensity at Q 1 is different from that at Q 2 even for J Q 2 /J Q 1 = 1, which indicates spontaneous fourfold rotational symmetry breaking like the low-field 1 Q states. Meanwhile, the S-SkX exhibits equivalent intensities at Q 1 and Q 2 so as to preserve fourfold rotational symmetry. Such a difference in terms of rotational symmetry is clearly found in the real-space spin configurations in figures 3(c) and (d). In the high-field region, the 2Q state appears, whose structure is mainly characterized by the dominant intensities of m x When setting J Q 2 /J Q 1 < 1 under the orthorhombic distortion, the S-SkX is suddenly replaced by the R-SkX owing to the breaking of fourfold rotational symmetry of the lattice structure; m Q 1 ̸ = m Q 2 . While decreasing J Q 2 /J Q 1 , one finds that there are no new phases in the phase diagram. The phase boundary between the 1Q proper-screw spiral and 1Q ′ states remains the same against J Q 2 /J Q 1 , since there is no intensity at Q 2 in these states. Similarly, the phase boundary between the 1Q ′ and the R-SkX is almost the same against J Q 2 /J Q 1 except for J Q 2 /J Q 1 ≳ 0.96. This indicates that the R-SkX in the low-field region is almost characterized by the anisotropic triple-Q structures at Q 1 , Q 3 , and Q 4 for J Q 2 /J Q 1 ≲ 0.96. Indeed, almost no (m Q 2 ) 2 is found in the R-SkX phase for small H in figure 2(b). The phase boundary between the R-SkX and 2Q states moves downward with decreasing J Q 2 /J Q 1 , since there is a contribution of m Q 2 in the high-field region of the R-SkX and its energy gain in the R-SkX becomes small for small J Q 2 /J Q 1 .
For J Q 2 /J Q 1 > 1, the S-SkX is also replaced by the R-SkX. In contrast to the situation with J Q 2 /J Q 1 < 1, the high-field 2Q state turns into the 2Q ′ , 1Q conical spiral, and 1Q fan states depending on J Q 2 /J Q 1 . These states are mainly characterized by the 1Q spin modulation at Q 2 . The 2Q ′ and 1Q conical spiral states have the intensities of m x Q 2 and m y Q 2 , while the 1Q fan state has the intensity of m x Q 2 . The difference between the 2Q ′ and 1Q conical spiral states is found in the additional modulations from the magnetic moments at Q 1 , Q 3 , and Q 4 ; the former has nonzero contributions of (m Q 1 ) 2 , (m Q 3 ) 2 , and (m Q 4 ) 2 , while the latter does not (see also figures 2(a)-(c)). The spin configuration of the 2Q ′ state is shown in figure 3(f). The R-SkX is no longer From the result of the phase diagram, one finds the essence to stabilize the R-SkX in the low-symmetric lattice structure. The important observation is that the R-SkX is robustly realized for J Q 2 /J Q 1 < 1. Especially, the R-SkX remains stable even for 1 , which indicates that m Q 2 plays a less role in the stabilization of the R-SkX; m Q 2 becomes negligibly small for J Q 2 /J Q 1 < 1 ( figure 2(b)), as discussed above. Thus, the multiple-Q feature in the R-SkX changes from the 2Q (or four-Q) structure with the intensities at Q 1 -Q 4 to the triple-Q structure with those at Q 1 , Q 3 , and Q 4 . The formation of the R-SkX with the triple-Q modulation is reasonable, since the contribution to the free energy in the coupling form of (S 0 · S Q 1 )(S −Q 3 · S Q 4 ) occurs like the triangular SkX in hexagonal systems owing to the relation of the triple-Q ordering vectors as Q 1 + (−Q 3 ) + Q 4 = 0. Indeed, it was shown that a similar R-SkX remains stable on a uniaxially-distorted triangular lattice [47]. The above tendency holds for J Q 2 /J Q 1 > 1. In this case, the triple-Q ordering vectors are composed of Q 2 , Q 3 , and Q 4 . The disappearance of the R-SkX for large J Q 2 /J Q 1 indicates that the ratio of the interaction channels between Q 2 and Q 3 is an important factor; the R-SkX is realized when the interactions at the triple-Q ordering vectors take relatively close values to each other, i.e. for J Q 3 /J Q 2 ≳ 0.812.
Although the S-SkX abruptly changes to the R-SkX under the orthorhombic lattice distortion, the fourfold symmetry tends to recover against the magnetic field in the R-SkX region. To quantify such a tendency, we evaluate the following quantity as where r squ = 0 means the presence of fourfold rotational symmetry. We show the contour plot of r squ in figure 2(d), where we only calculate r squ in the S-SkX and R-SkX regions. We also present the H dependence of r squ for several J Q 2 /J Q 1 in figure 4. The result clearly shows that the fourfold anisotropy is stronger for larger H. Thus, the shape of each skyrmion in real space looks like a square (rectangular) in the high(low)-field region. Finally, let us compare the present result with the experimental phase diagram in EuAl 4 [49][50][51][52][53], where the recent experiment suggested the importance of a spin-lattice coupling through a tetragonal-to-orthorhombic structural distortion [54,55]. Our present results indicate that the orthorhombic distortion, which makes the [100] and [010] directions inequivalent, gives a qualitatively similar result without distortion. Thus, the emergence of the R-SkX in EuAl 4 can be attributed to the synergy between the momentum-space frustrated interactions and orthorhombic distortion. In addition, we find a tendency that fourfold rotational symmetry recovers against the magnetic field, which means that the instability toward the S-SkX becomes larger for a larger magnetic field, although we could not obtain the S-SkX with r squ = 0 under the orthorhombic distortion. As shown in figure 4, r squ continuously changes as 0 < r squ < 1 for J Q 2 ̸ = J Q 1 . Accordingly, the uniform magnetization M z = (1/N) i S z i is also continuous, as shown in figure 4, which is different from the experimental situation where a discontinuous change of M z was found in the transition from the R-SkX to the S-SkX [55]. Thus, one needs to additionally consider the effects, such as the tilted ordering vector toward the [110] direction, in order to reproduce the transition between the R-SkX and S-SkX under the orthorhombic distortion. In addition, the effect of thermal fluctuations on the obtained phases is also an intriguing issue. Such additional effects on the phase diagram would be left as issues for future study.

Summary
To summarize, we have investigated the stability of the S-SkX and R-SkX appearing in the centrosymmetric tetragonal system against the orthorhombic distortion. By performing simulated annealing, we have shown that the R-SkX is robustly stabilized under distortion, while the S-SkX is suddenly replaced by the R-SkX. We found that twofold anisotropy in the R-SkX becomes weak in the high-field region, which indicates that the instability toward the S-SkX becomes large. Moreover, our result supported the appearance of the low-symmetric SkX under the orthorhombic distortion as observed in EuAl 4 , which will stimulate further exploration of the SkX in the low-symmetric lattice structures without a threefold or fourfold rotational axis.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.