Effect of In-Plane Magnetic Field on Skyrmions in a Centrosymmetric Triangular-Lattice System with Symmetric Anisotropic Exchange Interaction

: We report our numerical results on the stability of the skyrmion crystal phase in an ex-ternal magnetic field for both in-plane and out-of-plane directions in a centrosymmetric host. We analyze a spin model with the two-spin symmetric anisotropic exchange interaction that arises from relativistic spin–orbit coupling on a triangular lattice. By performing simulated annealing, we construct magnetic phase diagrams when the magnetic field is tilted from the out-of-plane field direction to the in-plane field direction. We find a different stability tendency of the skyrmion crystal phase according to the directions of the in-plane field, which provides a signal of the two-spin symmetric anisotropic exchange interaction for stabilizing the skyrmion crystal phase. Our results indicate that the mechanism of the skyrmion crystal phase triggered by the two-spin symmetric anisotropic exchange interaction can be experimentally tested by applying the in-plane magnetic field.

In contrast to noncentrosymmetric itinerant magnets, where the DM interaction plays an important role in stabilizing the SkX [97,98], effective multiple-spin interactions in the form of (S i • S j ) 2 and (S i • S j )(S k • S l ) [30,35], easy-axis single-ion anisotropy in the form of (S z i ) 2 [31,99], and two-spin symmetric anisotropic exchange interaction in the form of S x i S x j with the directional dependence [49,100] become the origin of the SkX in centrosymmetric itinerant magnets, where S i = (S x i , S y i , S z i ) is the localized spin at site i.Meanwhile, it is usually difficult to directly identify which factors are essential in SkX-hosting materials.
In the present study, we investigate different tendencies among the above mechanisms by examining the stability of the SkX under in-plane and out-of-plane external magnetic fields.Since the symmetry of the aforementioned spin interactions differs, one can obtain different stability tendencies depending on the magnetic field directions.In order to demonstrate this, we construct two magnetic phase diagrams for an effective spin model with the two-spin symmetric anisotropic exchange interaction on a two-dimensional triangular lattice: one corresponding to the field direction lying in the xz-plane and the other corresponding to the field direction lying in the yz-plane.We show that the stability region and behavior of the triple-Q SkX are different depending on the magnetic field direction.We also show that the instability toward multiple-Q states, other than the SkX, is qualitatively different for different magnetic field directions.We discuss the differences in the phase diagrams from other mechanisms, such as multiple-spin interactions and easy-axis single-ion anisotropy.Our results would be helpful for clarifying whether the two-spin symmetric anisotropic exchange interaction plays a significant role in stabilizing the SkX in centrosymmetric itinerant magnets.
The rest of this paper is organized as follows.In Section 2, we introduce an effective spin model in centrosymmetric itinerant magnets, where the effect of the two-spin symmetric anisotropic exchange interaction is incorporated.We also outline the numerical method based on simulated annealing.In Section 3, we discuss the numerical results by focusing on the stability of the SkX in an external magnetic field.First, we show that the SkX is induced by the two-spin symmetric anisotropic exchange interaction in an out-of-plane magnetic field.Then, we discuss the stability of the SkX by rotating the field direction in the xzand yz-planes.We compare the obtained results with those from models incorporating multiple-spin interactions and easy-axis single-ion anisotropy in Section 4. Section 5 is devoted to a summary of this paper.

