Double-Q Checkerboard Bubble Crystal in Centrosymmetric Tetragonal Magnets

Abstract

We report our numerical studies on the emergence of a double-Q checkerboard bubble crystal in centrosymmetric tetragonal magnets. The double-Q checkerboard bubble crystal is characterized by a fourfold-symmetric collinear spin configuration consisting of a superposition of two sinusoidal waves with the out-of-plane spin modulations along the [110] and [$\overline{1}$10] directions. The numerical calculations based on the simulated annealing for an effective spin model with the momentum-resolved easy-axis exchange interactions reveal that the double-Q checkerboard bubble crystal is energetically degenerate with the single-Q collinear state when the ordering wave vector lies on the quarter of the reciprocal lattice vector along the $⟨110⟩$ direction. We show that such a degeneracy is lifted by considering the biquadratic interaction. We also find that the double-Q checkerboard bubble crystal turns into another double-Q state characterized by the in-plane spin modulations by increasing an external magnetic field.

1. Introduction

Exploring unconventional magnetic orderings has been one of the central issues in condensed matter physics. Depending on the spin configurations, various intriguing physical phenomena occur, such as the anomalous Hall effect in ferromagnets [1,2,3,4,5,6,7,8,9,10] and collinear antiferromagnets [11,12,13,14,15,16,17,18,19,20], spin current generation in noncollinear antiferromagnets [21,22,23,24,25,26,27,28], and the nonreciprocal transport in noncoplanar antiferromagnets [29,30,31,32,33,34,35,36]. The findings of unconventional magnetic orderings and their related physical phenomena can not only provide fundamental understandings but also bring about future electronic and spintronic device applications [37,38,39].

The magnetic orderings can be described by a superposition of spin density waves, where the spin at site i, ${\mathit{S}}_i$, is generally expressed as

$$\begin{array}{c}\hfill {\mathit{S}}_i=\sum _{\nu =1}^n(e^{i{\mathit{Q}}_{\nu }·{\mathit{r}}_i}{\mathit{S}}_{{\mathit{Q}}_{\nu }}+e^{-i{\mathit{Q}}_{\nu }·{\mathit{r}}_i}{\mathit{S}}_{-{\mathit{Q}}_{\nu }}).\end{array}$$

(1)

Here, ${\mathit{S}}_{\mathit{q}}$ with the wave vector $\mathit{q}$ is the Fourier transform of ${\mathit{S}}_i$, and ${\mathit{r}}_i$ is the position vector at site i. For $n=1$, the spin state is referred to as the single-Q state, while it is referred to as the multiple-Q state for $n>1$ [40,41,42,43]. The multiple-Q spin configurations can describe a variety of complicated spin configurations, such as a skyrmion crystal [44,45,46,47,48,49,50,51], hedgehog crystal [52,53,54,55,56,57], meron–antimeron crystal [58,59,60,61,62], and vortex crystal [63,64,65,66,67]. For example, a double-Q superposition of spiral waves leads to the skyrmion crystal on a two-dimensional square lattice and a triple-Q superposition of spiral waves leads to the hedgehog crystal on a three-dimensional cubic lattice. Since the resultant multiple-Q spin configurations depend on the type of the constituent spin density waves, further exotic magnetic orderings can be expected.

In the present study, we explore a further multiple-Q state and its stabilization mechanism based on the spin model analysis. We perform the simulated annealing for the spin model on the two-dimensional centrosymmetric square lattice and construct the magnetic phase diagram at low temperatures. By changing an easy-axis anisotropic interaction and an external magnetic field, we find that the energy of the double-Q state is the same as that of the single-Q state within the bilinear exchange interaction. The infinitesimally small biquadratic interaction stabilizes two types of the double-Q states depending on the magnetic field. In the small field region, a double-Q checkerboard bubble crystal (2Q CBX), which consists of double-Q sinusoidal waves with the z-spin modulation, appears without the uniform magnetization, while a double-Q state with $xy$-spin modulations (2Q I state) appears in the high-field region. Our results indicate a possibility of further multiple-Q states by considering the interplay between the biquadratic interaction and magnetic field.

