Magnetic excitations of a trilayer antiferromagnetic Heisenberg model

Abstract

We investigate the squared sublattice magnetizations and magnetic excitations of a S=1/2 trilayer antiferromagnetic Heisenberg model with interlayer interaction J_{\mathrm{\perp }} and intralayer interaction J_{//} by employing stochastic series expansion quantum Monte Carlo (SSE-QMC) and stochastic analytic continuation (SAC) methods. Compared with the bilayer model, the trilayer model has one inner layer and two outer layers. The change in its symmetry can lead to special magnetic excitations. Our study reveals that the maximum of the magnetization of the outer sublattice corresponds to smaller ratio parameter g=J_{//}/J_{\mathrm{\perp }}, a finding that is verified using the finite-size extrapolation. As g decreases, the excitation spectra gradually evolve from a degenerate magnon mode with continua to low-energy and high-energy branches. Particularly when g is small enough, like 0.02, the high-energy spectrum further splits into characteristic doublon (\approx J_{\mathrm{\perp }}) and quarton (\approx 1.5J_{\mathrm{\perp }}) spectral bands. Moreover, the accuracy of the magnetic excitations is confirmed through the SpinW software package and the dispersion relations derived through the linear spin wave theory. Our results provide an important reference for experiments, which can be directly compared with experimental data from inelastic neutron scattering results to verify and guide the accuracy of experimental detection.

Graphical abstract

Keywords

quantum Monte Carlo / trilayer / magnetic excitation / magnetization / Heisenberg model

1 Introduction

Low-dimensional antiferromagnets with spin S=1/2 exhibit a variety of complex phenomena in the field of condensed matter physics as a result of their intense quantum fluctuations [1–3]. The complicated interplay between quantum mechanics and magnetic interactions is illustrated by the various kinds of ordered and disordered ground states that these systems can generate [4, 5]. The dimensionality crossover effect [6, 7] is one of the most intriguing features of these systems. The transition from one-dimensional (1D) chains [8, 9] to ladder structures [10] and two-dimensional (2D) systems [11] is characterized by variations in dimensionality that result in changes in magnetism. The impact of interlayer coupling must be taken into account as we transition from 2D systems [12, 13] to three-dimensional (3D) systems, whether it is a finite layering of 2D or a genuine three-dimensional structure. This complicates the issue and has, as a result, been less extensively studied.

As observed in bilayer systems, the stacking of two-dimensional layers offers additional degrees of freedom and may result in distinctive magnetic characteristics. At absolute zero, Néel ordered phases and magnetically disordered phases are identified, with the antiferromagnetic order−disorder transition point established at ( J_{\mathrm{\perp }}/J_{//}{)}_c= 2.5220(1) [14, 15]. Even for the more intricate dimer-diluted bilayer square lattice, by changing the interaction ratio J_{\mathrm{\perp }}/J_{//}, one can detect a Néel ordered phase, a gapless quantum glass phase, and a gapped quantum paramagnetic phase [16]. In addition, it has been seen that the bilayer materials {\mathrm{B}\mathrm{a}\mathrm{C}\mathrm{u}\mathrm{S}\mathrm{i}}_2{\mathrm{O}}_6, {\mathrm{C}\mathrm{s}}_3{\mathrm{C}\mathrm{r}}_2{\mathrm{B}\mathrm{r}}_9 and {\mathrm{M}\mathrm{n}\mathrm{B}\mathrm{i}}_2{\mathrm{T}\mathrm{e}}_4 have a spin-singlet dimerized ground state and a singlet-to-triplet dimer excitation [17, 18] that are different from those in monolayer materials. Under high pressures, ranging from 14.0 GPa to 43.5 GPa, the bilayer structure of {\mathrm{L}\mathrm{a}}_3{\mathrm{N}\mathrm{i}}_2{\mathrm{O}}_7 single crystals has been found to exhibit superconducting characteristics [19–21].

When there are three layers, the z-axis lacks symmetry, making calculations difficult. The trimer model, a simplified trilayer structure, is a good study starting point. The energy spectra of the 1D antiferromagnetic trimer chain highlights two key excitation modes: the doublon mode in the intermediate-energy regime and the quarton mode in the high-energy regime [22]. The same was also found in the 2D trimeric system [23]. The ground state of a trilayer is a trivial ordered state [24–26]; however, its excited states demonstrate a variety of physical properties by adjusting the interaction strength between layers, including compounds {\mathrm{B}\mathrm{i}}_2{\mathrm{S}\mathrm{r}}_2{\mathrm{C}\mathrm{a}}_2{\mathrm{C}\mathrm{u}}_3{\mathrm{O}}_{10+x}, {\mathrm{Y}\mathrm{B}\mathrm{a}}_2{\mathrm{C}\mathrm{u}}_3{\mathrm{O}}_{6+ϵ}, {\mathrm{H}\mathrm{g}\mathrm{B}\mathrm{a}}_2{\mathrm{C}\mathrm{a}}_2{\mathrm{C}\mathrm{u}}_3{\mathrm{O}}_8 and {\mathrm{C}\mathrm{r}\mathrm{I}}_3 [27–29]. The discovery of high-temperature superconductivity in the bilayer nickelate {\mathrm{L}\mathrm{a}}_3{\mathrm{N}\mathrm{i}}_2{\mathrm{O}}_7 [19] has sparked interest in {\mathrm{L}\mathrm{a}}_4{\mathrm{N}\mathrm{i}}_3{\mathrm{O}}_{10}, a compound with magnetic particles and a trilayer structure. {\mathrm{L}\mathrm{a}}_4{\mathrm{N}\mathrm{i}}_3{\mathrm{O}}_{10} can induce a trimeric lattice, with metal density waves leading to an unusual metal-to-metal transition [30, 31]; however, the results of its inelastic neutron scattering measurements have not yet been published.

Our research investigated the trilayer antiferromagnetic Heisenberg model shown in Fig.1. Both numerical simulation and theoretical analysis produced consistent findings. To account for intralayer and interlayer interactions, we set J_{\mathrm{\perp }}=1 and establish a tuning parameter g=J_{//}/J_{\mathrm{\perp }}. The squared sublattice magnetizations indicate that the outer layer has the biggest magnetization at g=1, whereas the inner layer has it at g=\frac{5}{3} at the thermodynamic limit. Moreover, we plotted the magnetic excitation spectra. As g decreases, the excitation spectrum evolves from a single magnon spin wave to a low-energy spin wave, intermediate-energy doublon flat band (\omega =1), and high-energy quarton flat band (\omega =1.5), which become distinguishable when g is reduced to 0.02. Studying trilayer magnetic excitations provides a theoretical foundation for trilayer magnetic materials and sheds light on high-temperature superconductivity.

The remaining parts follow this pattern. Section 2 outlines the trilayer Heisenberg model composition and the numerical and theoretical methodologies used for measurement. In Section 3, we present the numerical results of the model, including the squared sublattice magnetizations and magnetic excitation spectra. In Section 4, we analyze theoretically calculated dispersion relations. In the final Section 5, we provide our study findings and tentative plans for future research.

2 Model and method

2.1 Model Hamiltonian

The Hamiltonian of this S=1/2 trilayer antiferromagnetic Heisenberg model is given by

(1)H=J_{//}{\displaystyle \sum _{⟨i,j⟩}{\mathit{S}}_i\cdot {\mathit{S}}_j+J_{\mathrm{\perp }}{\displaystyle \sum _{⟨i,j{⟩}^{\mathrm{\prime }}}{\mathit{S}}_i\cdot {\mathit{S}}_j,}}

where S_i denotes the spin operator on each site i, ⟨i,j⟩ denotes nearest-neighbor sites on the intralayer bonds, and ⟨i,j{⟩}^{\mathrm{\prime }} denotes the interlayer bonds. We use J_{//} and J_{\mathrm{\perp }} to indicate intralayer and interlayer coupling strengths. This is seen when J_{\mathrm{\perp }}=0 (g=\mathrm{\infty }), the trilayer model degenerates into three 2D square lattices that have been studied in great detail, as shown in Fig.2(a). In contrast, when J_{//}=0 (g=0), only J_{\mathrm{\perp }} in the vertical direction remains, resulting in non-contact trimers in the horizontal direction, as shown in Fig.2(b). Trimers may now be considered a single, bigger spin, changing the model into a renormalized 2D plane with S=1/2.

2.2 Quantum Monte Carlo

The stochastic series expansion quantum Monte Carlo (SSE-QMC) approach makes physical quantity measurements easier and avoids systematic mistakes that might lower accuracy. The sign problem does not apply to our model since it only considers the closest neighbor interactions.

SSE-QMC works by Taylor expanding the partition function, as shown in the formula:

(2)Z=Tr\left\{{\mathrm{e}}^{-\beta H}\right\}= {\displaystyle \sum _{\alpha }{\displaystyle \sum _{n=0}^{\mathrm{\infty }}{\displaystyle \frac{{\left(\beta \right)}^n}{n!}}⟨\alpha \left|{\left(-H\right)}^n\right|\alpha ⟩.}}

After several updates, the system reaches equilibrium and is then amenable to representation by operators for physical observations like ground-state energy and spin correlation.

In addition to thermodynamic quantities, the SSE-QMC approach can calculate the dynamic spin structure factor. SSE-QMC sampling is used to obtain imaginary-time correlation functions. These functions are defined by

(3)G_q(\tau )=3⟨S_{-q}^z(\tau )S_q^z(0)⟩.

The spin operators S^x, S^y, and S^z are equivalent in the Hamiltonian, resulting in a factor of 3 in it. Since the trilayer model lacks periodic boundary conditions in the z-direction, Fourier expansion cannot be applied in this direction. We must compute G_q(\tau ) for each layer separately, then sum the layer values to get the G_q(\tau ) of the overall model.

The term S_q^z represents the Fourier transform of the z-component spin operator for each layer, which is expressed as

(4)S_q^z={\displaystyle \frac{1}{\sqrt{N}}}{\displaystyle \sum _{i=1}^N{\mathrm{e}}^{-\mathrm{i}{\mathit{r}}_i\cdot \mathit{q}}{S}_i^z.}

We use the method of stochastic analytic continuation (SAC) to reconstruct the spectral function S\left(q,\omega \right) from a sequence of imaginary-time points for G_q(\tau ):

(5)G_q(\tau )={\displaystyle \frac{1}{\pi }}{\displaystyle {\int }_{-\mathrm{\infty }}^{\mathrm{\infty }}\mathrm{d}\omega S(q,\omega ){\mathrm{e}}^{-\tau \omega }.}

We parametrize the spectrum using N_{\omega } delta functions in the continuum:

(6)S(q, \omega )= {\displaystyle \sum _{N_{\omega }}^{i=0}{a}_i\delta (\omega -{\omega }_i).}

2.3 Linear spin wave theory

The dispersion relations of the trilayer antiferromagnetic Heisenberg model are calculated using the linear spin-wave theory (LSWT).

To ensure reproducibility, we chose six lattices to construct a macromolecular cell (red circle in Fig.1).

For this model, the Holstein−Primakoff boson method can be used to calculate the dispersion relation [32]. We quantize the above Hamiltonian using Holstein−Primakoff bosons in the momentum space and disregard the higher-order components to achieve the low-energy approximations:

(7)H=E_{Cl} +{\displaystyle \sum _{k,⟨i,j⟩}{C}_{ii}{a}_{k,i}^{†}{a}_{k,i}+C_{ij}(a_{k,i}{a}_{k,j}+a_{k,i}^{†}{a}_{k,j}^{†}),}

where E_{Cl}=(- 12J_{//}-4 J_{\mathrm{\perp }})N{S}^2 is the ground state energy and i,j are spin indexes in the unit cell. Since the Hamiltonian contains off-diagonal terms, we diagonalize it using the extended Bogoliubov transformation:

(8)b_{k,i}={\displaystyle \sum _j{m}_{ij}{a}_{k,j}+m_{ij}^{\prime }{a}_{k,j}^{†}.}

Only one quasiparticle may result from the transformation due to the symmetry of the Hamiltonian:

(9)\begin{array}{rl}{b}_{k,i}=& m_{11}{a}_{k1}+m_{12}{a}_{k2}^{†}+m_{13}{a}_{k3}+m_{14}{a}_{k4}^{†}\\ & +m_{15}{a}_{k5}+m_{16}{a}_{k6}^{†}.\end{array}

The transformation coefficients m_{ki} can be obtained by setting the determinant zero from the equation of motion [32–34]

(10)\mathrm{i}\hslash {\dot{b}}_{k,i}=[b_{k,i},H ].

The diagonalized Hamiltonian is

(11)H=E_{Cl} +E_0+{\displaystyle \sum _{k,i}{\omega }_i(k)b_{k,i}^{†}{b}_{k,i},}

where {\omega }_i(k ) is the spin wave dispersion, and E_0 is the quantum zero-point energy correction.

3 Numerical results

In our QMC-SAC calculations, each layer of the trilayer Heisenberg model has dimensions L=L_x=L_y, with both L_x and L_y being even numbers, resulting in a total size of N=3L^2. After rigorous convergence testing, we set the inverse temperature to \beta = 6L, and each data set is equipped with 500 bins to reduce systematic error. These measures ensure the precision and reliability of the numerical calculations.

3.1 Sublattice magnetizations

In this section, we present the squared sublattice magnetizations of J_{\mathrm{\perp }} in different cases of the trilayer model, along with their extrapolated results.

The squared sublattice magnetizations m_s^2 are defined as

(12)m_s^2={\displaystyle \frac{3}{N^{{\mathrm{\prime }}^2}}} ⟨{\left({\displaystyle \sum _{i=1}^{N^{\mathrm{\prime }}}{\varphi }_i{S}_i^z}\right)}^2⟩,

where {\varphi }_i=\pm 1 signifies staggered phase factors, while the factor of 3 is attributed to the isotropic strength encompassing all three components of the spin. For the complete model, N^{\mathrm{\prime }}=3L^2, while for a single layer, N^{\mathrm{\prime }}=L^2.

Fig.3 shows the squared sublattice magnetizations versus g (with J_{\mathrm{\perp }}=1) for a system size of L=24. For any layer, m_s^2 first rises and then falls as g rises. After fitting three sets of data using a second-order polynomial function, we discovered that the inner maximum is around g=\frac{5}{3} and the outer maximum is about g=1. Additionally, the maximum of the trilayer is centered at g=\frac{10}{9}. Improved intralayer antiferromagnetic interactions enhance long-range antiferromagnetic order, with g increasing from 0 to roughly 1. The model shows better alignment in the three directions when g approaches 1, signifying the highest magnetic order. As the parameter g increases sufficiently, the model evolves into three almost isolated 2D square lattices. The reduction of dimensions leads to a decrease in long-range antiferromagnetism. The magnetizations of the outer layers are higher than those of the inner layer when g is between 0 and 1 (the black line is higher than the red line). However, the conclusion is reversed from 1 to 5. The inner layer spins have a coordination number of 6, with 4 horizontal and 2 vertical neighbors. In contrast, the outer layer spins have a coordination number of 5, with 4 horizontal and 1 vertical neighbors. As g increases, the inner layer spins are more significantly impacted.

To precisely determine the g corresponding to the extremum points for the three situations, we use finite-size extrapolation for sizes L = 16, 24, 30, 40, 50, and 60. The y-axis intercepts on Fig.3(b)−(d) show the thermodynamic limit values derived from second-order polynomial fittings. These match the pattern shown in Fig.3(a): the maximum magnetization in the outer layer [m_s^2= 0.1494(4)] is observed at g=1, while the maximum magnetization in the inner layer [m_s^2= 0.1567(5)] occurs at g=\frac{5}{3}, and the entire trilayer model reaches its maximum magnetization [m_s^2= 0.1510(2)] when g=\frac{10}{9}. As the model approaches the trimerization structure (g=0.2, light blue line), the squared sublattice magnetization of the outer layers [m_s^2=0.0809(3)] in the thermodynamic limit exceeds that of the inner layer [m_s^2= 0.0457(2)]. As g decreases to 0, the model gradually evolves into a vertical trimerization structure. Under these conditions, the inner layer spins pair with the coordinated spins of either the upper or lower layer to form a singlet, resulting in an unpaired spin in the outer layer. The unpaired spins yield greater magnetization in the outer sublattice.

The squared sublattice magnetization primarily reflects the ground-state properties of the model. To conduct a more comprehensive study of the model, we further analyze its excited states by calculating the dynamic structure factor. The results indicate that the maximum of the squared sublattice magnetization occurs at g\approx 1 [see Fig.3(a)], which also gives more stable spin wave excitations, as shown in Fig.4(b).

3.2 Magnetic excitations

In this subsection, we present the magnetic excitation spectra of the trilayer antiferromagnetic Heisenberg model via QMC-SAC calculations. We set system size L=24, inverse temperature \beta =6L=144. To visually represent the S\left(q,\omega \right), we follow a high-symmetry trajectory in the Brillouin zone: \mathrm{\Gamma }( 0,0) \to \mathrm{X}(\pi ,0)\to \mathrm{M}(\pi , \pi )\to\mathrm{\Gamma }(0,0 )\to T(0, \pi )\to X(\pi ,0). The above path is depicted in Fig.5. Fig.4 demonstrates how the dynamic spin structure factor S\left(q,\omega \right) varies with changes in g. We choose g at \mathrm{\infty } (J_{\mathrm{\perp }}=0), 1, 0.6, 0.4, 0.3, 0.1, 0.02, and 0 to illustrate the advancement. Due to the difference of three orders of magnitude between the minimum and maximum of S(q, \omega ), a piecewise function was applied to clearly see both areas of high and low intensity simultaneously in Fig.4. The low-intensity region has a linear distribution for values below U_0. A logarithmic scale is applied when the value surpasses this threshold, as denoted by U=U_0+{\log }_{10}(S(q,\omega ))-{\log }_{10}(U_0). After multiple attempts, we choose U_0=6. This outcome can be compared to the neutron scattering excitation spectra of material with a trilayer structure to confirm the accuracy of the theoretical results.

First, we must verify the findings of g=J_{//}/J_{\mathrm{\perp }}=\mathrm{\infty } as presented in Fig.4(a). The trilayer Heisenberg model separates into three unconnected single-layer 2D square lattices, as illustrated in Fig.2(a). The dispersion relation, indicative of the expansion of linear spin waves [32], has been determined by previous researchers and is expressed as follows:

(13){\omega }_{AF}=2 J_{//}S\sqrt{4-(\cos k_x+\cos k_y{)}^2.}

This equation aligns with the acoustic magnon band in our result for LSWT, as shown by the blue line in Fig.4(a). The spectra meet with the established pattern of linear spin waves, where the points \mathrm{\Gamma } and \mathrm{M} show gapless modes and the band top approaches \omega =2. The LSWT approximation causes band top matching differences. The LSWT cannot describe the high-energy spectra obtained by QMC-SAC. Since LSWT ignores several high-order factors throughout the calculation process, the theoretical results are slightly lower than actual results in the high-energy excitation region [23, 35]. The spectra display a gapless Goldstone mode at the point \mathrm{M}(\pi ,\pi ), distinguished by maximized spectral weight shown in a bright yellow color and indicating divergence. The uneven and arched excitation from \mathrm{T}(0,\pi ) to \mathrm{X}(\pi ,0) is due to the decay of spin waves into other excitations at T and X, which suppresses the high-energy excitations at these two points. This has been analyzed in many previous results [36, 37]. All findings fit with the SAC spectra observed in previous research [35]. It is important to acknowledge that the gapless point at \mathrm{M}(\pi ,\pi ) requires a very large \beta to achieve convergence results. Consequently, we did not include them in Fig.4. Nevertheless, we can infer their behavior from the points surrounding it.

Next, from Fig.4(a) and (b), J_{\mathrm{\perp }} decreases from \mathrm{\infty } to 1. The number of internal interactions in the model rises dramatically throughout this process, and its structure moves from a 2D plane to a quasi-two-dimensional one. This change results in a higher excitation energy and an expansion of the bandwidth. As g continues to decrease, the band top systematically declines, but the spectral weight increases. This singular continuous spectrum begins to split into a high-energy part and a low-energy branch as g\approx 0.3. The dispersion relation derived from the linear spin wave method (shown by the blue line in Fig.4) accurately corresponds to the low-energy spin wave spectrum in the bottom part. The high-energy continuous spectrum reflects the characteristics of the internal coupling within the trimer. The spectral weight is directly proportional to J_{\mathrm{\perp }}, resulting in a gradual increase, while the broadening concurrently reduces. This evidence indicates that high-energy quantum excitations can only be significantly observed when the value of J_{\mathrm{\perp }} is comparatively large. When g\approx 0.02, the high-energy part clearly splits into two flat bands, and three distinctly separated spectral lines are seen in Fig.4(g). The characteristic spectrum around \omega \approx 1, associated with the quasiparticle, is termed the doublon, while the spectrum around \omega \approx 1.5, corresponding to another type of quasiparticle, is called the quarton [22].

Finally, let us discuss another extreme case where g=J_{//}/J_{\mathrm{\perp }}=0, and only the vertical coupling strength within the unit cell J_{\mathrm{\perp }} remains. Currently, the model consists of isolated vertical trimers, each can be considered as an effective spin of S=1/2, as illustrated in Fig.2(b). The excitations within the trimers are localized and without any magnons. In Fig.4(h), the classical spin wave excitations disappear. Near \omega \approx 1 and \omega \approx 1.5, two continuous spectral bands exhibit minimal energy broadening and significant power weight. These bands represent the internal excitation states of the trimer. These two characteristic spectral lines perfectly match previous studies [22], not only demonstrating the accuracy of our SAC calculations but also elucidating the similarity between trimer chains and isolated vertical trimers. The ground state of a trimer is a doublet with energy E_0= -J_{\mathrm{\perp }}, whereas the first and second excited states are a doublet with energy E_1=0 and a quartet with energy E_2=J_{\mathrm{\perp }}/2, respectively. The two low-energy doublets include only singlets and unpaired spins, while the higher-energy quadruplet excitations include either a triplet pair with an unpaired spin or three unpaired spins collectively.

4 Analysis

As mentioned earlier, we have selected six lattice points as a macromolecular cell, so we can obtain three dispersion relationships: the acoustic mode ({\omega }_1), the lower energy optical mode ({\omega }_2), and the higher energy optical mode ({\omega }_3). The spin wave dispersion \omega is given by

(14)\{\begin{array}{l}{\omega }_1=S\sqrt{A_k-B_k},\\ {\omega }_2=S\sqrt{(4J_{//}+J_{\mathrm{\perp }}{)}^2-(2J_{//}\cos k_x+2J_{//}\cos k_y{)}^2},\\ {\omega }_3=S\sqrt{A_k+B_k},\end{array}

where

(15)\begin{array}{rl}& A_k=(4J_{//}{)}^2+\frac{J_{\mathrm{\perp }}^2}{2}+12J_{//}{J}_{\mathrm{\perp }}-4J_{//}^2(\cos k_x+\cos k_y{)}^2,\\ & B_k=4J_{\mathrm{\perp }}\sqrt{2J_{//}^2(\cos k_x+\cos k_y{)}^2+J_{//}^2+\frac{J_{\mathrm{\perp }}^2}{64}+\frac{3}{4}{J}_{//}{J}_{\mathrm{\perp }}}.\end{array}

The final results indicate that the model shows a certain degree of decoupling, with the inner and outer layers exhibiting identical dispersion relations. The three dispersion relationships represent acoustic and optical modes. The acoustic mode depicts the movement of cell centroids, whereas the optical mode describes the relative motion between spins within the cell.

In Fig.6, the theoretically derived dispersion relations show excellent agreement with the SpinW simulation results. In the limit of g\to \mathrm{\infty }(g=J_{//}/J_{\mathrm{\perp }}), the three dispersion relations converge to those described by Eq. (13) of the two-dimensional square lattice. As g decreases towards 1, the model becomes more three-dimensional. This change is accompanied by an increase in the influence of interlayer coupling, leading to an overall increase in excitation energies. At g=1, the highest energy mode ({\omega }_3) separates significantly from the other two lower-energy modes. As shown in Fig.6(a), a clear separation between the high-energy and low-energy regions is observed. Fig.6(b) illustrates that at g=0.3, a distinction emerges between the optical and acoustic modes. Notably, the spin wave dispersion becomes a flat band when g=0.

Theoretical acoustic modes ({\omega }_1) exhibit consistency with the numerical results, but the optical modes do not match the high-energy portion of the numerical calculations. Since the calculation only considers the lowest-order terms in the Hamiltonian, the resulting dispersion is a low-energy approximation, accurately describing the spin excitation spectra only in the low-energy regime.

The associated spin-wave velocities are

(16)\begin{array}{rl}& v_x=v_y=\\ & \sqrt{48J_{//}^2+2J_{\mathrm{\perp }}^2+48J_{//}{J}_{\mathrm{\perp }}-2J_{\mathrm{\perp }}\sqrt{192J_{//}^2+48J_{//}{J}_{\mathrm{\perp }}+J_{\mathrm{\perp }}^2}}.\end{array}

The spin-wave velocity exhibits a linear dependence on g. Consequently, as J_{//} increases, the wave velocity increases, and as J_{\mathrm{\perp }} increases, it decreases. These results are consistent with the theoretical prediction that strong J_{\mathrm{\perp }} favors the formation of trimer bound states, whereas strong J_{//} facilitates spin-wave propagation within the horizontal plane.

A key issue of the linear spin wave technique is to use operator characteristics to turn the Schrödinger equation into a linear system. Compared to other approaches, the linear spin wave method is easier to compute, simpler to program, and yields an exact analytical answer. However, it can only calculate a few physical quantities and is suitable only for low-energy approximations.

5 Summary and discussion

In this study, we investigated magnetic excitations of the trilayer antiferromagnetic Heisenberg model (Fig.1). We computed the squared sublattice magnetizations of each layer and found that the value of g corresponding to the maximum magnetization of the outer layer is smaller than that of the inner layer. We employed finite-size extrapolation to verify this discovery. Subsequently, we calculated the dynamic spin structure factor of the whole model using SSE-QMC and SAC methods to investigate the magnetic excitation spectra. We observed that the spin wave excitation spectrum separates into low-energy and high-energy sections at g\approx 0.3 as the parameter g decreases. The high-energy component further splits into doublon and quarton branches when g is small. To validate our numerical results, we derived the spin wave dispersions of the trilayer Heisenberg model from the LSWT. The resulting dispersion relations agree well with the SpinW simulations, exhibiting one low-energy acoustic branch and two high-energy optical branches. In the model, when g is sufficiently large, the LSWT results match well with the SSE-QMC and SAC calculations. However, as g approaches 0, doublon and quarton excitations become increasingly dominant, resulting in a worse alignment between them. Notably, the low-energy acoustic branch mainly coincides with the low-energy spin wave component of the excitation spectra, as predicted by the semiclassical spin-wave theory.

Going forward, our theoretical conclusions may be confirmed by inelastic neutron scattering measurements on {\mathrm{B}\mathrm{i}}_2{\mathrm{S}\mathrm{r}}_2{\mathrm{C}\mathrm{a}}_2{\mathrm{C}\mathrm{u}}_3{\mathrm{O}}_{10+x}, {\mathrm{H}\mathrm{g}\mathrm{B}\mathrm{a}}_2{\mathrm{C}\mathrm{a}}_2{\mathrm{C}\mathrm{u}}_3{\mathrm{O}}_8, and {\mathrm{L}\mathrm{a}}_4{\mathrm{N}\mathrm{i}}_3{\mathrm{O}}_{10} [38, 39]. In other trilayer materials, quantum excitations linked with spin waves should be observable when vertical interactions are substantial. Additionally, ultracold atoms in optical lattices, Rydberg atom systems, and quantum networks offer exciting options for researching the trilayer model beyond the current limits [40–42].

References

[1] N. D. Mermin and H. Wagner, Absence of ferromagnetism or antiferromagnetism in one- or two-dimensional isotropic Heisenberg models, Phys. Rev. Lett. 17(22), 1133 (1966)

[2] H. J. Deiseroth, K. Aleksandrov, C. Reiner, L. Kienle, and R. K. Kremer, Fe₃GeTe₂ and Ni₃GeTe₂ – two new layered transition-metal compounds: Crystal structures, HRTEM investigations, and magnetic and electrical properties, Eur. J. Inorg. Chem. 2006(8), 1561 (2006)

[3] F. Lu, L. Yuan, J. Zhang, B. Li, Y. Luo, and Y. Li, The observation of quantum fluctuations in a kagome Heisenberg antiferromagnet, Commun. Phys. 5(1), 272 (2022)

[4] K. Kim, S. Y. Lim, J. U. Lee, S. Lee, T. Y. Kim, K. Park, G. S. Jeon, C. H. Park, J. G. Park, and H. Cheong, Suppression of magnetic ordering in XXZ-type antiferromagnetic monolayer NiPS₃, Nat. Commun. 10(1), 345 (2019)

[5] K. Lee, A. H. Dismukes, E. J. Telford, R. A. Wiscons, J. Wang, X. Xu, C. Nuckolls, C. R. Dean, X. Roy, and X. Zhu, Magnetic order and symmetry in the 2D semiconductor CrSBr, Nano Lett. 21(8), 3511 (2021)

[6] Z. Wang, M. Gibertini, D. Dumcenco, T. Taniguchi, K. Watanabe, E. Giannini, and A. F. Morpurgo, Determining the phase diagram of atomically thin layered antiferromagnet CrCl₃, Nat. Nanotechnol. 14(12), 1116 (2019)

[7] M. F. DiScala, D. Staros, A. de la Torre, A. Lopez, D. Wong, C. Schulz, M. Barkowiak, V. Bisogni, J. Pelliciari, B. Rubenstein, and K. W. Plumb, Elucidating the role of dimensionality on the electronic structure of the van der Waals antiferromagnet NiPS₃, Adv. Phys. Res. 3(4), 2300096 (2024)

[8] J. des Cloizeaux and J. J. Pearson, Spin-wave spectrum of the antiferromagnetic linear chain, Phys. Rev. 128(5), 2131 (1962)

[9] S. Nagler, D. Tennant, R. Cowley, T. Perring, and S. Satija, Spin dynamics in the quantum antiferromagnetic chain compound KCuF₃, Phys. Rev. B 44(22), 12361 (1991)

[10] M. Reigrotzki, H. Tsunetsugu, and T. Rice, Strongcoupling expansions for antiferromagnetic Heisenberg spin-one-half ladders, J. Phys.: Condens. Matter 6(43), 9235 (1994)

[11] G. Chen, H. Q. Ding, and Goddard , Elementary excitations for the two-dimensional quantum Heisenberg antiferromagnet, Phys. Rev. B 46(5), 2933 (1992)

[12] J. Richter and J. Schulenburg, The spin-1/2 J₁–J₂ Heisenberg antiferromagnet on the square lattice: Exact diagonalization for N = 40 spins, Eur. Phys. J. B 73(1), 117 (2010)

[13] L. Wang and A. W. Sandvik, Low-energy excitations of two-dimensional diluted Heisenberg quantum antiferromagnets, Phys. Rev. B 81(5), 054417 (2010)

[14] W. H. Zheng, Various series expansions for the bilayer S= 1/2 Heisenberg antiferromagnet, Phys. Rev. B 55, 12267 (1997)

[15] L. Wang, K. Beach, and A. W. Sandvik, High-precision finite-size scaling analysis of the quantum-critical point of S= 1/ 2 Heisenberg antiferromagnetic bilayers, Phys. Rev. B 73, 014431 (2006)

[16] N. S. Ma, A. W. Sandvik, and D. X. Yao, Mott glass phase in a diluted bilayer Heisenberg quantum antiferromagnet, J. Phys.: Conf. Ser. 640, 012045 (2015)

[17] Y. Sasago, K. Uchinokura, A. Zheludev, and G. Shirane, Temperature-dependent spin gap and singlet ground state in BaCuSi₂O₆, Phys. Rev. B 55(13), 8357 (1997)

[18] B. Leuenberger, A. Stebler, H. Güdel, A. Furrer, R. Feile, and J. Kjems, Spin dynamics of an isotropic singlet-ground-state antiferromagnet with alternating strong and weak interactions: An inelastic-neutron-scattering study of the dimer compound Cs₃Cr₂Br₉, Phys. Rev. B 30(11), 6300 (1984)

[19] H. Sun, M. Huo, X. Hu, J. Li, Z. Liu, Y. Han, L. Tang, Z. Mao, P. Yang, B. Wang, J. Cheng, D. X. Yao, G. M. Zhang, and M. Wang, Signatures of superconductivity near 80 k in a nickelate under high pressure, Nature 621(7979), 493 (2023)

[20] Z. Luo, X. Hu, M. Wang, W. Wu, and D. X. Yao, Bilayer two-orbital model of La₃Ni₂O₇ under pressure, Phys. Rev. Lett. 131(12), 126001 (2023)

[21] J. Li, Y. Li, S. Du, Z. Wang, B. L. Gu, S. C. Zhang, K. He, W. Duan, and Y. Xu, Intrinsic magnetic topological insulators in van der Waals layered MnBi₂Te₄-family materials, Sci. Adv. 5(6), eaaw5685 (2019)

[22] J. Q. Cheng, J. Li, Z. Xiong, H. Q. Wu, A. W. Sandvik, and D. X. Yao, Fractional and composite excitations of antiferromagnetic quantum spin trimer chains, npj Quantum Mater. 7, 3 (2022)

[23] Y. Y. Chang, J. Q. Cheng, H. Shao, D. X. Yao, and H. Q. Wu, Magnon, doublon and quarton excitations in 2D S= 1/2 trimerized Heisenberg models, Front. Phys. (Beijing) 19(6), 63202 (2024)

[24] Y. Hu, F. Shi, N. Jia, Y. Liu, H. Wu, and A. Du, Exchange bias and its propagation in ferromagnetic/antiferromagnetic/ferromagnetic trilayers, J. Appl. Phys. 114(15), 153901 (2013)

[25] S. Chatterjee, T. Wang, E. Berg, and M. P. Zaletel, Intervalley coherent order and isospin fluctuation mediated superconductivity in rhombohedral trilayer graphene, Nat. Commun. 13(1), 6013 (2022)

[26] L. Weber, A. Honecker, B. Normand, P. Corboz, F. Mila, and S. Wessel, Quantum Monte Carlo simulations in the trimer basis: First-order transitions and thermal critical points in frustrated trilayer magnets, SciPost Phys. 12, 054 (2022)

[27] G. Roth, B. Renker, G. Heger, M. Hervieu, B. Domenges, and B. Raveau, On the structure of non-superconducting YBa₂Cu₃O_6+ϵ, Z. Phys. B 69(1), 53 (1987)

[28] V. Shamray, A. Mikhailova, and A. Mitin, Crystal structure and superconductivity of Bi-2223, Crystallogr. Rep. 54(4), 584 (2009)

[29] D. Wang and B. Sanyal, Systematic study of monolayer to trilayer CrI₃: Stacking sequence dependence of electronic structure and magnetism, J. Phys. Chem. C 125(33), 18467 (2021)

[30] J. Zhang, D. Phelan, A. Botana, Y. S. Chen, H. Zheng, M. Krogstad, S. G. Wang, Y. Qiu, J. Rodriguez-Rivera, R. Osborn, S. Rosenkranz, M. R. Norman, and J. F. Mitchell, Intertwined density waves in a metallic nickelate, Nat. Commun. 11(1), 6003 (2020)

[31] Q. Qin, J. Wang, and Y. Yang, Frustrated superconductivity and intrinsic reduction of T_c in trilayer nickelate, The Innovation Mater. 2(4), 100102 (2024)

[32] D. X. Yao and E. W. Carlson, Magnetic excitations in the high-T_c iron pnictides, Phys. Rev. B 78(5), 052507 (2008)

[33] Y. Tan and D. X. Yao, Spin waves and phase transition on a magnetically frustrated square lattice with long-range interactions, Front. Phys. (Beijing) 18(3), 33309 (2023)

[34] E. Carlson, D. Yao, and D. Campbell, Spin waves in striped phases, Phys. Rev. B 70(6), 064505 (2004)

[35] H. Shao, Y. Q. Qin, S. Capponi, S. Chesi, Z. Y. Meng, and A. W. Sandvik, Nearly deconfined spinon excitations in the square-lattice spin-1/2 Heisenberg antiferromagnet, Phys. Rev. X 7(4), 041072 (2017)

[36] B. Dalla Piazza, M. Mourigal, N. B. Christensen, G. Nilsen, P. Tregenna-Piggott, T. Perring, M. Enderle, D. F. McMorrow, D. Ivanov, and H. M. Rønnow, Fractional excitations in the square-lattice quantum antiferromagnet, Nat. Phys. 11(1), 62 (2015)

[37] N. Headings, S. Hayden, R. Coldea, and T. Perring, Anomalous high-energy spin excitations in the high-T_c superconductor-parent antiferromagnet La₂CuO₄, Phys. Rev. Lett. 105(24), 247001 (2010)

[38] L. Abbas, M. Abbas, and O. Mahdi, Effect of sintering temperature on Bi_2−xCu_xPb_0.3Sr₂Ca₂Cu₃O_10+δ superconductors, J. Al-Nahrain Univ. Sci. 15(2), 113 (2012)

[39] F. W. Ali and K. A. Jasim, Study of some physical properties of the superconducting compound PbBa₂Ca₂Cu₃O_8+δ, Ibn AL-Haitham J. Pure Appl. Sci. 37(2), 168 (2024)

[40] C. S. Adams, J. D. Pritchard, and J. P. Shaffer, Rydberg atom quantum technologies, J. Phys. At. Mol. Opt. Phys. 53(1), 012002 (2020)

[41] G. Wirth, M. Olschlager, and A. Hemmerich, Evidence for orbital superfluidity in the p-band of a bipartite optical square lattice, Nat. Phys. 7(2), 147 (2011)

[42] C. Gross and I. Bloch, Quantum simulations with ultra-cold atoms in optical lattices, Science 357(6355), 995 (2017)

RIGHTS & PERMISSIONS

Higher Education Press