Model and Method
We consider the following spin model on a two-dimensional triangular lattice, which is given by where spin moments, which is related to S i by Fourier transformation; and ν is the index for independent wave vectors in the first Brillouin zone.The first term represents the interaction defined in momentum space with the coupling constant J, which consists of the two-spin isotropic exchange interaction and symmetric anisotropic exchange interaction.I αβ Q ν for α, β = x, y, z represents the Q νdependent anisotropic form factor, which originates from the interplay between spin-orbit coupling and the crystalline electric field in the triangular-lattice structure, as demonstrated in Ref. [101].The momentum-resolved bilinear exchange interaction can be mapped into the real-space one through Fourier transformation.The second term represents the Zeeman coupling term through an external magnetic field H = (H x , H y , H z ).
The spin model with the momentum-resolved interaction in Equation ( 1) is derived from the Kondo lattice model with classical spin by considering a situation where the Kondo coupling between itinerant electron spins and localized spins (J K ) is small compared to the bandwidth of itinerant electrons [30,101,102].By performing perturbative analysis in terms of the Kondo coupling, one finds that the interaction in the first term in Equation (1) corresponds to the lowest-order contribution in terms of the Kondo coupling in the expansion, i.e., J ∝ J 2 K , which is known as the Ruderman-Kittel-Kasuya-Yosida interaction [103][104][105], whose magnitude is determined by the bare susceptibility of itinerant electrons depending on the band structure and Fermi level in the system.We ignore the effect of multiple-spin interactions, which appear as higher-order contributions in the expansion, by assuming that the higher-order effect is negligible, although a situation where the bare susceptibility shows distinct peak structures at Q ν often leads to instability toward the SkX [30].
In Equation (1), we further simplify the model by taking into account the momentumresolved interactions at particular wave vectors that provide the dominant contribution to the ground-state energy.In other words, we consider the situation where the bare susceptibility of itinerant electrons shows maxima at particular wave vectors Q 1 -Q 3 , which are often attained by the nesting of the Fermi surface.Specifically, we set ), and , where the lattice constant of the triangular lattice is taken as unity.It is noted that Q 1 , Q 2 , and Q 3 are connected by the threefold rotational symmetry of the triangular lattice.The following results are not qualitatively affected by a different choice of Q unless the ordering wave vectors lie at the Brillouin zone boundary.
For these wave vectors, I αβ Q ν is given so as to satisfy the sixfold rotational symmetry as follows [101]: This symmetric anisotropic form factor I A corresponds to bond-dependent interaction, such as the compass and Kitaev interactions [106][107][108][109][110], which become the origin of unconventional topological spin textures [111][112][113][114] and nonreciprocal magnon excitations [110].Furthermore, this type of symmetric anisotropic exchange interaction on the triangular lattice induces the SkX even in the absence of multiple-spin interactions or easy-axis single-ion anisotropy [49,100].We focus on the behavior of the SkX induced by I A against the magnetic field direction in Section 3 and compare the results among different mechanisms in Section 4.
For the model in Equation ( 1), we construct two low-temperature magnetic phase diagrams depending on the magnetic field direction: one is the phase diagram obtained by varying the out-of-plane field H z and the in-plane x-directional field H x with H y = 0, and the other is the phase diagram obtained by varying H z and the y-directional field H y with H x = 0. We set the other model parameters as J = 1 and I A = 0.1.The phase diagrams are calculated by performing simulated annealing based on the standard Metropolis algorithm with local spin updates for the system size N = 24 2 under periodic boundary conditions.The procedure is as follows.By taking a random spin configuration as an initial spin state, we gradually reduce the temperature from high temperatures, T 0 = 1-5, to the final temperature, T = 0.0001, using the ratio T n+1 = 0.999999T n , where T n is the temperature at the nth-step.We perform the local spin updates in real space at each temperature up to T. At T, we further perform Monte Carlo sweeps around 10 5 -10 6 in order to optimize the spin configuration.We independently perform the above procedure for different sets of magnetic fields.We also perform the simulations from the spin configurations obtained at low temperatures in the vicinity of the phase boundaries in order to avoid meta-stable spin configurations.
For the spin configurations obtained through simulated annealing, we calculate the spin and scalar chirality quantities to identify the magnetic phase.The order parameter corresponds to the Q ν component of magnetic moments m η Q ν with the spin component η = x, y, z, which is given by where r i is the position vector at site i.We also calculate (m The squared quantity (m Q ν ) 2 corresponds to the susceptibility in the Q ν channel.The uniform magnetization is given by The scalar chirality is given by where R denotes the position vectors at the centers of triangles and µ = (u, d) represents upward and downward triangles, respectively, in the triangular lattice.The site indices i, j, k represent the vertices on the triangle at R in counterclockwise order.The magnetic states with nonzero χ sc correspond to the topologically nontrivial states since they generate an auxiliary magnetic field and act on the itinerant electrons via the spin Berry phase mechanism.

Results
We show the stability of the SkX and other multiple-Q states in the different magnetic field directions.In Section 3.1, we show that the SkX appears in the ground state under the out-of-plane magnetic field with the aid of the two-spin symmetric anisotropic exchange interaction.Then, we show the effect of the in-plane magnetic field on the SkX.The results under the x-directional magnetic field are shown in Section 3.2, whereas those under the y-directional magnetic field are shown in Section 3.3.

Out-of-Plane Field
First, we investigate the instability toward the SkX in the model in Equation ( 1) under the out-of-plane magnetic field, i.e., H x = H y = 0 and H z ̸ = 0.In the case of I A = 0, the system exhibits the single-Q conical spiral state under H z , whose spin configuration is given by S i = (sin θ cos Q ν • r i , sin θ sin Q ν • r i , cos θ) with cos θ = H z /2.Thus, there is no instability toward multiple-Q states in the absence of I A .
Meanwhile, the multiple-Q states, including the SkX, appear by introducing I A [100]. Figure 1 shows the H z dependence of the magnetic and scalar chirality quantities at I A = 0.1.For better readability, we sort (m 1b and (m z Q ν ) 2 in Figure 1c when the states with different Q ν are energetically degenerate.There are six phases except for the fully polarized state stabilized for H z ≳ 2.2; the spin configuration of the fully polarized state is given by S i = (0, 0, 1).We describe the details of the obtained phases one by one below.The nonzero scalar chirality and magnetic moments in each phase are summarized in Table 1.
Phase I.This state is stabilized at zero and low fields.At zero fields, this state exhibits an anisotropic double-Q modulation; the dominant contribution is given by (m and the subdominant contribution is given by (m xy Q 2 ) 2 , as shown in Figure 1b,c.Specifically, the spin configuration is represented by a superposition of the proper-screw spiral wave in the Q 1 component and the sinusoidal wave in the Q 2 component.Such a feature can be seen in the real-space spin configuration, as shown in Figure 2a, where it is noted that the dominant and subdominant components are given by (m Q 3 ) 2 and (m Q 1 ) 2 , respectively.Owing to the double-Q superposition, this spin configuration is noncoplanar, and the density wave in terms of the scalar chirality occurs in the subdominant (m Q ν ) 2 component, as shown in Figure 3a.On the other hand, there is no uniform scalar chirality in this state, as shown in Figure 1d.
When H z is applied, the magnetization M z continuously increases, as shown in Figure 1a.Accordingly, (m xy Q 3 ) 2 becomes nonzero, as shown in Figure 1b, which indicates that Phase I is characterized by an anisotropic triple-Q state.No net scalar chirality is induced in this triple-Q spin texture, as shown in Figure 1d.
Phase II.This phase appears next to Phase I by increasing H z , as shown in Figure 1a.The phase transition is of second order, where the magnetization smoothly changes, as shown in Figure 1a.The spin configuration in this phase is similar to that in Phase I, as shown in Figure 2a,b.Meanwhile, there is a small difference in (m xy Q ν ) 2 for ν = 2, 3; in contrast to Phase I, the subdominant intensity in (m Reflecting this difference, the scalar chirality in real space exhibits a checkerboard-type modulation, as shown in Figure 3b, which indicates that the scalar chirality is characterized by the double-Q modulation.Similar to Phase I, this state does not exhibit a net scalar chirality.
SkX II.When the magnetic field is increased in Phase II, the SkX II appears with a jump in the magnetization, as shown in Figure 1a, which indicates a first-order phase transition.This state is described by a superposition of the triple-Q spiral waves in the Q 1 , Q 2 , and Q 3 components.As shown in Figure 1b,c, the triple-Q superposition is anisotropic; (m Q 1 ) 2 is equal to (m Q 2 ) 2 , but neither is equal to (m Q 3 ) 2 .This feature appears in the real-space spin configuration in Figure 2c, where the shape of the skyrmion core, denoted as S z i = −1, seems elliptic rather than circular; the skyrmion core is located on the interstitial site owing to the discrete lattice structure without easy-axis single-ion anisotropy.The deviation from the circular shape of the skyrmion core is caused by the breaking of the threefold rotational symmetry, which is found in the relation (m The scalar chirality in real space is also characterized by an anisotropic triple-Q structure, as shown in Figure 3c.In contrast to Phase I and Phase II, this state exhibits a nonzero scalar chirality in Figure 1d, which signals the topologically nontrivial spin texture leading to the topological Hall effect.
SkX.With a further increase in H z , the SkX II turns into the SkX, as shown in Figure 1a.This state is also characterized by the triple-Q spiral states, although the intensities in the triple-Q components are equivalent to each other, i.e., (m 2 in contrast to the situation in the SkX II (see Figure 1b,c).Owing to the isotropic triple-Q structure, the real-space spin configuration keeps the threefold rotational symmetry of the triangular lattice, as shown in Figure 2d.Accordingly, the scalar chirality is also distributed in a threefold-symmetric way, as shown in Figure 3d.The net scalar chirality becomes nonzero in Figure 1d.
The appearance of both the SkX and SkX II is owing to the two-spin symmetric anisotropic exchange interaction I A .Since I A > 0 (I A < 0) favors the proper-screw (out-ofplane cycloidal) spiral wave, the resultant triple-Q SkX favors the Bloch (Néel) SkX.The Bloch helicity means that the spins surrounding the skyrmion core rotate in the xy-plane with a helicity of ±π/2, whereas the Néel helicity means that the spins are radially directed in an inward or outward way with a helicity of 0 or π.Indeed, we obtain the Bloch-type helicity in both the SkX and SkX II phases for I A = 0.1 > 0, as shown in Figure 2c,d.It is noted that the introduction of I A lifts the degeneracy with the anti-type SkX with the opposite sign of the scalar chirality.A similar situation also occurs in the square SkX in the tetragonal lattice system [115].
Phase III.This phase appears upon increasing H z in the SkX with a jump in the magnetization, as shown in Figure 1a.The phase transition between Phase III and SkX is of first order.The spin configuration in this state is characterized by the sinusoidal waves in two out of three (m xy Q ν ) 2 , as shown in Figure 1b, which means the appearance of the double-Q fan state.On the other hand, there are no z-spin modulations in the Q ν component, as shown in Figure 1c.The in-plane double-Q spin texture leads to the periodic alignment of the vortex and anti-vortex, as shown by the real-space spin configuration in Figure 2e.The distribution of the scalar chirality is described by a superposition of the Q 1 − Q 2 and Q 3 components when the magnetic modulations are represented by the Q 1 and Q 2 components, as shown in Figure 3e.Owing to the equivalence of the vortex and anti-vortex, the net scalar chirality becomes zero in the whole system, as shown in Figure 1d.One finds that the scalar chirality shows a jump from the finite values to zero in the transition from the SkX, as shown in Figure 1d.
Phase IV.This phase is stabilized in the high-field region close to the saturation field, as shown in Figure 1a.As shown in Figure 1b,c, the spin configuration in this state is described by an in-plane triple-Q fan structure with nonzero (m The real-space spin and scalar chirality configurations are shown in Figures 2f and 3f, respectively, where the threefold symmetric distributions are seen in both configurations.In the real-space picture, one finds two types of vortices: one is the vortex with a winding number of −2, whereas the other is the vortex with a winding number of +1.Since the number of the vortex with the winding number of +1 is twice that of the one with the winding number of −2, the net scalar chirality becomes zero in the whole system.A similar spin texture is also found in the frustrated triangular-lattice magnets [116,117].This phase continuously turns into a fully polarized state with the increase in the magnetic field.

In-Plane Field along the x-Direction
We investigated the effect of the in-plane magnetic field on the SkX and other multiple-Q states obtained in Section 3.1.In this section, we examine the effect of the x-directional magnetic field H x while keeping H y = 0. Figure 4 shows the low-temperature phase diagram when H x and H z are varied.The contours in Figure 4a-d represent (χ sc ) 2 in Figure 4a, (m Q 1 ) 2 in Figure 4b, (m Q 2 ) 2 in Figure 4c, and (m Q 3 ) 2 in Figure 4d.
Since (χ sc ) 2 becomes nonzero only for the SkX and SkX II, this quantity signals the appearance of the skyrmion spin texture.As shown in Figure 4a, both SkX phases remain stable when the magnetic field is tilted from the out-of-plane direction to the in-plane x-direction up to around 54 • measured from the z-axis.Then, both SkX phases vanish when H x is further increased, which indicates that I A does not favor the SkX in an in-plane magnetic field along the x-direction.This is consistent with the results that the easy-plane single-ion anisotropy in the form of A(S z i ) 2 for A > 0 is required in order to stabilize the SkX in the in-plane magnetic field [118].
Since the magnetic field is tilted from the z-direction, the threefold rotational symmetry of the SkX is broken; the Q 1 component is no longer equivalent to the Q 2 and Q 3 components.Accordingly, the skyrmion core is elongated along the y-direction, as shown by the real-space spin configuration in Figure 5a.The scalar chirality distribution is also modulated so as to break the threefold rotational symmetry, although such a change is small, as shown in Figure 5b.Indeed, (χ sc ) 2 takes similar values in the entire region where the SkX and SkX II phases appear, as shown in Figure 4a.Meanwhile, the skyrmion size under the in-plane field is almost the same as that under the out-of-plane field since the present SkX consists of a superposition of spin-density waves at the same ordering wave vectors.We discuss the effect of the x-directional magnetic field on the other phases by examining the data at H x = 0.4 in Figure 6a-d and H x = 0.8 in Figure 6e-h.In the low-field region for H z ≲ 0.4 at H x = 0.4, the behavior of (m xy Q ν ) 2 is similar to that at H x = 0, as shown in Figures 1b and 6b.Meanwhile, (m z Q 2 ) 2 and (m z Q 3 ) 2 become nonzero for 0.1 ≲ H z ≲ 0.4 in Figure 6c, which is different from the result in Figure 1c.Thus, we denote this phase as Phase II' instead of Phase II.By increasing H x , (m Q 2 ) 2 and (m Q 3 ) 2 vanish while keeping nonzero (m Q 1 ) 2 for 0 ≲ H z ≲ 0.6, as shown in Figure 6f,g, which indicates that the magnetic structure is characterized by the single-Q spiral state; we denote it as Phase I'.The choice of the spiral waves in the Q 1 component is attributed to the energy gain by both I A and H x since the spiral plane lies on the yz-plane for the Q 1 component.
In the high-field region for 1.3 ≲ H z ≲ 2.2 at H x = 0.4, the magnitudes of the triple-Q ordering wave vectors become nonzero with different intensities in both the xy and z components, as shown in Figure 6b,c.We denote this phase as Phase III'.Such a feature holds for larger H x , as shown in the case of H x = 0.8 in Figure 6f,g.This anisotropic triple-Q phase is replaced with the single-Q state stabilized in the low-field region by further increasing H x , as shown in Figure 4b-d.The net scalar chirality remains zero in the regions of Phase I', Phase II', and Phase III', as shown in Figure 6d,h.

In-Plane Field along the y-Direction
Next, we discuss the results obtained by rotating the magnetic field on the yz-plane, i.e., H x = 0. Figure 7a-d show the phase diagram in the plane of H y and H z .Similar to the results in Figure 4, we show the contours of (χ sc ) 2 in Figure 7a, (m Q 1 ) 2 in Figure 7b, (m Q 2 ) 2 in Figure 7c, and (m Q 3 ) 2 in Figure 7d.
The SkX remains stable up to H y ∼ 0.9, as shown in Figure 7a.The stability region of the SkX under H y is similar to that under H x , although one can observe a slight difference between them; the SkX region for the former case is wider than that for the latter.The real-space spin and scalar chirality configurations in the presence of H y are shown in Figures 8a and 8b, respectively.Owing to the y-directional magnetic field, the spins tend to align the positive y-direction, which breaks the threefold rotational symmetry.In contrast to the result under H x , the SkX and SkX II are not distinguished from each other when H y is introduced, which means that they smoothly merge into a single phase, as shown in Figure 9b,c.
In the low-and high-field regions, the instabilities toward the magnetic states are qualitatively different from those in Section 3.2.In the low-field region, the dominant intensity of (m Q ν ) 2 is found in the Q 2 (or Q 3 ) component in order to gain the Zeeman energy with H y , as shown in Figure 9b,c,f,g.We denote this phase as Phase I".In addition, the magnetic state with dominant double-Q intensities (m Q 2 ) 2 = (m Q 3 ) 2 appears in the high-field region, where we denote it as Phase III", as shown in Figure 9b,c,f,g.It is noted that (m Q 1 ) 2 is slightly induced in this state.This phase remains stable with the increase in H y ; the double-Q feature can be seen in the vicinity of H z = 0 and H y = 2, as shown in Figure 7c,d.Thus, one can find that the importance of I A can be tested by examining the magnetic phases against H x and H y for H z = 0.When I A plays an important role, the magnetic phase corresponds to the single-Q spiral state for H x ̸ = 0, whereas it corresponds to the triple-Q state with strong double-Q spin modulations for H y ̸ = 0.

Comparison with Other Mechanisms
Finally, let us briefly compare the present results with other results obtained under different interactions and magnetic anisotropy.In the present model, we consider the situation where the SkX is induced by the two-spin symmetric anisotropic exchange interaction I A .Meanwhile, there are other mechanisms to induce the SkX in centrosymmetric itinerant magnets, such as the positive biquadratic interaction in the form of (S q • S −q ) 2 [30] and the easy-axis single-ion anisotropy in the form of −A(S z i ) 2 [31,100].Although the phase diagrams under the out-of-plane magnetic field are similar to each other, those under the in-plane field are expected to be different, which is easily understood from the symmetry aspect.
In the case of the positive biquadratic interaction, the effect of the x-directional field is the same as that of the y-directional field since this isotropic interaction is invariant under spin rotation.Similarly, the easy-axis single-ion anisotropy also does not bring about a difference in terms of H x and H y owing to the spin rotational symmetry around the z-axis.Thus, the mechanism based on the two-spin symmetric anisotropic exchange interaction can be distinguished from the other two mechanisms by investigating the magnetic states under the in-plane magnetic fields along the two directions.

Summary
We have numerically investigated the effect of the in-plane magnetic field on the SkX in centrosymmetric itinerant magnets with the two-spin symmetric anisotropic exchange interaction.By performing simulated annealing for the effective spin model with the momentum-resolved interaction on the triangular lattice and constructing two magnetic phase diagrams in the in-plane xand y-directional magnetic fields, we show the similarity and difference of the SkX stability according to the field direction.Furthermore, we show that the different instabilities toward the magnetic phases are caused by the different inplane field directions, which enable us to distinguish the microscopic mechanism of the SkX induced by the two-spin symmetric anisotropic exchange interaction from the others by the multiple-spin interaction and single-ion anisotropy.One of the candidate materials is the breathing kagome compound Gd 3 Ru 4 Al 12 [50,51,[119][120][121], where the SkX has been identified in experiments [48], and its emergence has been accounted for by the two-spin symmetric anisotropic exchange interaction [49].

Figure 1 .
Figure 1.Out-of-plane magnetic field H z dependence of (a) the magnetization M η , (b) the Q ν component of the squared in-plane magnetic moment (m xy Q ν ) 2 , (c) the Q ν component of the squared out-of-plane magnetic moment (m z Q ν ) 2 , and (d) the squared scalar chirality (χ sc ) 2 for I A = 0.1 and H x = H y = 0. FP stands for the fully polarized state.

Figure 2 .
Real-space spin configurations of (a) Phase I at H z = 0, (b) Phase II at H z = 0.4, (c) SkX II at H z = 0.5, (d) SkX at H z = 0.8, (e) Phase III at H z = 1.35, and (f) Phase IV at H z = 1.9.The contour shows the z component of the spin moment, and the arrows represent the xy components.

Table 1 .
Scalar chirality χ sc and Q ν components of the magnetic moments m Q ν for ν = 1-3 in each phase.

Figure 5 .
Figure 5. (a) Real-space spin configurations of the SkX at H x = 0.75 and H z = 0.8.(b) Real-space scalar chirality configurations corresponding to (a).

Figure 8 .Figure 9 .
Figure 8.(a) Real-space spin configurations of the SkX at H y = 0.8 and H z = 0.8.(b) Real-space scalar chirality configurations corresponding to (a).