Bose-Einstein condensate of Dirac magnons: pumping and collective modes

We explore the formation and collective modes of Bose-Einstein condensate of Dirac magnons (Dirac BEC). While we focus on two-dimensional Dirac magnons, an employed approach is general and could be used to describe Bose-Einstein condensates with linear quasiparticle spectrum in various systems. By using a phenomenological multicomponent model of pumped boson population together with bosons residing in Dirac nodes, the formation and time evolution of condensates at Dirac points is investigated. The condensate coherence and its multicomponent nature are manifested in the Rabi oscillations whose period is determined by the gap in the spin wave spectrum. A Dirac nature of the condensates could be also probed by the spectrum of collective modes. It is shown that the Haldane gap provides an efficient means to tune between the gapped and gapless collective modes as well as controls their stability.

Magnon Dirac materials [27] such as van der Waals crystals Cr 2 Ge 2 Te 6 [28] and the transition metal trihalides family AX 3 (where A = Cr and X = F,Cl,Br,I), e.g., CrI 3 [23][24][25]29] received attention recently. These materials are layered ferromagnets with typical ferromagnetic transition in the temperature range of 10-50 K. Transition metal trihalides were known for a long time [30,31] and can be viewed as a magnetic analog of ABC stacked graphite. In three dimensions (3D), a magnonic Weyl state was predicted in pyrochlores [32][33][34]. Linear crossing in the magnon spectrum was indeed observed in the antiferromagnet Cu 3 TeO 6 [35] and the multiferroic ferrimagnet Cu 2 OSeO 3 [36]. It is known that a Dirac-like magnon spectrum is generated for localized magnetic moments on a 2D honeycomb lattice [23,27]. Wide range of properties within this emergent class of materials: interaction effects [23,37] that renormalize the slope of the dispersion, topological properties [36,[38][39][40][41][42][43], strain effects [44][45][46] and non-Hermitian dynamics [47] were investigated observed in these materials. While some properties of fermionic and bosonic Dirac materials are similar, there are, of course, crucial differences. The key is the difference in particle statistics: there is no Pauli principle for bosons. Therefore, no Fermi surface exists in bosonic Dirac materials. At low temperatures the chemical potential for number conserving bosons is small and the only relevant excitations are near the center of the Brillouin zone at the minimum of dispersion.
For the excitations like magnons, the equilibrium chemical potential is zero. Magnons are excitations that are absent in the ground state in equilibrium. Thus no BEC of magnons is possible in equilibrium. Nevertheless, one can ask if it possible to create the coherent states or even Bose-Einstein condensate of Dirac bosons and what would be the proper description? We propose that the answer to this question is positive and suggest a realization of coherent Dirac bosons and condensates by creating steady-state nonequilibrium population of magnons at the Dirac nodes.

II. MODEL
In this section, we present a phenomenological approach to Dirac BECs. The approach is similar to that for conventional BECs [51] albeit incorporates the pseudospin and a linear Dirac dispersion relation. We start with the following minimal ansatz for the free energy density (for a schematic derivation in the case of magnons, see Appendix A): g a n a,ζ n a,ζ + 1 2 g b n b,ζ n b,ζ + g ab n a,ζ n b,ζ + 1 2 ζ=± u ab ψ † a,ζ ψ † a,ζ ψ b,ζ ψ b,ζ + u * ab ψ † b,ζ ψ † b,ζ ψ a,ζ ψ a,ζ . (1) The four-component wave function of Dirac BEC is where ζ = ± denotes the K or K ′ valley in 2D or chirality in 3D, ψ a,ζ and ψ b,ζ are the condensate wave functions for the different pseudospins, i.e., A and B sublattices for 2D hexagonal lattices, σ is the vector of the Pauli matrices in the pseudospin space, τ z is the Pauli matrix in the valley space, n i,ζ = ψ † i,ζ ψ i,ζ with i = a, b is the density for a certain pseudospin, n tot = ζ=± (n a,ζ + n b,ζ ) is the total condensate density, c 0 is the position of the Dirac nodes in energy, and µ is the effective chemical potential of condensates. Finally, the term quantified by ∆ denotes the Haldane gap in 2D. (Note that term ∝ ∆ plays a different role in 3D where it shifts the Weyl nodes of opposite chiralities in momentum space and is called the chiral shift [72]. In what follows, however, we focus only on 2D case.) Furthermore, we note that the part of the free energy linear in the density is similar to that for other Dirac (quasi-)particles. In particular, the part with spatial derivatives is similar to that for Dirac fermions. Moreover, there are two pseudospin or sublattice degrees of freedom. The fact that kinetic energy terms and Hilbert space of the excitations is very similar to the case of known Dirac excitations allows us to dub the condensate as the Dirac BEC.
In addition to the conventional for Dirac systems linear terms, we included also the nonlinear interaction terms quantified by constants g a , g b , and g ab . While the former two correspond to the interaction strength of boson densities of the same pseudospin, the latter describes the inter-pseudospin interaction. In addition, the Josephson coupling term quantified by u ab is considered. This term is reminiscent to the spin-orbital coupling term in conventional spin-1 BECs [73]. Cross-condensate Josephson terms describe the internal "Josephson effect" where a pair of magnons of type a converts into a pair of magnons of type b. As we will show below, the Josephson coupling term mixes the phases of the condensate wave functions corresponding to different pseudospins and might significantly affect the dynamics of the condensates.
Finally, we notice that several other terms are possible in the free energy, e.g., terms corresponding to the inter-valley mixing ∝ n a,−ζ n a,ζ . For the sake of simplicity, we consider the effect of only a few parameters on Dirac condensates. A more detailed discussion is left for the future investigation.
The 2D Dirac dispersion relation (see Appendix A for the discussion of the Dirac spin wave dispersion relation) and BECs corresponding to different valleys and pseudospins are shown schematically in Fig. 1. Here red shaded regions correspond to the accumulated magnons (or other Dirac bosonic quasiparticles) and spikes represent condensates. Note that unlike fermionic systems, energy is calculated with respect to the bottom of the bands. Two Dirac points are situated at k = ±k D in momentum space and at ω a = ω b = c 0 in energy. As one can see from Fig. 1(b), the presence of the ∆ term opens the gap in the spectrum.
By varying the free energy (1) with respect to the fields ψ † a,ζ and ψ † b,ζ , we derive the following Gross-Pitaevskii equations for Dirac BECs at A and B sublattices: These are nonlinear Dirac equations where bosonic fields have four degrees of freedom: two pseudospins and two valley or chirality indices. It is interesting to notice the nontrivial topology of the multicomponent BEC. If all coupling terms are ignored, the corresponding state would correspond to the U (1) a+ × U (1) a− × U (1) b+ × U (1) b− symmetry. Therefore, one would expect the rich phase diagram with multiple phases that could emerge in Dirac BECs, where relative phases of these symmetries can be broken. For simplicity and to make an explicit case, we will assume that relative oscillations of the condensates in both valleys are suppressed an there are only two flavor symmetries, which are explicitly broken to U (1) a × U (1) b . In other words, we consider a symmetric response of the valleys. Therefore, we fix ζ = + and omit the valley index when it is not necessary. Nevertheless, we mention some non-trivial phases of the non-interacting BECs in Appendix B.
Unlike fermions, where there is a Fermi surface that can be tuned to intersect Dirac points, bosons occupy the lowest energy state. Therefore, in order to investigate the properties of the Dirac BEC, one first needs to achieve a sufficient quasiparticle density at the Dirac points. As we discuss in the next section, Dirac Bose-Einstein condensation can be achieved in a steady state created by pumping.

III. GENERATION OF DIRAC BOSE-EINSTEIN CONDENSATES VIA PUMPING
In this section, we consider pumping and time evolution of Dirac BECs. Since the pumping scheme depends on the nature of Bose-Einstein condensates and, in general, is different for polaritons, magnons, Cooper pairs, etc., we focus on the case of 2D Dirac magnons. (For pumping and dynamical instabilities in fermionic Dirac matter, see, e.g., Ref. [48].) As we discussed in the introduction, Dirac magnons can be realized in, e.g., CrI 3 , see Refs. [23][24][25]. Note that the Haldane gap for Dirac magnons can be generated via the pseudo-Zeeman effect (for a discussion of the strain-induced pseudo-Zeeman effect in graphene, see Ref. [74]). In addition, the gap in the Dirac spectrum could appear due to the Kitaev interaction. For a recent study in CrI 3 , see, e.g., Ref. [25].

A. Three-component model
To describe the pumped Dirac BEC, we introduce a phenomenological three-component model (or five-component if the valley index is taken into account). In addition to the magnon populations at the A and B sublattices for the K point, we include the pumped magnon population with the energy ω 3 (see also Fig. 1). The latter population could be achieved by the parametric pumping similarly to the case of conventional magnetic materials [71]. The system reads as ∂ t n 3 + Γ 3 n 3 = Γ 3P (t) − ζ=± (P 3a Γ a n a,ζ + P 3b Γ b n b,ζ ) n 3 .
Compared to Eqs. (3) and (4), we have added the dissipation terms iΓ a ψ a and iΓ b ψ b that correspond to the magnonphonon interactions. The terms with P a3 and P b3 describe to the inflow of magnons from the pumped population and act as the source terms for Dirac magnons. On the other hand, the terms with P 3a and P 3b correspond to the depletion of the pumped magnon population. The latter is described by the rate equation (7) for the magnon density n 3 . As in the case of the magnon BEC, Γ 3 denotes the decay rate. Finally, the term Γ 3P (t) describes the direct inflow of magnons caused by the pump. For the sake of simplicity, we consider a simple step-like profile for the pump power Here the pump starts at t = t Begin and ends at t = t End . Further, P max is power of the pumping field and θ(x) is the step function. The square-root dependence on the pump power is in agreement with the experimental data for YIG [65,75] (see also Refs. [76,77] for the theoretical discussion). While only a simple step-like profile in Eq. (8) is considered in this work, the time profiles of magnon BECs can be modulated by applying more complicated timedependent pump.

B. Time evolution of Dirac magnons
In this subsection, we present the numerical results for the pumped magnon populations and discuss the role of the Haldane gap and the Josephson coupling terms. Unless otherwise stated, we use the following set of the model parameters: Here t 0 is the characteristic interaction timescale. Since no experimental data is present for the Dirac magnon BEC, we use model parameters. In addition, we note that, in general, P i3 = P 3i with i = a, b. For example, a strong disparity between the scattering efficiency of gaseous population to magnons at the bottom of the energy dispersion and vice versa is observed for YIG in Ref. [68]. Further, because of the space separation between the lattice site, the inter-sublattice interaction constant g ab is expected to be much weaker than the intra-sublattice ones g a and g b . Therefore, we focus on the case |g ab | ≪ g a , g b . The time profiles of magnon densities are shown in Figs. 2 and 3 for several values of the pump power. As one can see, the pumped magnon population n 3 quickly reaches its maximum value determined by the interplay of the pumping power and the decay rate. Magnon populations at the Dirac points occur later in time as well as require pump power to exceed a certain threshold (see Fig. 3). It is interesting to note that the densities of the magnon BECs reach their maximum values at the onset of the profile. In addition, since we chose equal inflow and outflow terms in Eqs. (5) through (7), there is a noticeable depletion of the pumped population. After turning off the pump, both condensates and pumped population decay during the time τ ∼ 1/Γ. It is notable, also, that the total number of magnons is not conserved and shows a spike when the BECs start to form. This could be qualitatively explained by the fact that the Dirac points accumulate more magnons during the rise phase than they can support in the steady state. These magnons are released during the equilibration phase.
Further, we investigate the effect of the Haldane gap ∆. The corresponding results are shown in Fig. 2 at a few values of the gap for P max = 0.5 P 0 . As one can see, the gap leads to different densities of magnon BECs at the A and B sublattices. Therefore, it can be used to effectively control the pseudospin distribution of Dirac magnon BECs.
Finally, let us study the effect of the Josephson coupling term u ab on the time evolution of magnon densities. The corresponding results for a few values of u ab are shown in Fig. 4. As one can see, the Josephson coupling term ∝ u ab leads to the oscillations of the Dirac BEC. Oscillations of n a and n b always occur in antiphase and the pumped magnon population remains constant. The observed oscillations are, in fact, the Rabi oscillations between magnons BECs with different pseudospins. To elucidate the nature of the inter-condensate oscillation and to show that they have the same scaling as the Rabi oscillation, we plot the corresponding period as a function of ∆ in Fig. 5. The period indeed scales as 1/∆ for large values of ∆. The case of small ∆ is more complicated, however. For example, the period becomes a nonmonotonic function and depends, e.g., on the value of u ab and other interaction parameters. A nontrivial dynamics becomes manifested at ∆ t −1 0 . Thus, the gap not only induces different magnon densities at the A and B sublattices, it can activate the Rabi oscillations between the condensates. We believe that the experimental observation of the Rabi oscillation of the magnon densities could serve as a definitive proof of their coherence. Such observations can be done via the Brillouin light scattering technique used for conventional magnons in YIG. The dynamics of the Dirac magnon populations can be presented by introducing the vector in the density space (n a , n b , n 3 ) /n tot . The tip of such vector is always confined to the one-eight of a unit sphere. The evolution of densities then resembles the dynamics of two-level quantum system whose pure states correspond to the surface of the Bloch sphere. The animation of dynamics of the Dirac magnon populations is given in the Supplemental Material. The snapshot after the pump is turned off (t = 300 t 0 ) is shown in Fig. 6.

IV. COLLECTIVE MODES
In the previous section, we considered pumping of uniformly distributed magnons. To provide a deeper insight into the Dirac nature of magnon condensate and uncover properties related to the linear spectrum, one should consider coordinate-dependent probes. Collective modes are, perhaps, among the most well-known, powerful, and versatile of them.
Before we proceed with more detailed analysis, let us briefly discuss the structure of collective modes. There are two condensate densities n a and n b . Hence we would expect two sound modes where the coupling is independent of We see a significant scatter in the period due to nonlinear dynamic at ∆ t −1 0 . In this case, other interaction parameters also become important. the relative phase. Two linear modes will remain linear with renormalized velocities once the interaction is turned on. This is what we indeed find in the regime n a = 0 and n b = 0. However, we also find the state where the condensate breaks the pseudospin symmetry and one of the components vanishes in the ground (or, more precisely, quasi-ground) state, e.g., n b = 0. In that case, we find only one gapless mode corresponding to the phase mode of the condensate. The second mode is gaped and corresponds to the Higgs or amplitude mode, where n a and n b populations oscillate.
In order to investigate the spectrum of collective excitations, we employ the Bogolyubov approach [49][50][51] for the magnon BEC. In this case, one looks for the solution as ψ a = ψ a,0 (r) + e −iµt u a e −iωt+ik·r + v * a e iω * t−ik·r , FIG. 6: The trajectory of the systems in the density space (na, n b , n3) /ntot is shown by the red line. Since the population densities are normalized to the total density, the dynamics is always constrained to the one-eight of a unit sphere. When the pump is turned on, only the pumped population n3 contributes to the total density. Later, the condensates develop and the density vector start to oscillate closer to the equator. Finally, the condensates decay and the system returns to its initial state.
Here the values of parameters given in Eqs. (9) and (10)  Here ψ a,0 (r) and ψ b,0 (r) denote ground-state solutions, which are, in general, nonuniform. They satisfy the timeindependent systems of the Gross-Pitaevskii equations. The perturbations are described by u a,b and v * a,b terms. Finally, ω is the frequency and k is the wave vector of collective modes.

A. Ground state
Let us start with the ground-state solutions ψ a,0 (r) and ψ b,0 (r). For simplicity, we consider a uniform ground state. By using Eqs. (11) and (12), the Gross-Pitaevskii equations (3) and (4) lead to 0 = (c 0 − µ + ∆) ψ a,0 + g a n a ψ a,0 + g ab n b ψ a,0 + u ab ψ * a,0 (ψ b,0 ) 2 , It is convenient to separate the absolute value and the phase of the ground-state wave functions, i.e., ψ a,0 = √ n a e iθa and ψ b,0 = √ n b e iθ b . As follows from Eqs. (13) and (14), the phase of the ground-state solutions is not fixed for u ab = 0. Therefore, the state has U (1) a × U (1) b symmetry. If the Josephson coupling term is nonzero, we have where u ab = |u ab |e iθu ab . Further, we find an effective chemical potential that allows for a solution of the Gross-Pitaevskii equations. We obtain µ = c 0 + ∆ + g a n a + g ab n b + u ab n b , Note that even and odd values of l correspond to positive and negative values of u ab in the equation above. The system (16) and (17) has a solution for (One can fix n b and solve for n a , which is physically equivalent). The effective chemical potential for the density (18) reads µ = c 0 + ∆ + g a n a + (g ab + u ab ) 2∆ + (g a − g ab − u ab ) n a g b − g ab − u ab .
There is also an additional solution for Eqs. (13) and (14), where n b = 0 and µ = c 0 + ∆ + g a n a .
We note that solution (18) exists as long as n b > 0. Otherwise, only the ground state with n b = 0 and µ given in Eq. (20) is possible. Therefore, even for simple uniform and static equations for the Dirac BEC, we find two types of solutions: one with nonvanishing magnon densities n a = 0 and n b = 0 and the solution that breaks the pseudospin symmetry with n a = 0 but n b = 0. The latter solution is an example of nematic instability similar to nematic states in fermion superconductors [78] and in Bose nematics (see, e.g., Ref. [79]).

B. Free energy
In order to clarify which of the two possible solutions for the steady-state homogenous Gross-Pitaevskii equations, discussed in Sec. IV A, is the true ground state, we calculate the free energy of the system. By using Eq. (1), we derive the following expression for n b given in Eq. (18): In the case n b = 0, we have a simpler expression In addition to comparing the free energies, one should check whether n b remains positive for the state with n b = 0. We present the phase diagram of the system for a few combinations of parameters in Fig. 7. Unless otherwise stated, we use the values of parameters given in Eqs. (9) through (10). The state with the lower free energy is a true (quasi-)ground state. As one can see, the phase diagram is nontrivial where states with both n b = 0 and n b = 0 are possible. In general, however, the ground state with n b = 0 is favored for repulsive interactions (g a > 0 and g b > 0). Furthermore, it is clear that the Haldane gap parameter could be used to tune the ground state (see Figs. 7(c) and 7(d)).

C. Dispersion relations
In this section, we use the solutions found in Sec. IV A and linearize the Gross-Pitaevskii equations (3) and (4). By using relations (11) and (12), one obtains the following equations for u a,b and v a,b : −ωv a = (c 0 + ∆ − µ) v a − v(k x + ik y )v b + 2g a,1 n a v a + g ab,1 n b v a + g ab,1 ψ * a,0 ψ b,0 v b + g a,1 (ψ * a,0 ) 2 u a + g ab,1 ψ * a,0 ψ As before, we omit the valley index ζ assuming that the valleys have the same dynamics. Note that complex conjugation was performed in Eqs. (25) and (26). Equating the determinant of the system (23) through (26) to zero, the spectrum of collective modes is determined. In what follows, we focus on the case u ab = 0. The presence of the Josephson coupling term complicates the expressions for the frequencies. For example, numerical calculation suggest that u ab enhances the imaginary parts of the frequencies in the case of the synphase condensates (θ a = θ b = θ u ab = 0). For a more detailed discussion and a few dispersion relations at u ab = 0, see Appendix C.
In the case of the nontrivial ground state solution given in Eq. (18), the frequencies are Note that both branches are gapless. In addition, we fix the phase of the ground state solutions as θ a = θ b = 0. Even at ∆ = 0, the dispersion relation (27) is rather interesting. For example, the mode is stable (i.e., there is no large imaginary part) for g 2 ab − g a g b (g ab − g a ) (g ab − g b ) < 0. In the case g a = g b , the stability criterion is simple g 2 a ≥ g 2 ab , which agrees with phase-separation criterion for conventional BECs [51]. For the ground state with n b = 0, one of the collective modes (e.g., ω + ) can be gaped. The full spectrum in this case is The gap between ω + and ω − branches reads Note that the gap is determined both by the Haldane gap term and interaction terms. In a simple case ∆ = 0, the modes at n b = 0 are stable at g 2 a ≥ g 2 ab . We present the spectra (27) and (28) in Figs. 8 and 9, respectively. It is notable that the rotational symmetry is spontaneously broken in the interacting BEC of magnons at n b = 0. Indeed, this is evident from the expression in Eq. (27), where an explicit dependence on k y is present. It can be checked that the anisotropy is set by the relative phase of the solutions ψ a,0 and ψ b,0 . For example, one could replace k y → k x for θ b − θ a = π/2. Furthermore, both collective modes are gapless. These modes can be identified with two Nambu-Goldstone modes related to the independent oscillations of the condensates at A and B sublattices, i.e., spontaneous breakdown of the symmetries ψ a,0 → e iθa ψ a,0 and ψ b,0 → e iθ b ψ b,0 .
The spectrum at n b = 0, on the other hand, contains one gapless and one gapped modes. Indeed, the presence of a gapless Nambu-Goldstone mode is guaranteed because U (1) a symmetry ψ a,0 → e iθa ψ a,0 is spontaneously broken in this case. The other mode can be gapped. It corresponds to the Higgs mode of the Dirac BEC. This mode is also known as the amplitude mode because it corresponds to the oscillations of the amplitude of the order parameter. On the other hand, the gapless Nambu-Goldstone mode is related to the oscillations of the phase rather than the absolute value of the BEC's wave function.
Further, we notice that the spectrum of collective modes might contain an imaginary part. The presence of the imaginary part at small k y for n b = 0 and large k for n b = 0 signifies the dynamical instability of the system. Interestingly, only one of the collective modes (at least at g ab = 0) has a nonzero imaginary part for n b = 0. On the other hand, both modes could have a nonvanishing imaginary part at n b = 0. Furthermore, we note that the dispersion relation of the collective modes becomes trivial, i.e., ω ± = vk, at 2∆ + n a (g a − g ab ) = 0. In this case, n b = 0. gap ∆ can be used to tune the dispersion relation of the collective modes at k ŷ. For example, as one can see from Fig. 8, one of the modes (ω + ) becomes unstable. The Haldane gap can be also used to stabilize collective modes for n b = 0 (see Fig. 9). In particular, ω + and ω − are stable at large values of ∆. Furthermore, we identify three different phases determined by the combination of parameters 2∆ + (g b − g ab )n a at n b = 0. For ∆ < −∆ cr , where ∆ cr = |(g b − g ab )n a |/2, we have one gapped mode and one gapless mode. The latter mode has a nonzero imaginary part and is unstable. These modes merge at ∆ = ∆ cr and have a sound-like dispersion relation ω ± = vk. One of the modes remain gapless while the other acquires a gap for −∆ cr < ∆ < ∆ cr . The modes are unstable for sufficiently large values of momentum when the spectral branches merge. Finally, the degeneracy is lifted and modes become stable at ∆ ≥ ∆ cr . For the parameters used, ω + is gapped and ω − is gapless in this case. Thus, the Haldane gap appears as an efficient tuning parameter that could be used to access different regimes for collective modes.

V. SUMMARY AND OUTLOOK
We investigated the properties of Dirac Bose-Einstein condensate and outline the possible routes to achieve this state in bosonic systems. The proposed framework is general and can be applied for investigating various physical systems. We, however, focus on the case of 2D Dirac magnons. To achieve the Bose-Einstein condensation at the Dirac points, a parametric pumping of magnons was considered. In particular, we introduced a phenomenological model consisting of the Gross-Pitaevskii equations for the Dirac BEC amended with a rate equation for the pumped magnon population. We found that if the pump power reaches a certain threshold determined by the coupling terms between the pumped and Dirac magnon populations, magnons can populate the Dirac nodes. Depending on the values of the coupling terms, a noticeable depletion might be observed for the pumped magnon populations when the BECs form. The profiles of the magnon densities are nonmonotonic. Indeed, they reach a peak value at the onset. Then a plateau is developed for a sufficiently large pumping power. Finally, both condensate and pumped magnon population decay when the pump is turned off. The multicomponent condensate coherence is manifested in the Rabi oscillations of the Dirac magnon populations allowed by the Haldane-like gap term and the Josephson coupling. These features of the magnon populations are readily accessible by the Brillouin light scattering technique used for conventional magnons in YIG.
Dirac nature of magnon BEC is manifested also in the spectrum of collective modes. Even for a simple model at hand, where only the dynamics of a single valley is considered, we found that the spectrum is rather rich. In particular, two uniform (quasi-)ground states are possible, where the magnon density is distributed among the pseudospin (sublattice) degrees of freedom or is accumulated at a certain sublattice. In the first case, the rotation symmetry is spontaneously broken leading to a directional dependence of the dispersion relation. The direction is determined by the phase of the ground state. We found two branches of gapless collective modes in this case that can be identified with the Nambu-Goldstone modes. Depending on the values of parameters, however, one or even both modes could become unstable for small values of wave vector. The case of a maximally anisotropic distribution of magnons among the sublattices is qualitatively different. In general, the spectrum of one of the branches is gapped and can be identified with the Higgs or amplitude mode. For the certain range of parameters, the modes in this case could be also unstable. It is worth noting that while the case of a maximally anisotropic distribution of magnons seems to be unlikely, it could have lower energy. Independent on the ground state, the dynamics of collective modes could be effectively controlled via the Haldane gap term. The latter might be intrinsically present in a system or generated and tuned via the effective pseudo-Zeeman effect.
With the results presented in this paper, we hope that the search for the condensates in bosonic Dirac materials will broaden. We reiterate the Dirac equation as a mathematical structure that can be equally applied to fermions and bosons with specific symmetries [80]. The experimental realization of Dirac magnons and BECs would be a new direction to pursue. Indeed, the majority of the studies in the field of magnon condensates are related to YIG that does not have Dirac points in the magnon spectrum. Trihalides are the promising class of materials that can support the Dirac magnons. A suitable class of materials for investigating Dirac BEC is, however, yet to be identified. We thus expect that the suite of interesting questions about experimental realization, topology, hydrodynamics, edge states, and collective excitations in Dirac BECs would need to be further addressed.
Multicomponent nature of Dirac BEC means that one would have at least four-component hydrodynamics with both normal and superfluid components for the Dirac quasiparticles. Since the condensates are coherent, an interesting problem is the interference of coexisting condensates from different valleys. In a state with different phases of the condensates at different valleys, one would expect spatial interference effects that can be observed by using the Brillouin light scattering, just as for conventional magnon BEC [66]. It would be interesting also to investigate the formation of the condensates in finite samples as well as the corresponding edge (surface) modes.
Finally, let us discuss a few limitations of this study. The employed model contains a minimal number of terms allowed by the symmetries of the system. Clearly, several other terms could be considered. Among them, we mention the terms corresponding to the inter-valley interactions and nonlinear decay terms. While their effects could be interesting, additional terms greatly expand the parameter space of the problem. Therefore, the corresponding studies will be reported elsewhere. Further, the pumping model proposed in this study is phenomenological. The derivation of the corresponding microscopic model is an interesting question that is outside the scope of this study.
are the primitive translation vectors of the hexagonal lattice and a is the spacing between lattice sites. Further, the position vectors for the sublattice B are given by where δ 1 = (a 1 − a 2 )/3 is the relative position of site B with respect to site A in the unit cell. The relative positions of the other two sites of type B are given by δ 2 = δ 1 + a 2 and δ 3 = δ 1 − a 1 .
To bosonise Hamiltonin (A1), we use the standard Holstein-Primakoff transformation truncated to the first order in the inverse spin 1/S. By introducing the bosonic annihilation (creation) operatorsâ i andb i (â † i andb † i ) corresponding to the two sublattices A and B, respectively, we obtain Similar expressions can be written for the sublattice B. By using these relations and performing the Fourier transform where N is the total number of sites and k is momentum, we derive the following Hamiltonian in the leading order approximation (see also Ref. [23]): where and The spectrum of Hamiltonian (A9) reads as As is easy to check, the two branches in Eq. (A12) touch at two non-equivalent isolated points in the Brillouin zone: x, where plus and minus signs correspond to K and K ′ Dirac points or valleys, respectively. Herê r = r/|r| is the unit vector. In the vicinity of K and K ′ points, i.e., at k = k ± D + δk, where δk is small, the expression for ǫ ± k takes the following simple form: where ζ = ± is the valley degree of freedom and v = 3JSa/2 is an analog of the Fermi velocity in graphene. Then, Hamiltonian in Eq. (A9) reads Here σ = (σ x , σ y ) is the vector of the Pauli matrices acting in the sublattice space. In terms of the four-component wave functionΨ we have the following Hamiltonian in sublattice and valley space: This magnon Hamiltonian for has a structure of a massless Dirac Hamiltonian. Therefore, we call such quasiparticles "Dirac magnons". Further, we add the term describing interaction with an effective magnetic field Note that we assume that B i could be different at different sublattices. Effectively, this might stem from slightly different gyromagnetic ratios g at the A and B sublattices. It is worth noting also that the pseudo-Zeeman effect due to strains was predicted for graphene in Ref. [74]. For the sake of definiteness, we assume that B i ẑ. Performing the Holstein-Primakoff transformation and omitting the constant term ∼ gµ B BS, we obtain the following term in the Hamiltonian: As one can see, the usual magnetic field simply shifts 3JS → 3JS + gµ B B. Therefore, it can be used to tune the effective chemical potential of magnons. The effect of the sublattice-dependent fieldB = (B a − B b )/2 is qualitatively different. This term opens the gap in the spectrum and is similar to the Zeeman effect albeit in the pseudospin space. Indeed, the energy spectrum (A12) reads as The structure of the effective fieldB on Dirac magnons allow us to identify it with the Haldane gap term ∆ used in the main text. The energy spectrum of Dirac magnons atB = 0 andB = 0 is shown in the left and right panels of Fig. 1, respectively. The next-order Holstein-Primakoff transformation gives the following interaction term [23]: Let us derive the interaction terms in the vicinity of the Dirac points and neglect the gradient terms. We use   For simplicity we focus on the possible phase lockings between the condensates at two non-equivalent Dirac points. The two amplitudes n a and n b at the two different valleys can be equal to each other or might be different. In the former case, we designate the condensates as symmetric and in the latter case, they are asymmetric. For each of these amplitude classes, we then classify the condensates based on the phase relations. If the phases of condensate components are equal in magnitude at the two valleys as θ a/b,+ = θ a/b,− , we call them co-rotating. In the other case, when the phases are equal in magnitude but differ in sign, the condensates are called as anti co-rotating. However, for suitable interaction channels between the multicomponent condensates, a completely different phase locking can appear as θ a,+ = −θ a,− and θ b,+ = θ b,− (or θ a,+ = θ a,− and θ b,+ = −θ b,− ). In this case, we define it as anti co-rotating "A" ("B") and co-rotating "B" ("A").
Appendix C: Collective modes with Josephson coupling term Let us briefly discuss the effect of the Josephson coupling term u ab on collective modes considered in Sec. IV. The coupling terms complicate the dynamics of the system and does not allow for a simple analytical solution. Therefore, the frequencies of the collective modes are obtained by solving the system of equations (23) to (26) numerically with c 0 = g a = g b = t −1 0 , g ab = 0, and the ground state defined in Sec. IV A. The corresponding results for u ab = 0.2 t −1 0 are shown in Figs. 10 and 11 at n b = 0 and n b = 0, respectively. As one can see, the imaginary part is enhanced compared to the case u ab considered in Sec. IV C. Moreover, since we no longer have independent phases for the ground state wave functions at n b = 0, one of the modes can become gapped (see Fig. 10(a)). Except a single point ∆ = ∆cr, one of the modes is gapped and the other is gapless.