The rest of this paper is organized as follows: In Section 2, we construct an effective spin model on the square lattice to induce the 2Q CBX by taking into account the momentum-resolved spin interactions. We also introduce the simulated annealing, which enables us to obtain the thermodynamic spin configurations in an unbiased way. In Section 3, we discuss the stability of the 2Q CBX by changing the model parameters. First, we show that the energy of the 2Q CBX is degenerate with that of the 1Q collinear state within the easy-axis anisotropic bilinear exchange interaction. Then, we show that the biquadratic interaction favors the 2Q CBX rather than the 1Q collinear state. Finally, we present the magnetic phase diagram against the external magnetic field. Section 4 is devoted to the conclusion of the present paper.

2. Model and Method

We analyze the following spin model with the momentum-resolved interaction on the two-dimensional square lattice, whose Hamiltonian is given by

$$\begin{array}{cc}\hfill \mathcal{H}=& -2\sum _{\nu }\left(J{\mathsf{\Gamma }}_{{\mathit{Q}}_{\nu }}+{\displaystyle \frac{K}{N}}{\mathsf{\Gamma }}_{{\mathit{Q}}_{\nu }}^2\right)-H\sum _i{S}_i^z,\hfill \end{array}$$

(2)

with

$$\begin{array}{c}\hfill {\mathsf{\Gamma }}_{{\mathit{Q}}_{\nu }}={\mathit{S}}_{{\mathit{Q}}_{\nu }}·{\mathit{S}}_{-{\mathit{Q}}_{\nu }}+I^z{S}_{{\mathit{Q}}_{\nu }}^z{S}_{-{\mathit{Q}}_{\nu }}^z,\end{array}$$

(3)

where ${\mathit{S}}_{{\mathit{Q}}_{\nu }}$ denotes the ${\mathit{Q}}_{\nu }$ component of the spin and $\nu$ is the index of the ordering wave vectors ${\mathit{Q}}_{\nu }$. ${\mathit{S}}_{{\mathit{Q}}_{\nu }}=(S_{{\mathit{Q}}_{\nu }}^x,S_{{\mathit{Q}}_{\nu }}^y,S_{{\mathit{Q}}_{\nu }}^z)$ is related to the localized spin at site i, ${\mathit{S}}_i=(S_i^x,S_i^y,S_i^z)$, via the Fourier transformation. We set the spin length as unity, i.e., $|{\mathit{S}}_i|=1$.

The first term represents the bilinear and biquadratic exchange interactions with the coupling constants J and K, respectively; N denotes the total number of spins in the system. We take into account the interactions at particular wave vectors by supposing the nesting of the Fermi surface at ${\mathit{Q}}_{\nu }$; we set ${\mathit{Q}}_1=(\pi /2,\pi /2)$ and ${\mathit{Q}}_2=(-\pi /2,\pi /2)$; and the lattice constant is set to be unity. It is noted that the coupling constants J and K are common for ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$, since they are connected by the fourfold rotational symmetry of the square lattice. The prefactor 2 is due to the contribution from $-{\mathit{Q}}_1$ and $-{\mathit{Q}}_2$. In the nesting regime, the bilinear exchange interaction corresponds to the Ruderman–Kittel–Kasuya–Yosida (RKKY) interaction [68,69,70], while the biquadratic exchange interaction corresponds to the higher-order RKKY interaction [71,72]. Microscopically, the RKKY (higher-order RKKY) interaction is obtained in the lowest-order (second-lowest-order) perturbative expansions with respect to the spin–charge coupling in the Kondo lattice model. The bilinear exchange interaction tends to favor the single-Q state at ${\mathit{Q}}_{\nu }$, while the biquadratic interaction tends to favor the multiple-Q state, such as the vortex crystal [73], bubble crystal [74], and the higher-order skyrmion crystal [75]. In addition, we introduce the dimension-less easy-axis anisotropic form factor in ${\mathsf{\Gamma }}_{{\mathit{Q}}_{\nu }}$ as $I^z>0$, while we neglect other symmetry-allowed magnetic anisotropic form factors, such as the bond-dependent magnetic anisotropy [76]. Such magnetic anisotropy in both two-spin and four-spin interactions originates from the interplay between the spin–orbit coupling and crystalline electric field; in particular, the easy-axis anisotropy arises from the inequivalence between the $xy$- and z-spin components. Although we choose the same ratio of the interactions between the $xy$ and z components for the bilinear and biquadratic interactions, the following results are not altered even for different interaction ratios in the $xy$ and z components of the biquadratic interaction, since the biquadratic interaction is introduced to lift the degeneracy between the single-Q and double-Q states detailed in Section 3 in the present study. The second term stands for the effect of the external magnetic field along the out-of-plane direction; H represents the magnitude of the magnetic field.

