Insights on magnon topology and valley-polarization in 2D bilayer quantum magnets

The rich and unconventional physics in layered 2D magnets can open new avenues for topological magnonics and magnon valleytronics. In particular, two-dimensional (2D) bilayer quantum magnets are gaining increasing attention due to their intriguing stacking-dependent magnetism, controllable ground states, and topological excitations induced by magnetic spin–orbit couplings (SOCs). Despite the substantial research on these materials, their topological features remain widely unexplored to date. The present study comprehensively investigates the magnon topology and magnon valley-polarization in honeycomb bilayers with collinear magnetic order. We elucidate the separate and combined effects of the SOC, magnetic ground-states, stacking order, and inversion symmetry breaking on the topological phases, magnon valley transport, and the Hall and Nernst effects. The comprehensive analysis suggests clues to determine the SOC’s nature and predicts unconventional Hall and Nernst conductivities in topologically trivial phases. We further report on novel bandgap closures in layered antiferromagnets and detail their topological implications. We believe the present study provides important insights into the fundamental physics and technological potentials of topological 2D magnons.

The interest in 2D magnetic materials naturally extends to their magnetic excitations known as spin waves or magnons [10,25,26,34,[40][41][42]. SOCs in 2D magnets induce topological magnons [27, 30-32, 34, 35, 37-40, 43-55] with exotic features and great potentials for practical applications. Topological 2D magnets display a finite magnon Hall conductivity that can be harnessed in magnonic nanodevices. Equally important are their topologically protected edge and domain-wall modes, which are robust against structural or magnetic disorders and can form ideal waveguides (or quantum wires) for long-range and coherent magnon transport. Honeycomb CrI 3 ferromagnet, for example, represents a model platform to investigate 2D magnetic interactions and topological excitations. Recent studies on monolayer CrI 3 [25][26][27] reported large Kitaev couplings and topological magnonic bands with Chern numbers C = ±1 in the Heisenberg-Kitaev model [27]. The presence of DMI is also promoted in CrI 3 [24,26,34], leaving the Stacking layers of 2D ferromagnets is a systematic way to realize new topological quantum materials with richer topology compared to monolayers. Research on layered 2D ferromagnets is exceptionally active due to their intriguing stacking-dependent magnetic behavior and controllable ground states [56][57][58][59][60][61]. CrI 3 bilayers were extensively studied in this context. Their stacking-dependent interlayer coupling leads to ferromagnetic (FM) and layered antiferromagnetic (LAFM) spin configurations in AB and monoclinic stackings, respectively. Similar stacking-dependent magnetism was recently reported in CrBr 3 bilayers [62,63], with LAFM state in the R 33 stacking configuration. AB and monoclinic CrBr 3 bilayers, however, favor FM interlayer coupling. Controlled switching between LAFM and FM ground states is also possible in these quantum materials by various experimental techniques [7,59,60,64,65]. For illustration, the FM and LAFM ground states are presented schematically in figures 1(a) and (b), respectively. The monoclinic and R 33 stackings are illustrated in supplementary figures 1(a) and (b) (https://stacks.iop.org/NJP/23/053022/ mmedia), respectively.
The SOCs, stacking-dependent magnetism, and controllable magnetic ground states in 2D bilayer magnets constitute novel degrees of freedom that can induce a rich and controllable magnon topology. Advances in manipulating the magnon Hall conductivity, edge magnons, and valley transport in topological 2D bilayer magnets can lead to a new generation of magnonic and spintronic devices [66][67][68]. Magnon topology in AB-stacked honeycomb FM and LAFM bilayers was analyzed in reference [31], adopting a Heisenberg model with DMI. In both configurations, the magnon spectrum was found topological (Zeeman effect is necessary in the LAFM case) with a single phase. In a recent study [39], we have extended the analysis on AB-stacked honeycomb FM bilayers, including layer-dependent electrostatic doping [65] (ED). Five distinct topological phases were predicted as a result of the interplay between the model parameters. Moreover, recent theoretical studies [69,70] initiated the discussion on topological magnon valley transport in AB-stacked FM and LAFM bilayers. Nevertheless, a comprehensive investigation on the magnon topology and valley transport in honeycomb collinear FM and LAFM bilayers is still missing.
An essential step toward low-dimensional topological magnonics is to understand the implications of quantum interactions, symmetry breaking scenarios, magnetic ground states, and stacking orders on the magnonic topology and transport in 2D layered magnets. The present study serves this ambitious goal by analyzing the magnon topology in an extended group of honeycomb bilayer models with collinear magnetic order. The models are generated from all relevant combinations of Kitaev and/or DM interactions, inversion symmetry breaking terms, collinear FM and LAFM states, and standard stacking orders (AB, monoclinic, and R 33 , in particular). The comprehensive study reports several unconventional findings and elucidates the magnon topology's strong dependence on the stacking order, magnetic ground state, and SOC type. Finite Hall and Nernst conductivities are proved possible in FM and LAFM phases with zero Chern numbers. LAFM bilayer models reveal anomalous bandgap closures at unconventional nodal points (UNPs) formed away from the high symmetry points. These play a primary role in shaping the topological phase diagrams (TPDs) for LAFM bilayers. We further discuss opportunities to realize magnon valley transport in all models and suggest routes to determine the SOC's nature in several models based on the Hall and Nernst conductivities' profiles. The study aligns with current efforts to manipulate topological magnons in quantum 2D magnets toward their implementation in low-dimensional spintronics and magnonics.   Closes the gap between 2 and 3 at −K Cyan Closes the gap between 2 and 3 at UNPs away from the high symmetry points Before proceeding, we will set conventions and definitions adopted throughout the paper. The strength of the Kitaev interaction, DMI, and inversion symmetry breaking term will be denoted K, D, and U respectively. The bold notation ±K is used for the Brillouin zone (BZ) corners (or valleys). We break the inversion symmetry by introducing ED in the FM case and Zeeman effect due to an external magnetic field in the LAFM case. By default, all models include intralayer and interlayer Heisenberg exchange interactions. The short labels for our models are listed in table 1. The four magnonic bands are termed 1 , . . . , 4 in ascending energy order, and the color code for the corresponding band closure manifolds is defined in table 2. Nonadiabatic paths represented by black-dashed curves are introduced in selected TPDs. These paths help investigate the nonadiabatic evolution of the thermal Hall and Nernst conductivities. The conductivities are temperature-dependent and reach saturation at sufficiently high temperatures. Our figures present the maximum (saturation) values of these conductivities. Of note, the numerical results presented in the study are representative samples of our extensive numerical analysis. All stated conclusions are general and independent of the specific parameters used in the plots.

Magnonic Hamiltonian
The real-space Hamiltonian describing the various interaction in the collinear FM or LAFM honeycomb bilayer can be expressed as (1) The coefficients J, K, J ⊥ , and D mn determine the strength of the intralayer exchange, Kitaev interaction, interlayer exchange, and DMI, respectively. A stands for the easy axis anisotropy parameter. The inversion symmetry breaking parameter U l will be defined shortly. The vector S li denotes the spin on a site i in layer l (l = 1, 2 for the bottom and top layers, respectively). In the FM state, the spins in both layers can be expressed as S = S x , S y , S z . The LAFM configuration follows from the FM case by a π rotation of the top layer along the y-axis. Consequently, the Hamiltonian in equation (1) stays valid for the LAFM bilayers, with the spins in layer 2 expressed as S = −S x , S y , −S z .
The first, second, and third terms correspond to the nearest neighbor (NN) intralayer exchange, intralayer Kitaev, and interlayer exchange interactions, respectively. The indices i and j in these three terms are summed over the NN. The Kitaev contribution in equation (1) deserves further illustration as follows. For an A spin at site i, the Kitaev unit vectors areγ i1 = − The fourth term in H corresponds to the DM interaction, and the indices m and n are summed over the six next NNs. These are situated at relative position vectors, figure 1(e)). We have introduced the vector S D ln = S lny , −S lnx , 0 which transforms the DMI contribution to a scalar-product rather than a cross-product [39]. The coefficient D mn = ±D, with the sign determined in the conventional way from the local geometry of the honeycomb lattice [30].
The last term in H accounts for the easy-axis magnetic anisotropy, while the fifth term is introduced to break the inversion symmetry. For LAFM bilayers, the fifth term is induced by an external magnetic field h = Uẑ, setting U l = ±U for layers 1 and 2, respectively. For the FM bilayer case, we adopt the approach suggested in reference [69] to break the inversion symmetry via layer-dependent ED. In particular, the ED technique can efficiently tune the localized moments in each layer [65,71], generating opposite potentials (±U) when the two monolayers are doped equally with opposite charges. An alternative approach is to induce layer-dependent magneto-crystalline anisotropy in the FM bilayer. This can be achieved when the FM bilayer is sandwiched between a substrate and a superstrate composed of different materials. The sandwiched structure ensures a distinct atomic environment for each layer and allows independent manipulation of the monolayers' magnetic properties, including the magneto-crystalline anisotropy In supplementary notes 1 and 2, we develop the semi-classical linear spin-wave approach [72][73][74][75][76][77] and derive the FM and LAFM momentum p space Hamiltonians for bilayers stacked in the AB, monoclinic, and R 33 configurations. The Holstein-Primakov method yields identical results. To our knowledge, the semi-classical approach is not yet developed for FM and LAFM bilayers with Kitaev interaction, which made us prefer it over the Holstein-Primakov formalism. The Hamiltonians are derived as 8 × 8 matrices and produce four physical bands. We will use the notation [C 4 , C 3 , C 2 , C 1 ] to group the Chern numbers of the bands [ 4 , 3 , 2 , 1 ], respectively. Details on the Berry curvatures and Chern numbers calculation are presented in the methods section.
Before proceeding, we fix the value of the easy-axis anisotropy parameter A. This parameter is essential to stabilize the magnetic order in real 2D magnets. Nevertheless, as long as magnons are concerned, the effect of A is a uniform upward shift in the magnon spectra. A consequently lifts the degeneracy between positive and negative solutions, which is imperative for topological calculations. The specific value of A is, however, irrelevant to the topological results. We will fix the easy axis anisotropy at A= 0.25J in the numerical calculations, chosen arbitrarily as per the above argument.

Topological analysis in AB-stacked FM bilayers
The D model for AB-stacked FM bilayers was studied in reference [31], revealing a single topological phase with Chern numbers [0, −2, 0, 2]. The rich magnon topology in the FM D + U model was discussed in our recent work [39]. The model hosts five topological phases with Chern The corresponding topological phase transitions are accompanied by bandgaps closure at the ±K BZ corners exclusively. The U model for these bilayers was extensively analyzed in references [69,70], reporting valley-polarized magnons in a remarkable analogy with electronic 2D materials. These three models are hence excluded here.
We start with the K model for AB-stacked FM bilayers. The magnon spectrum in the present case is gapped and topological for any nonzero value of K, with Chern numbers [0, −2, 0, 2]. This is the same topological phase encountered in the D model [31]. An illustrative example of the gapped four-band magnon spectrum is presented in figure 2(a) for J ⊥ = 0.3J and K = 2J. The bands are plotted along the high symmetry axes KΓ, ΓM, and MK of the conventional honeycomb BZ (not shown). Including the DMI (K + D model) augments the bandgaps without affecting the above conclusions regarding the magnon topology.
In the K + U model, the magnon band spectrum is gapped except at specific 3D manifolds in the (K, U, J ⊥ ) parametric space. These manifolds close the gaps at the BZ corners (±K) and lead to topological phase transitions. To illustrate, the TPD is presented in figure 2 [39]. Consequently, the DM and Kitaev interactions in AB-stacked FM bilayers cannot be easily differentiated in terms of their topological consequences. Region VI, however, is gapped but topologically trivial with zero Chern numbers.
Interestingly, the Berry curvatures in VI do not vanish and are peaked at the BZ corners (figures 2(c) and (d)) with opposite signs near ±K. This sign behavior is absent in the topological phases. The Berry curvatures' profile in VI is similar (but not identical) to the U model, which was proved to support topological valley transport [69,70]. We will shortly illustrate the consequences of the nonzero Berry curvatures on the standard and valley Hall conductivities in VI.
The rich topology in the K + U model motivates the analysis of the topological conductivities' nonadiabatic evolution. The black-dashed path in figure 2  with weak ED (figure 1(c)). The region with weak ED might be particularly interesting for experimental studies. The Hall conductivity in this region is found to present an exotic profile, with significantly abrupt jumps and asymptotic-like behavior near the boundaries. Importantly, κ xy and α xy do not necessarily vanish in the trivial phase VI (figures 2(e) and (f), and supplementary figure 1(d)). For sizable Kitaev interactions, the absolute values of the Berry curvature's peaks differ at ±K (figure 2(c)), resulting in net Hall and Nernst conductivities. For weaker (but nonnegligible) Kitaev interactions, the peaks' absolute values are equal at ±K (figure 2(d)), in analogy with the U model [69,70]. Consequently, the standard conductivities κ xy and α xy vanish, while the valley Hall and Nernst conductivities become finite [69].
In the most general K + U + D model for AB-stacked FM bilayers, the parametric space is 4D, generated by the Hamiltonian parameters (K, U, J ⊥ , D). The DMI added to the Kitaev interaction is found to reconstruct the gaps closure manifolds without inducing new topological phases compared to the previous K + U model. Nevertheless, the DMI wipes out the non-topological phase rendering the model topological for any parametric combination away from the gap closure manifolds. Moreover, the behavior of κ xy and α xy in the K + U + D model does not show any relevant difference compared to the K + U model: the Before closing this section, we point out the important implications of the magnon spectrum's nonreciprocity ( i p = − i − p ) in models with both ED and SOCs. The independent gap closures at ±K (figure 2(b)) are a direct consequence of this bands' nonreciprocity. Moreover, in the K + U model, the bands' nonreciprocity plays an essential role in realizing the Hall and Nernst conductivities in the trivial phase VI. Bands' nonreciprocity holds in all models with sizable ED and SOCs, regardless of the magnetic ground state and stacking order.

Topological analysis in AB-stacked LAFM bilayers
In the AB-stacked LAFM configuration, a perpendicular external magnetic field is indispensable to lift the degeneracy and induce bandgaps. Models with U = 0 are hence irrelevant to the present study. Like the FM case, the bandgaps in AB-stacked LAFM models with DMI and/or Kitaev interactions can close at ±K, resulting in numerous topological phases. Nevertheless, the magnon spectra in LAFM bilayers have an essential aspect that is absent in FM bilayers. Namely, they display UNPs that close the bandgaps away from ±K, specifically between 2 and 3 bands. These play a crucial role and enrich the magnon topology in LAFM bilayers.
Reference [69] discussed the U model for AB-stacked LAFM bilayers. Here, we complement the study and highlight the role of the UNPs in this model. The magnon spectrum is gapped at ±K for any nonzero value of U. Nevertheless, the interplay between J ⊥ and U can close the gap between 2 and 3 at the UNPs. Examples of the UNPs and the corresponding bandgap closures are presented in the left panel of figure 3(a). In the right panel, we plot bands 2 and 3 over the shaded region centered at the +K corner of the BZ. The bandgap between 2 and 3 reaches minimal values along the contour highlighted in red. The bandgap closes at the UNPs near the intersection between the red contour and the high symmetry axes KΓ and MK (left panel of figure 3(a)). This discussion is valid for all LAFM models presenting UNPs. Figure 3(b) plots the minimal bandgap in the U model as a function of the normalized external field (U/JS) and interlayer interaction (J ⊥ /J). The contour plot reveals a gapped phase in weakly coupled LAFM bilayers. The Berry curvatures are nonzero in the gapped phase, with antisymmetric peaks near the ±K valleys. Hence, the standard Hall and Nernst conductivities vanish, but valley-polarized magnons can be generated in the gapped region, as pointed out in reference [69].
The D + U model in AB-stacked LAFM bilayers was first studied in reference [31], excluding the important consequences of the UNPs. Figure 3(c) presents the minimal bandgap's contour plot in the dU space for an exaggerated value of the interlayer exchange, J ⊥ = −0.5J. The figure equally serves as a TPD. The TPD preserves its main features for more realistic (weaker) J ⊥ in van der Waals (vdWs) magnets. The numerical results reveal a topological phase I, an ungapped phase II with gap closures at UNPs, and a gapped but trivial phase III with zero Chern numbers. Phase III disappears for weaker interlayer exchange (|J ⊥ | < 0.2J as a rough estimation).
The ±K gap closure manifolds still exist ( figure 3(c)), notably in phase II, but do not contribute to the magnon topology. The green manifold in figure 3(c) coincides with the pink one, simultaneously closing the ( 1 , 2 ) and ( 3 , 4 ) gaps at +K (see the color code in table 2). Illustrative examples of the gap closure scenarios are presented in figures 3(d) and (e) at the intersection between the green, pink, and orange manifolds.
The Berry curvatures (not shown) are nonzero in the trivial phase III and have opposite-sign peaks near ±K. Nevertheless, the absolute values of the peaks are significantly different at ±K. Hence, the trivial phase can support standard Hall and Nernst conductivities, but magnon valley transport seems unlikely in this phase.
Replacing the DMI by Kitaev interaction enriches the magnon topology drastically. Unlike the D + U model, both conventional and unconventional gap closure manifolds conspire to shape the topology in the K + U model. The TPD is presented in figure 4 The drastic difference between the magnon topologies in the D + U model (1 topological phase) and the K + U model (6 topological phases) is remarkable. The magnon topology in the K + U model is tunable via experimentally controllable parameters, like the external magnetic field and interlayer interactions [59]. While topological phase transitions are unlikely in the D + U model, they can be experimentally realized in the K + U counterpart. This can help determine the nature of the SOC in AB-stacked 2D magnets with LAFM ground state.
Furthermore, the rich magnon topology in the K + U model induces novelties in the profiles of κ xy and α xy . These are plotted in figures 4(c) and (d) along the nonadiabatic directional path of figure 4(a), specified by U = 0.4 cos 4 (u) , K = 2 sin(u) + cos 2 (u) + 1, − π 2 < u < π 2 . The topological phase transitions induce abrupt jumps and can even lead to asymptotic-like behaviors at some boundaries. Experimental observation of any sharp variation in κ xy and α xy shall constitute an unambiguous signature of the Kitaev interaction since topological phase transitions are absent in the D + U model. More importantly, the topological phase transitions in the K + U model can reverse the signs of κ xy and α xy , revealing a unique feature of AB-stacked LAFM Kitaev magnets. The sign reversal is best manifested in the Hall conductivity, showing two asymptotic-like behaviors of opposite signs at the (V-IV) and (IV-III) boundaries. As will be demonstrated shortly, these signatures are absent in samples with relevant values of the DMI. Consequently, experimental observation of conductivity sign reversal eliminates the chances for sizable DMI and uncovers the SOC's exclusive Kitaev nature.
Before proceeding to the K + U + D model, we note that the signs of κ xy and α xy in the D + U model's single topological phase can be reversed only by reversing the direction of the external magnetic field [31]. However, the observed sign reversal in the K + U model with preserved magnetic field direction is of fundamental nature and related to the topological phase transitions. This sign reversal shall be detectable in experiments with appropriate control on the external magnetic field's strength and/or the interlayer exchange.
The phase diagram in the K + U + D model for AB-stacked LAFM bilayers is found to be significantly different compared to the K + U model. The DMI enhances the gaps between 2 and 3 , which can only close in bilayers with tiny DMI and exaggerated interlayer exchange. The reduced number of band closure manifolds decreases the number of topological phases, leaving only 2 of them in bilayers with sizable DMI and/or weak vdWs interlayer exchange. The TPD is presented in figure 4(b) for J ⊥ = −0.2J and D = 0.05J. Similar to the FM case, the K + U + D model does not allow topologically trivial phases. The Chern for phases I and II, respectively. These match the first two topological phases reported in the previous K + U model. Topological phase transitions between I and II cannot reverse the signs of the Hall and Nernst conductivities. Consequently, sign reversal of κ xy and α xy is not allowed in AB-stacked LAFM bilayers with relevant DMI and/or weak (vdWs) interlayer exchange. We hence restate our conclusion that sign reversals unambiguously reveal the exclusive Kitaev nature of the SOC.

Topological analysis in monoclinic FM and LAFM bilayers
The previous section elucidated the significant topological difference between FM and LAFM ground states in AB-stacked bilayers. In this section, the stacking-dependent topology is exposed.
We start with monoclinic FM bilayer models. The U model is very similar to the AB-stacked case and supports topological magnon valley transport for any U. The K, D, and K + D models are gapped with a single topological phase [1, 1, −1, −1], which is different from the [0, −2, 0, 2] phase in the AB-stacking. The magnonic Hamiltonian hence evolves nonadiabatically between the AB and the monoclinic stacking. The remaining models (K + U, D + U, and K + U + D) are topologically identical to the K and D models. The rich magnon topology reported in the AB-stacked FM bilayers is indeed absent in the monoclinic case.
We arrive at LAFM monoclinic bilayers with some novel results. Unlike the AB-stacked case, the monoclinic U model is not gapped, and magnon valley transport is not possible. The bandgaps in the U model close at ±K and/or the UNPs, depending on the value of J ⊥ and U. The K and D models are topologically trivial, where 1 and 2 bands are degenerate throughout the BZ (similarly for 3 and 4 ). Next, the TPD for the D + U model is plotted in figure 5(a) with J ⊥ = −0.1J. Phase I is topological, with Chern numbers [−1, 1, 1, −1], while phase II is ungapped. The bandgaps close at UNPs in phase II. Phase II can be realized experimentally via ED. This ungapped phase is absent in the K + U and K + U + D models, presenting phase I only. Consequently, the realization of an ungapped phase in LAFM monoclinic bilayers reveals the exclusive DMI nature of the SOC. We also note that topological magnon valley transport is absent in all monoclinic LAFM models.

Topological analysis in R 33 stacked FM and LAFM bilayers
To our knowledge, magnon topology in R 33 stacked bilayers was not discussed in any previous work. The results for the FM U model are very similar to the AB and monoclinic cases. For FM bilayers with SOC, the R 33 configuration presents new results that are absent in the AB and monoclinic cases. Figure 5 −1, 1, −1]. The nonadiabatic behaviors of κ xy and α xy (not shown) in all R 33 models with SOCs are comparable to figures 2(e) and (f). The topology of the K + U, D + U, and K + U + D models can be tuned by U, and the significant topological difference between these models might be useful to explore the nature of the SOC. Finally, the TPDs and related conclusions in R 33 stacked LAFM models are very similar to their monoclinic counterparts.