We examine the instability toward the multiple-Q states by numerically solving the model in Equation (2). Although the model includes only a few interactions in momentum space, it gives an almost accurate internal energy for the ground-state spin configurations when the interactions at the other wave vectors are smaller than those at ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$. In other words, the effective momentum-resolved spin model enables us to efficiently search for the multiple-Q states [77,78]. Such a simplification is especially useful for the examination of the instability toward the multiple-Q states in the ground state, since the ground-state energy is determined by the dominant multiple-Q (or single-Q) modulations, and the contribution from the higher-harmonic wave vectors like ${\mathit{Q}}_1+{\mathit{Q}}_2$ is negligible [79]. Indeed, the simplification of the model can capture the essential interactions in materials hosting multiple-Q magnetic orderings, such as EuPtSi [80], GdRu₂Si₂ [81], and EuNiGe₃ [82], where the experimental complicated phase diagrams are well reproduced. Meanwhile, the simplification of the model does not capture unconventional multiple-Q states, whose ordering wave vectors are not symmetry-related. In addition, the model with the momentum-resolved interaction is not appropriate for the analysis of the impurity effect. In the following, we set $J=1$ as the energy unit of the model, and change $I^z$ and H. For K, we set $K=0$ or $K=0.05$, since infinitesimally small K is enough to lift the degeneracy between the single-Q and multiple-Q states.

The magnetic phase diagram is calculated in the plane of $I^z$ and H based on the simulated annealing. For the system size with $N={16}^2$ under the periodic boundary conditions, we perform the Monte Carlo simulations with standard metropolis local updates in real space from a high temperature, $T_0=1.5$, to the final temperature, $T=0.0001$. In the simulations, we flip the spin in real space according to the Metropolis algorithm. Then, we evaluate the Fourier transform of the spins only for the ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$ components in order to evaluate the energy of the model in Equation (2). Starting from a random spin configuration at $T_0$, we gradually reduce the temperature with a rate of $T_{n+1}=0.999999T_n$ while performing the above procedure, where $T_n$ is the nth-step temperature. When the temperature reaches T, we perform ${10}^5$–${10}^6$ Monte Carlo sweeps for measurements. In addition, the simulations are performed from the spin configurations obtained at low temperatures in the vicinity of the phase boundaries in order to avoid the meta-stable solution. We independently perform the above procedure for a set of $I^z$ and H.

We calculate the following physical quantities in order to identify magnetic phases obtained by the simulated annealing. The spin structure factor is given by

$$\begin{array}{cc}\hfill S_s\left(\mathit{q}\right)& ={\displaystyle \frac{1}{N}}\sum _{i,j}({\mathit{S}}_i·{\mathit{S}}_j)e^{i\mathit{q}·({\mathit{r}}_i-{\mathit{r}}_j)},\hfill \end{array}$$

(4)

where $\mathit{q}$ is the wave vector in the first Brillouin zone. The net magnetization along the magnetic field direction is given by

$$M^z={\displaystyle \frac{1}{N}}\sum _i{S}_i^z.$$

(5)

3. Results

We discuss the multiple-Q instability in the model in Equation (2). First, we discuss the result without the biquadratic interaction in Section 3.1, where the single-Q and double-Q states are energetically degenerate with each other. Then, we introduce the biquadratic interaction to lift such a degeneracy in Section 3.2.

3.1. Without Biquadratic Interaction

Let us discuss the case without the biquadratic interaction, i.e., $K=0$. For $I^z=0$, the single-Q state with the up–up–down–down spin configuration is realized in the ground state. We deduce the spin ansatz satisfying the local spin length constraint for the spin configurations obtained by the simulated annealing. In the case of $H=0$, one of the expressions of this spin configuration is given by

$$\begin{array}{c}\hfill S_i^x=S_i^y={\displaystyle \frac{1}{\sqrt{2}}}(cos{\mathit{Q}}_1·{\mathit{r}}_i-sin{\mathit{Q}}_1·{\mathit{r}}_i),S_i^z=0,\end{array}$$

(6)

where we assume that the spin oscillates along the [110] direction. This spin configuration leads to the lowest energy of the model Hamiltonian, $-J=-1$, under the constraint of the spin length $|{\mathit{S}}_i|=1$.

Meanwhile, we find that the double-Q spin configuration gives the same lowest energy as the above single-Q state. By considering the superposition of two spin density waves at ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$ while the local spin length is kept, the spin ansatz for the double-Q state is given as follows:

$$\begin{array}{cc}\hfill S_i^x& =c_{1}cos{\mathit{Q}}_1·{\mathit{r}}_i-\sqrt{1-c_1^2}sin{\mathit{Q}}_2·{\mathit{r}}_i,\hfill \\ \hfill S_i^y& =-c_{1}sin{\mathit{Q}}_1·{\mathit{r}}_i+\sqrt{1-c_1^2}cos{\mathit{Q}}_2·{\mathit{r}}_i,\hfill \\ \hfill S_i^z& =0,\hfill \end{array}$$

(7)

where $0<c_1<1$ is the numerical coefficient. The same energy between the single-Q and double-Q states is unnderstood from the fact that the spin pattern in Equation (7) satisfies the local spin constraint as $|{\mathit{S}}_i|=1$, owing to the special value of the ordering wave vectors ${\mathit{Q}}_1=(\pi /2,\pi /2)$ and ${\mathit{Q}}_2=(-\pi /2,\pi /2)$. In other words, all the intensities in the spin structure factor lie at ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$ rather than the higher-harmonic wave vector like ${\mathit{Q}}_1+{\mathit{Q}}_2$ and other wave vectors. This indicates no energy loss in terms of the exchange energy compared to the single-Q state, which leads to the energy degeneracy between the single-Q and double-Q states. Such a situation also holds for nonzero H. A similar energy degeneracy has been found at particular ordering wave vectors like $\mathit{Q}=(\pi ,0)$ [72,83] and $\mathit{Q}=(5\pi /8,3\pi /8)$. Thus, the magnetic states with arbitrary $c_1$ are obtained after the simulated annealing. We show the real-space spin configurations in the left panels of Figure 1a,b. Although the spin configurations in Figure 1a,b seem to be completely different from each other, both are described by the expression in Equation (7) and their energy is given by $-J[1-{\left(M^z\right)}^2]+H{M}^z$, which takes the same value as the single-Q state. The corresponding spin structure factors are shown in the right panels of Figure 1a,b; all the intensities are found at ${\mathit{Q}}_1$, ${\mathit{Q}}_2$, and $\mathit{q}=\mathit{0}$.

When $I^z$ is introduced, the energy degeneracy between the single-Q state and the double-Q state also holds. At zero field, both single-Q and double-Q states are characterized by the collinear spin configurations with the z-spin component in order to gain the energy by $I^z$. The single-Q spin configuration is given by

$$\begin{array}{c}\hfill S_i^x=S_i^y=0,S_i^z=cos{\mathit{Q}}_1·{\mathit{r}}_i-sin{\mathit{Q}}_1·{\mathit{r}}_i,\end{array}$$

(8)

and the double-Q spin configuration is given by

$$\begin{array}{c}\hfill S_i^x=S_i^y=0,S_i^z=sin{\mathit{Q}}_1·{\mathit{r}}_i+cos{\mathit{Q}}_2·{\mathit{r}}_i.\end{array}$$

(9)

We show the real-space spin configuration of the double-Q state in the left panel of Figure 1c and that of the single-Q state in the left panel of Figure 1d. Their spin structure factors are shown in the right panels of Figure 1c,d. For the double-Q state, the up and down spins form the $2\times 2$ squares, and they are aligned in a way that resembles a checkerboard. We refer to this state as the 2Q CBX. A similar CBX has also been found in the tetragonal model with the different ordering wave vectors; in contrast to the present CBX state, another CBX state in the previous study was characterized by the quadruple-Q superposition at ${\mathit{Q}}_1^{\prime }=(\pi /4,3\pi /4)$, ${\mathit{Q}}_2^{\prime }=(3\pi /4,-\pi /4)$, ${\mathit{Q}}_3^{\prime }=(3\pi /4,\pi /4)$ and ${\mathit{Q}}_4^{\prime }=(-\pi /4,3\pi /4)$. In this way, various types of the single-Q and double-Q states are energetically degenerate within the bilinear exchange interaction for ${\mathit{Q}}_1=(\pi /2,\pi /2)$ and ${\mathit{Q}}_2=(-\pi /2,\pi /2)$.