Discussion
The stacking-dependent magnetism, SOCs, and electric control of magnetic order in 2D bilayer magnets have recently received exceptional attention given their potentials in spintronics. Nevertheless, the implications of these degrees of freedom on the magnon topology remain widely unexplored. The present study analyzed the magnonic topology and valley-polarization in 39 honeycomb bilayer models that span all relevant combinations of SOCs, magnetic configurations, stacking orders, and inversion symmetry breaking. The study is kept general, presenting numerical results over extensive ranges of the fundamental parameters to cover the broadest spectrum of vdWs magnets. The main conclusions are summarized in tables 3 and 4, revealing the remarkable topological dependence on all degrees of freedom included in our study.
The work is further motivated by the open question regarding the nature of SOCs underlying 2D magnets [24]. We show that the standard Hall and Nernst conductivities carry signatures that can be matched with the type of SOC in several bilayer models. Indeed, magnetic materials' topological nature is best explored through their Hall and Nernst responses. Nevertheless, our study concludes that these conductivities are strongly related to the SOC's strength and bands' nonreciprocity, and they can survive in bilayers with zero Chern numbers.
Understanding the relation between bandgap closures and the structure of the TPDs is another primary objective. We revealed fundamental differences between the TPDs in FM and LAFM bilayers. The TPDs in FM bilayers present conventional bandgap closure manifolds associated with Dirac cones located at the ±K BZ corners. However, in LAFM bilayers, the TPDs are tailored primarily by the UNPs, introducing a novel bandgap closure scenario away from the high symmetry points. These contribute additional boundaries to the TPDs and enrich the topology.
We end with a brief discussion on the advantages of 2D magnetic bilayers over their monolayer version. Bilayers present additional degrees of freedom, detailed in our study, which can help tune their magnonic bands and topological response. Monolayers are restricted to a single topological phase in the presence of Kitaev and/or DM interactions. Consequently, monolayers' topology cannot be tuned, and the information extracted from their topological response is limited. Bilayers are hence superior from fundamental and applied perspectives.

Methods
The Berry curvatures and Chern numbers are determined using the approach developed in reference [78]. We first discretize the BZ with steps δp in the x and y directions and numerically calculate the phase variations generated by infinitesimal displacements in the BZ, Here, | i denotes the eigenstate for the band i. The Wilson loop W i ( p) can then be deduced as