3.2. With Biquadratic Interaction

In order to lift the degeneracy between the single-Q and double-Q states, we take into account the effect of the biquadratic interaction K. Since positive K tends to lead to the energy loss for the single-Q state compared to the multiple-Q state [71,72], the stabilization of the double-Q state is expected. Meanwhile, negative K tends to stabilize the single-Q state rather than the multiple-Q state. Accordingly, we focus on the effect of positive $K=0.05$ in the following analysis.

Figure 2 shows the magnetic phase diagram while varying $I^z$ and H at $K=0.05$. There are mainly two magnetic phases except for the fully polarized state appearing for $H\gtrsim 2$. For $I^z=0$, the double-Q spin configuration in Equation (7) is realized irrespective of H; we denote it as 2Q I in the phase diagram in Figure 2. The real-space spin configuration and spin structure factor of the 2Q I state are presented in the left and right panels of Figure 3a, respectively. In order to gain the energy by K, the intensities of the spin structure factor at ${\mathit{Q}}_1$ and ${\mathit{Q}}_2$ become almost equal; the isotropic double-Q state has a lower energy than the anisotropic double-Q state. The magnetization in this state exhibits the linear development against H, and continuously turns into the fully polarized state, as shown in Figure 4a.

For $I^z\ne 0$, the 2Q CBX is stabilized in the low-field region, which indicates that the energy of the 2Q CBX is lower than that of the single-Q collinear state in the presence of K. The real-space spin configuration and spin structure factor are shown in Figure 3b, which are the same as those in Figure 1c. When H increases, the magnetization remains zero in the 2Q CBX, as shown in the case of $I^z=0.1$ in Figure 4b. With a further increase in H, the 2Q CBX changes into the 2Q I state with a jump of M at $H\simeq 0.6$. This first-order phase transition is regarded as the spin–flop transition, where the spin moments are flopped from the parallel direction to the perpendicular direction in terms of the magnetic field. For larger $I^z$, the region where the 2Q CBX is stabilized becomes wider, and then, the magnetization jump becomes larger, as shown in the case of $I^z=0.3$ in Figure 4c and $I^z=0.48$ in Figure 4d.

Finally, let us comment on the intermediate magnetic phase (2Q II) appearing in the narrow region sandwiched by the 2Q CBX and 2Q I state for large $I^z$, as shown in Figure 2. This state is characterized by the isotropic double-Q superposition in both $xy$- and z-spin components like the magnetic skyrmion; the isotropic double-Q feature is found in the spin structure factor in the right panel of Figure 3c. Meanwhile, the spin configuration does not exhibit the scalar spin chirality corresponding to the topological charge, which means that the 2Q II state is topologically trivial in contrast to the magnetic skyrmion. We show the real-space spin configuration of the 2Q II state in the left panel of Figure 3c; there is no vortex with the nonzero winding number. The magnetization in this state takes almost the median value between the 2Q CBX and 2Q I states, as shown in Figure 4d, which indicates that two first-order phase transitions occur when the 2Q II state appears.

4. Conclusions

We have investigated the instability toward a new type of multiple-Q states in centrosymmetric tetragonal magnets by focusing on the position of the ordering wave vectors. We have shown that the single-Q and double-Q states are energetically degenerate within the easy-axis anisotropic bilinear exchange interaction and the magnetic field when the ordering wave vectors lie on $(\pi /2,\pi /2)$ and $(-\pi /2,\pi /2)$. By performing the simulated annealing and constructing the magnetic phase diagram, we have found that two types of the double-Q states, i.e., 2Q CBX and 2Q I state, are mainly realized in the presence of the biquadratic interaction. In particular, the low-field phase corresponds to the 2Q CBX state characterized by the double-Q collinear spin texture without a net magnetization. We have also found that another double-Q state denoted as $2Q$ II appears in the narrow region sandwiched by the low-field 2Q CBX and high-field 2Q I states. The present results indicate that the position of the ordering wave vectors might lead to further intriguing multiple-Q states distinct from the magnetic skyrmion and other well-known multiple-Q states.