\extrafloats

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

P. O. Sukhachov [email protected] Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden Department of Physics, Yale University, New Haven, CT 06520, USA S. Banerjee [email protected] Theoretical Physics III, Center for Electronic Correlations and Magnetism, Institute of Physics, University of Augsburg, 86135 Augsburg, Germany A. V. Balatsky [email protected] Department of Physics, University of Connecticut, Storrs, CT 06269, USA Nordita, KTH Royal Institute of Technology and Stockholm University, Roslagstullsbacken 23, SE-106 91 Stockholm, Sweden

(January 4, 2021)

Abstract

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 at Dirac nodes, the formation and time evolution of condensates of Dirac bosons 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.

I Introduction

Discrete parity and time-reversal symmetries as well as non-Bravais lattices lead to relativistic-like Dirac energy spectrum in what is called fermionic Dirac and Weyl semimetals Katsnelson:rev-2007 ; Geim-Novoselov:rev-2007 ; CastroNeto-Geim:rev-2009 ; Geim:rev-2009 ; Katsnelson:book-2012 ; Hasan-Kane:rev-2010 ; Qi-Zhang:rev-2010 ; Bernevig:book-2013 ; Franz:book-2013 ; Wehling-Balatsky:rev-2014 ; Yan-Felser:2017-Rev ; Hasan-Huang:rev-2017 ; Armitage-Vishwanath:2017-Rev . The same attributes can be arranged in the case where the lattice is populated with bosons. Thus, one can easily extend the notion of Dirac and Weyl materials to include bosonic Dirac materials that host multicomponent elementary bosonic excitations described effectively by the Dirac or Weyl Hamiltonian. Dirac bosons were realized in photonic Jacqmin-Amo:2014 ; Lu-Soljacic:2015 ; Wang-Jiang:2016 ; Yang-Zhang:2018 ; Guo-Zhang:2018 ; Yang-Chen:2019 and phononic Xiao-Chan:2015 ; Serra-Garcia-Huber:2018 ; Cai-Liu:2020 crystals and magnets Pershoguba-Balatsky:2018 ; Chen-Dai:2018 ; Lee-Hammel:2020 . Honeycomb arrays of superconducting grains were also predicted to host bosonic Dirac excitations Banerjee-Balatsky:2016 . Photonic and phononic nodal excitations are seen in specifically crafted metamaterials that reveal Dirac points when illuminated by electromagnetic or acoustic waves. For example, a linear graphene-like dispersion relation of polaritons was observed in the photoluminescence spectrum of a honeycomb micropillar lattice in Ref. Jacqmin-Amo:2014 .

It is known that a Dirac-type quasiparticle spectrum is generated for localized magnetic moments on a two-dimensional (2D) honeycomb lattice Fransson-Balatsky:2016 ; Pershoguba-Balatsky:2018 . The corresponding magnon Dirac materials Fransson-Balatsky:2016 such as van der Waals crystals Cr₂Ge₂Te₆ Gong-Zhang:2017 and the transition metal trihalides family $AX_{3}$ (where $A=\mbox{Cr}$ and $X=\mbox{F,Cl,Br,I}$ ), e.g., CrI₃ Huang-Xu:2017 ; Pershoguba-Balatsky:2018 ; Chen-Dai:2018 ; Lee-Hammel:2020 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 Gossard-Remeika:1961 ; Davis-Narath:1964 and can be viewed as a magnetic analog of $ABC$ stacked graphite. In three dimensions (3D), a magnonic Weyl state was predicted in pyrochlores Li-Chen:2016 ; Mook-Mertig:2016 ; Owerre:2018 . Linear crossing in the magnon spectrum was indeed observed in the antiferromagnet Cu₃TeO₆ Yao-Li:2018 and the multiferroic ferrimagnet Cu₂OSeO₃ Zhang-Mokrousov:2020 .

A wide range of properties was observed within this emergent class of Dirac magnon materials. We mention interaction effects Pershoguba-Balatsky:2018 ; Banerjee-Balatsky:2020 that renormalize the slope of the dispersion, topological properties Zhang-Li:2013 ; Mook-Mertig:2014 ; Owerre:2016 ; Nakata-Loss:2017 ; Kovalev-Li:2017 ; Owerre:2017 ; Zhang-Mokrousov:2020 ; Wang-Troncoso:2020 , strain effects Ferreiros-Vozmediano:2018 ; Nayga-Vojta:2019 ; Liu-Shi:2020 , and non-Hermitian dynamics McClarty-Rau:2019 . 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.

For the excitations like magnons or polaritons, the equilibrium chemical potential is zero because these excitations are absent in the ground state in equilibrium. Thus, no Bose-Einstein condensate (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 at the Dirac nodes. Such condensates could naturally be viewed as Dirac BECs: they would have multiple component and the nodal spectrum described by an effective Hamiltonian resembling that of Dirac (quasi-)particles.

Dirac nodes are usually located at higher energies where, unlike fermions, a nonvanishing population cannot be achieved via equilibrium doping or gating. Therefore, one would need to produce excitations (i.e., pump the material) to populate the Dirac points, e.g., see the discussion on pumping for magnons below. (For population of the Dirac nodes via light pulses for fermions, see Ref. Pertsova-Balatsky:2020 .)

Bosonic nature of excitations and pumping allow for an accumulation of Dirac bosons in the state with the same energy and ultimately opens up the possibility of Bose-Einstein condensation at the Dirac point. This uncovers a treasure trove of various nontrivial effects discussed for conventional BECs Dalfovo-Stringari:1999-rev ; Ozeri-Davidson:2005-rev ; Pitaevskii-Stringari:book . Recently, properties of BECs with Dirac points in the energy spectrum were studied in Refs. Haddad-Carr:2011a ; Haddad-Carr:2011b ; Haddad-Carr:2015a ; Haddad-Carr:2015b ; Haddad-Carr:2015c ; Haddad-Carr:2015d ; Haddad-Carr:2015e (see also a vast literature on the nonlinear Dirac equation, e.g., Refs. Cuevas-Maraver-Lan:2016 ; Cuevas-Saxena:2018 ; Poddubny-Smirnova:2018 ). Among the most interesting features, we notice various types of vortices and solitons Haddad-Carr:2015a ; Haddad-Carr:2015b ; Haddad-Carr:2015c ; Haddad-Carr:2015d . As plausible systems for the realization of Dirac BECs, cold atoms and honeycomb optical lattices were proposed. In all of these studies, the prior equilibrium BEC was subjected to a perturbing potential to form the Dirac nodes. To the best of our knowledge, however, the possibility and properties of Dirac magnon BECs were not investigated before. Obviously, the ultimate proof of the existence of these condensates would be given by a successful experiment. Here, we lay out the theoretical framework that would facilitate the realizations of the Dirac BEC. The framework should be understood as a model suitable for Bose-Einstein condensates where bosons can be described by a relativistic-like Dirac equation.

Magnons represent a convenient and tested platform for realizing quasiparticle coherent states and BEC in solids. For example, magnon BEC in yttrium–iron–garnet (YIG) films was recently experimentally observed and analyzed Demokritov-Slavin:2006 ; Demidov-Slavin:2007 ; Dzyapko-Slavin:2007 ; Demidov-Slavin:2008 ; Nowik-Boltyk-Demokritov:2012 ; Serga-Hillebrands:2014 ; Bozhko-Hillebrands:2019 ; Mohseni-Pirro:2020 ; Mohseni-Pirro:2020b . (While we focus on magnons, the BEC of polaritons is also possible and was observed in, e.g., Refs. Kasprzak-Dang:2006 ; Deng-Yamamoto:2010 ; Deveaud:2015 ; Sun-Nelson:2017 ; Proukakis-Littlewood:book-2017 .) Among the most notable properties, we mention spatial coherence, supercurrents, and sound modes of the condensate. In order to achieve a magnon BEC, one needs to accumulate magnons in the minima of dispersion relation, which is usually achieved by a certain pumping protocol. Often used is the parametric pumping technique Lvov:book . The central process in this scheme is the splitting of a photon of an external oscillating electromagnetic field into two magnons with arbitrary opposite momenta but with the half-frequency. These magnons subsequently accumulate at the minima of the dispersion and condense if the pumping intensity exceeds certain threshold. In contrast to the case of Dirac magnons, the spin-wave spectrum in YIG is “trivial”, i.e., it does not exhibit any Dirac crossings and is well described by a conventional parabolic dispersion. Moreover, the vast literature investigating the properties of driven BEC in YIG relies on the standard single-component models for magnon BEC, which cannot be applied to Dirac quasiparticles.

For the case under consideration, we assume the nodal structure of spin excitations (Dirac nodes) and aim to describe the formation of long-lived and ultimately condensed magnons at Dirac nodes. We derive a minimal model for a Dirac BEC from the magnetic Hamiltonian, and investigate the time evolution of pumped condensates and the properties of collective modes. While we focus on the case of 2D Dirac magnons, the model is general and can be applied for investigating various Dirac BECs. Since no Fermi surface exists for bosons, we apply pumping to generate a macroscopic occupancy of magnons and ultimately a magnon BEC. This process is phenomenologically described by a multicomponent system, which consists of nonlinear Gross-Pitaevskii equations for a Dirac magnon BEC and a rate equation for a pumped magnon population. As in the case of conventional magnon BECs, the condensate appears after reaching a certain pump power. The effective Haldane gap term, allowed by the Dirac nature of magnons, lifts the degeneracy between the magnon densities for different pseudospins (or sublattices). Furthermore, the gap leads to the Rabi oscillations of the magnon densities when the Josephson coupling terms are present. If experimentally observed, such oscillations would be a definitive signature of the condensate coherence. The nontrivial nature of Dirac BEC is manifested also in the spectrum of collective modes. In our analysis, we used the standard Bogolyubov approach and considered a homogeneous ground state. Interestingly, there are two possibilities for the latter: (i) magnon densities for both pseudospins are nonzero and (ii) one of the densities vanishes. The first state spontaneously breaks the rotational symmetry leading to two gapless Nambu-Goldstone modes. On the other hand, the second ground state is characterized by a gapped Higgs mode and a gapless mode. Depending on the model parameters, collective modes could be unstable. For example, the Haldane gap can be used to access different regimes for collective modes. To the best of our knowledge, the investigation of Dirac bosons, in particular, Dirac magnons, is in its infancy and a Dirac magnon BEC remains yet to be experimentally realized. We hope that our findings will stimulate the corresponding experimental search and the development of the field of Dirac BECs.

The paper is organized as follows. We introduce the model and key notions in Sec. II. While we focus on the case of magnons, the model is general and could be used to describe various Dirac BECs. Pumping and time evolution of Dirac magnons are phenomenologically considered in Sec. III. Section IV is devoted to collective modes in the Dirac BECs. The results are discussed and summarized in Sec. V. A schematic derivation of the free energy ansatz used in the main text is given in Appendix A. A few non-trivial phases of the Dirac BEC are summarized in Appendix B. The additional results for the population dynamics at different interaction constants are presented in Appendix C. Appendix E contains dispersion relations of collective modes for a nonzero Josephson coupling. The model is extended to include the magnon condensate at the $\Gamma$ point and the results for the population dynamics are presented in Appendix D. Through this study, we set $\hbar=k_{B}=1$ .

II Model

In this section, we present a phenomenological approach to Dirac BECs. The approach is similar to that for conventional BECs Pitaevskii-Stringari:book albeit it 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):

	$\displaystyle F$	$\displaystyle=$	$\displaystyle(c_{0}-\mu)n_{\rm tot}+\Delta\sum_{\zeta=\pm}(n_{a,\zeta}-n_{b,\zeta})-iv\Psi^{{\dagger}}(\tau_{z}\otimes\bm{\sigma}\cdot\bm{\nabla})\Psi$		(1)
		$\displaystyle+$	$\displaystyle\sum_{\zeta=\pm}\left[\frac{1}{2}g_{a}n_{a,\zeta}n_{a,\zeta}+\frac{1}{2}g_{b}n_{b,\zeta}n_{b,\zeta}+g_{ab}n_{a,\zeta}n_{b,\zeta}\right]+\frac{1}{2}\sum_{\zeta=\pm}\left(u_{ab}\psi_{a,\zeta}^{}\psi_{a,\zeta}^{}\psi_{b,\zeta}\psi_{b,\zeta}+u_{ab}^{}\psi_{b,\zeta}^{}\psi_{b,\zeta}^{*}\psi_{a,\zeta}\psi_{a,\zeta}\right).$		(1)

The four-component wave function of Dirac BEC is

\Psi=\left(\psi_{a,+},\psi_{b,+},\psi_{b,-},\psi_{a,-}\right)^{T},

(2)

where $\zeta=\pm$ denotes the $K$ or $K^{\prime}$ valley in 2D or chirality in 3D, $\psi_{a,\zeta}$ and $\psi_{b,\zeta}$ are the condensate wave functions for the different pseudospins, e.g., $A$ and $B$ sublattices for 2D hexagonal lattices, $\bm{\sigma}$ is the vector of the Pauli matrices in the pseudospin space, $\tau_{z}$ is the Pauli matrix in the valley space, $n_{i,\zeta}=\psi_{i,\zeta}^{*}\psi_{i,\zeta}$ with $i=a,b$ is the density for a certain pseudospin, $n_{\rm tot}=\sum_{\zeta=\pm}\left(n_{a,\zeta}+n_{b,\zeta}\right)$ is the total condensate density, $c_{0}$ is the position of the Dirac nodes in energy, and $\mu$ is the effective chemical potential of condensates. Finally, the term quantified by $\Delta$ denotes the Haldane gap in 2D. (Note that the term $\propto\Delta$ plays a different role in 3D where it shifts the Weyl nodes of opposite chiralities in momentum space and is called the chiral shift Gorbar:shift-2009 . In what follows, however, we focus only on the 2D case.) Furthermore, we notice 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 resembles 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 are 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 Lin-Spielman:2011 . 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 $\propto n_{a,-\zeta}n_{a,\zeta}$ . 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 and brown 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=\pm k_{\rm D}$ in momentum space and at $\omega_{a}=\omega_{b}=c_{0}$ in energy. As one can see from Fig. 1(b), the presence of the $\Delta$ term opens the gap in the spectrum.

Refer to caption — Figure 1: Schematic representation of the Dirac spectrum and BECs. The condensates of Dirac magnons occur at the Dirac points located at $k=\pm k_{\rm D}$ in momentum space. Panels (a) and (b) correspond to gapless $\Delta=0$ and gapped $\Delta\neq 0$ spectra, respectively. Each of the condensates contains two components related to the pseudospin (sublattice) degree of freedom. The corresponding densities are denoted as $n_{a,\zeta}$ and $n_{b,\zeta}$ with $\zeta=\pm$ . The density $n_{3}$ corresponds to the pumped population of bosons with energy $\omega_{3}$ (see Sec. III for the description of the pumping mechanism). Finally, $n_{0}$ denotes the population at the $\Gamma$ point (see also Appendix D).

By varying the free energy (1) with respect to the fields $\psi_{a,\zeta}^{*}$ and $\psi_{b,\zeta}^{*}$ , we derive the following Gross-Pitaevskii equations for Dirac BECs at $A$ and $B$ sublattices:

	$\displaystyle i\partial_{t}\psi_{a,\zeta}$	$\displaystyle=$	$\displaystyle\left(c_{0}-\mu+\Delta\right)\psi_{a,\zeta}-iv\left[\left(\zeta\partial_{x}-i\partial_{y}\right)\psi_{b,\zeta}\right]+g_{a}n_{a,\zeta}\psi_{a,\zeta}+g_{ab}n_{b,\zeta}\psi_{a,\zeta}+u_{ab}\psi_{a,\zeta}^{*}\left(\psi_{b,\zeta}\right)^{2},$		(3)
	$\displaystyle i\partial_{t}\psi_{b,\zeta}$	$\displaystyle=$	$\displaystyle\left(c_{0}-\mu-\Delta\right)\psi_{b,\zeta}-iv\left[\left(\zeta\partial_{x}+i\partial_{y}\right)\psi_{a,\zeta}\right]+g_{b}n_{b,\zeta}\psi_{b,\zeta}+g_{ab}n_{a,\zeta}\psi_{b,\zeta}+u_{ab}^{}\psi_{b,\zeta}^{}\left(\psi_{a,\zeta}\right)^{2}.$		(4)

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+}\times U(1)_{a-}\times U(1)_{b+}\times 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 and there are only two flavor symmetries, which are explicitly broken to $U(1)_{a}\times U(1)_{b}$ . In other words, we consider a symmetric response of the valleys. Therefore, we fix $\zeta=+$ 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 and a few other details depends on the nature of Bose-Einstein condensates and, in general, are 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. Pertsova-Balatsky:2020 .) As we discussed in the introduction, Dirac magnons can be realized in, e.g., CrI₃ (see Refs. Pershoguba-Balatsky:2018 ; Chen-Dai:2018 ; Lee-Hammel:2020 ). 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 for fermions in graphene, see Ref. Manes-Vozmediano:2013 ). In addition, the relatively large gap in the Dirac spectrum could appear due to the Kitaev interaction. For a recent study in CrI₃, see, e.g., Ref. Lee-Hammel:2020 .

III.1 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 $\omega_{3}$ (see also Fig. 1). The extension of this model to the case where the magnon BEC is formed at the $\Gamma$ point is presented in Appendix D. These populations could be achieved by the parametric pumping similarly to the case of conventional magnetic materials Lvov:book or any other suitable pumping scheme. The results below do not depend on the details of the pumping. The system reads as

	$\displaystyle i\partial_{t}\psi_{a}+i\Gamma_{a}\psi_{a}=\left(c_{0}-\mu+\Delta\right)\psi_{a}-iv\left[\left(\partial_{x}-i\partial_{y}\right)\psi_{b}\right]+g_{a}n_{a}\psi_{a}+g_{ab}n_{b}\psi_{a}+u_{ab}\psi_{a}^{*}\left(\psi_{b}\right)^{2}+iP_{a3}\Gamma_{a}n_{3}\psi_{a},$		(5)
	$\displaystyle i\partial_{t}\psi_{b}+i\Gamma_{b}\psi_{b}=\left(c_{0}-\mu-\Delta\right)\psi_{b}-iv\left[\left(\partial_{x}+i\partial_{y}\right)\psi_{a}\right]+g_{b}n_{b}\psi_{b}+g_{ab}n_{a}\psi_{b}+u_{ab}^{}\psi_{b}^{}\left(\psi_{a}\right)^{2}+iP_{b3}\Gamma_{b}n_{3}\psi_{b},$		(6)
	$\displaystyle\partial_{t}n_{3}+\Gamma_{3}n_{3}=\Gamma_{3}\tilde{P}(t)-\sum_{\zeta=\pm}\left(P_{3a}\Gamma_{a}n_{a,\zeta}+P_{3b}\Gamma_{b}n_{b,\zeta}\right)n_{3}.$		(7)

Compared to Eqs. (3) and (4), we have added the dissipation terms $i\Gamma_{a}\psi_{a}$ and $i\Gamma_{b}\psi_{b}$ that correspond to the magnon-phonon 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, $\Gamma_{3}$ denotes the decay rate. Finally, the term $\Gamma_{3}\tilde{P}(t)$ describes the direct inflow of magnons caused by the pump. (Note that explicit dynamics of noncondensed magnons Hahn-Kopietz:2020 at the Dirac nodes is not included in the above system. In the case under consideration, it is implicitly taken into account in the inflow terms.) Therefore, the formation of the condensates is a multistep process where the pumping first creates a magnon gas which decays into the Dirac magnons.

For the sake of simplicity, we consider a simple step-like profile for the pump power

\tilde{P}(t)=\sqrt{P_{\rm max}}\,\theta(t-t_{\rm Begin})\theta(t_{\rm End}-t).

(8)

Here, the pump starts at $t=t_{\rm Begin}$ and ends at $t=t_{\rm End}$ . Further, $P_{\rm max}$ is the power of the pumping field and $\theta(x)$ is the step function. The square-root dependence on the pump power is in agreement with the experimental data for YIG Demidov-Slavin:2008 ; Demokritov-Slavin:2008 (see also Refs. Rezende:2009 ; Rezende:2010 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 time-dependent pump.

III.2 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:

	$\displaystyle\Gamma_{a}=\Gamma_{b}=\Gamma_{3}=0.25\,t_{0}^{-1},\quad c_{0}=t_{0}^{-1},\quad\mu=0,\quad\Delta=0,\quad g_{a}=g_{b}=t_{0}^{-1},\quad g_{ab}=u_{ab}=0,$		(9)
	$\displaystyle P_{3a}=P_{3b}=P_{a3}=P_{b3}=1,\quad t_{\rm Begin}=10\,t_{0},\quad t_{\rm End}=200\,t_{0}.$		(10)

Here, $t_{0}$ is the characteristic interaction timescale. Since no experimental data are present for the Dirac magnon BEC, we use model parameters. In addition, we note that, in general, $P_{i3}\neq 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. Bozhko-Hillebrands:2019 . Further, because of the space separation between the lattice sites, 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}|\ll g_{a},g_{b}$ . The results for nonvanishing $g_{ab}$ are presented in Appendix C. Qualitatively, however, they are similar to the case $|g_{ab}|\ll g_{a},g_{b}$ . Finally, we note that the interaction constants $g_{a}$ , $g_{b}$ , and $g_{ab}$ are irrelevant for the condensate dynamics in the absence of the phase-mixing terms [e.g., the terms with $u_{ab}$ in Eqs. (5) and (6)]. This follows from the fact that the density dynamics is determined by the imaginary part of the Gross-Pitaevskii equations (5) and (6).

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 $\tau\sim 1/\Gamma$ . 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 $\Delta$ . The corresponding results are shown in Fig. 2 at a few values of the gap for $P_{\rm 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 $\propto u_{ab}$ leads to the oscillations of the Dirac BEC. Oscillations of $n_{a}$ and $n_{b}$ always occur in the 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 $\Delta$ in Fig. 5. The period indeed scales as $1/\Delta$ for large values of $\Delta$ . The case of small $\Delta$ 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 $\Delta\lesssim t_{0}^{-1}$ .

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. In the systems with small distance between the Dirac nodes in the momentum space, such observations could be, in principle, done via the Brillouin light scattering technique used for conventional magnons in YIG. In the case of relatively large wave vectors, different techniques might be required to resolve magnons belonging to different Dirac nodes or sublattices. Among them, we mention time-resolved magneto-optic Kerr effect measurements and nitrogen-vacancy (NV) centers (see Refs. Prananto-An:2020 ; Kruglyak-Grundler:2010 for reviews).

The dynamics of the Dirac magnon populations can be presented by introducing the vector in the density space $\left(n_{a},n_{b},n_{3}\right)/n_{\rm 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 SM . 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 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}\neq 0$ and $n_{b}\neq 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 gapped and corresponds to the Higgs or amplitude mode, where the populations at the $A$ and $B$ sublattices oscillate.

In order to investigate the spectrum of collective excitations, we employ the Bogolyubov approach Dalfovo-Stringari:1999-rev ; Ozeri-Davidson:2005-rev ; Pitaevskii-Stringari:book for the magnon BEC. In this case, one looks for the solution as

	$\displaystyle\psi_{a}$	$\displaystyle=$	$\displaystyle\psi_{a,0}(\mathbf{r})+e^{-i\mu t}\left[u_{a}e^{-i\omega t+i\mathbf{k}\cdot\mathbf{r}}+v_{a}^{}e^{i\omega^{}t-i\mathbf{k}\cdot\mathbf{r}}\right],$		(11)
	$\displaystyle\psi_{b}$	$\displaystyle=$	$\displaystyle\psi_{b,0}(\mathbf{r})+e^{-i\mu t}\left[u_{b}e^{-i\omega t+i\mathbf{k}\cdot\mathbf{r}}+v_{b}^{}e^{i\omega^{}t-i\mathbf{k}\cdot\mathbf{r}}\right].$		(12)

Here, $\psi_{a,0}(\mathbf{r})$ and $\psi_{b,0}(\mathbf{r})$ denote ground-state solutions, which are, in general, nonuniform. They satisfy the time-independent systems of the Gross-Pitaevskii equations. The perturbations are described by $u_{a,b}$ and $v_{a,b}^{*}$ terms. Finally, $\omega$ is the frequency and $\mathbf{k}$ is the wave vector of collective modes.

IV.1 Ground state

Let us start with the ground-state solutions $\psi_{a,0}(\mathbf{r})$ and $\psi_{b,0}(\mathbf{r})$ . For simplicity, we consider a uniform ground state. By using Eqs. (11) and (12), the Gross-Pitaevskii equations (3) and (4) lead to

	$\displaystyle 0=\left(c_{0}-\mu+\Delta\right)\psi_{a,0}+g_{a}n_{a}\psi_{a,0}+g_{ab}n_{b}\psi_{a,0}+u_{ab}\psi_{a,0}^{*}(\psi_{b,0})^{2},$		(13)
	$\displaystyle 0=\left(c_{0}-\mu-\Delta\right)\psi_{b,0}+g_{b}n_{b}\psi_{b,0}+g_{ab}n_{a}\psi_{b,0}+u_{ab}^{}\psi_{b,0}^{}(\psi_{a,0})^{2}.$		(14)

It is convenient to separate the absolute value and the phase of the ground-state wave functions, i.e., $\psi_{a,0}=\sqrt{n_{a}}e^{i\theta_{a}}$ and $\psi_{b,0}=\sqrt{n_{b}}e^{i\theta_{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}\times U(1)_{b}$ symmetry. If the Josephson coupling term is nonzero, we have

\theta_{u_{ab}}=2(\theta_{a}-\theta_{b})+\pi l,\quad l=0,1,2,3,\ldots,

(15)

where $u_{ab}=|u_{ab}|e^{i\theta_{u_{ab}}}$ . Further, we find an effective chemical potential that allows for a solution of the Gross-Pitaevskii equations. We obtain

	$\displaystyle\mu$	$\displaystyle=$	$\displaystyle c_{0}+\Delta+g_{a}n_{a}+g_{ab}n_{b}+u_{ab}n_{b},$		(16)
	$\displaystyle\mu$	$\displaystyle=$	$\displaystyle c_{0}-\Delta+g_{b}n_{b}+g_{ab}n_{a}+u_{ab}n_{a}.$		(17)

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

n_{b}=\frac{2\Delta+\left(g_{a}-g_{ab}-u_{ab}\right)n_{a}}{g_{b}-g_{ab}-u_{ab}}.

(18)

(One can fix $n_{b}$ and solve for $n_{a}$ , which is physically equivalent). The effective chemical potential for the density (18) reads as

\displaystyle\mu=c_{0}+\Delta+g_{a}n_{a}+\left(g_{ab}+u_{ab}\right)\frac{2\Delta+\left(g_{a}-g_{ab}-u_{ab}\right)n_{a}}{g_{b}-g_{ab}-u_{ab}}.

(19)

There is also an additional solution for Eqs. (13) and (14), where $n_{b}=0$ and

\mu=c_{0}+\Delta+g_{a}n_{a}.

(20)

We note that solution (18) exists as long as $n_{b}>0$ . Otherwise, only the ground state with $n_{b}=0$ and $\mu$ 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}\neq 0$ and $n_{b}\neq 0$ and the solution that breaks the pseudospin symmetry with $n_{a}\neq 0$ but $n_{b}=0$ . The latter solution is an example of nematic instability similar to nematic states in fermion superconductors Fradkin-Mackenzie:2010 and in Bose nematics (see, e.g., Ref. Jian-Zhai:2011 ).

IV.2 Free energy

In order to clarify which of the two possible solutions for the steady-state homogenous Gross-Pitaevskii equations, discussed in Sec. IV.1, is the true ground state, we calculate the free energy of the system. By using Eq. (1), we derive the following expression:

	$\displaystyle F_{n_{b}\neq 0}$	$\displaystyle=$	$\displaystyle\frac{1}{\left(g_{ab}+u_{ab}-g_{b}\right)^{2}}\Bigg{\{}g_{b}\left\{n_{a}^{2}\left[(g_{ab}+u_{ab})^{2}+2g_{a}(g_{ab}+u_{ab})-g_{a}^{2}\right]-4g_{a}n_{a}\Delta-4\Delta^{2}\right\}$		(21)
		$\displaystyle-$	$\displaystyle n_{a}(g_{ab}+u_{ab})^{2}\left[n_{a}\left(2(g_{ab}+u_{ab})-g_{a}\right)-4\Delta\right]-g_{a}g_{b}^{2}n_{a}^{2}\Bigg{\}},$		(21)

where $n_{b}$ given in Eq. (18) was used. In the case $n_{b}=0$ , we have a simpler expression

\displaystyle F_{n_{b}=0}=-3g_{a}n_{a}^{2}.

(22)

In addition to comparing the free energies, one should check whether $n_{b}$ remains positive for the state with $n_{b}\neq 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) and (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}\neq 0$ and $n_{b}=0$ are possible. In general, however, the ground state with $n_{b}\neq 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)].

IV.3 Dispersion relations

In this section, we use the solutions found in Sec. IV.1 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}$ :

$\displaystyle\omega u_{a}$	$\displaystyle=$	$\displaystyle\left(c_{0}+\Delta-\mu\right)u_{a}+v(k_{x}-ik_{y})u_{b}+2g_{a}n_{a}u_{a}+g_{ab}n_{b}u_{a}+g_{ab}\psi_{a,0}\psi_{b,0}^{*}u_{b}+g_{a}\psi_{a,0}^{2}v_{a}+g_{ab}\psi_{a,0}\psi_{b,0}v_{b}$	(23)
	$\displaystyle+$	$\displaystyle u_{ab}\left[2\psi_{b,0}\psi_{a,0}^{*}u_{b}+(\psi_{b,0})^{2}v_{a}\right],$	(23)
$\displaystyle\omega u_{b}$	$\displaystyle=$	$\displaystyle\left(c_{0}-\Delta-\mu\right)u_{b}+v(k_{x}+ik_{y})u_{a}+2g_{b}n_{b}u_{b}+g_{ab}n_{a}u_{b}+g_{ab}\psi_{b,0}\psi_{a,0}^{*}u_{a}+g_{b}\psi_{b,0}^{2}v_{b}+g_{ab}\psi_{a,0}\psi_{b,0}v_{a}$	(24)
	$\displaystyle+$	$\displaystyle u_{ab}^{}\left[2\psi_{a,0}\psi_{b,0}^{}u_{a}+(\psi_{a,0})^{2}v_{b}\right],$	(24)
$\displaystyle-\omega v_{a}$	$\displaystyle=$	$\displaystyle\left(c_{0}+\Delta-\mu\right)v_{a}-v(k_{x}+ik_{y})v_{b}+2g_{a}n_{a}v_{a}+g_{ab}n_{b}v_{a}+g_{ab}\psi_{a,0}^{}\psi_{b,0}v_{b}+g_{a}(\psi_{a,0}^{})^{2}u_{a}+g_{ab}\psi_{a,0}^{}\psi_{b,0}^{}u_{b}$	(25)
	$\displaystyle+$	$\displaystyle u_{ab}^{}\left[2\psi_{b,0}^{}\psi_{a,0}v_{b}+(\psi_{b,0}^{*})^{2}u_{a}\right],$	(25)
$\displaystyle-\omega v_{b}$	$\displaystyle=$	$\displaystyle\left(c_{0}-\Delta-\mu\right)v_{b}-v(k_{x}-ik_{y})v_{a}+2g_{b}n_{b}v_{b}+g_{ab}n_{a}v_{b}+g_{ab}\psi_{a,0}\psi_{b,0}^{}v_{a}+g_{b}(\psi_{b,0}^{})^{2}u_{b}+g_{ab}\psi_{a,0}^{}\psi_{b,0}^{}u_{a}$	(26)
	$\displaystyle+$	$\displaystyle u_{ab}\left[2\psi_{a,0}^{}\psi_{b,0}v_{a}+(\psi_{a,0}^{})^{2}u_{b}\right].$	(26)

As before, we omit the valley index $\zeta$ 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 phase-locked condensates ( $\theta_{a}=\theta_{b}=\theta_{u_{ab}}=0$ ). For a more detailed discussion and a few dispersion relations at $u_{ab}\neq 0$ , see Appendix E.

In the case of the nontrivial ground state solution given in Eq. (18), the frequencies are

\omega_{\pm}=\sqrt{v^{2}k^{2}\pm\frac{2v|k_{y}|}{g_{ab}-g_{b}}\sqrt{\left(g_{ab}-g_{b}\right)\left(g_{ab}^{2}-g_{a}g_{b}\right)\left[2\Delta+\left(g_{a}-g_{ab}\right)n_{a}\right]}}.

(27)

Note that both branches are gapless. In addition, we fix the phase of the ground state solutions as $\theta_{a}=\theta_{b}=0$ . Even at $\Delta=0$ , the dispersion relation (27) is rather interesting. For example, the mode is stable (i.e., there is no large imaginary part) for $\left(g_{ab}^{2}-g_{a}g_{b}\right)\left(g_{ab}-g_{a}\right)\left(g_{ab}-g_{b}\right)<0$ . In the case $g_{a}=g_{b}$ , the stability criterion is simple $g_{a}^{2}\geq g_{ab}^{2}$ , which agrees with phase-separation criterion for conventional BECs Pitaevskii-Stringari:book .

For the ground state with $n_{b}=0$ , one of the collective modes (e.g., $\omega_{+}$ ) can be gapped. The full spectrum in this case is

	$\displaystyle\omega^{2}_{\pm}$	$\displaystyle=$	$\displaystyle v^{2}k^{2}+\frac{\left[2\Delta+n_{a}(g_{a}-g_{ab})\right]^{2}}{2}$		(28)
		$\displaystyle\pm$	$\displaystyle\frac{1}{2}\sqrt{\left\{2v^{2}k^{2}+\left[2\Delta+n_{a}(g_{a}-g_{ab})\right]^{2}\right\}^{2}-4\left\{v^{4}k^{4}+v^{2}k^{2}g_{a}n_{a}\left[2\Delta+(g_{a}-g_{ab})n_{a}\right]\right\}}.$		(28)

The gap between $\omega_{+}$ and $\omega_{-}$ branches reads as

\left|\omega_{+}(k=0)-\omega_{-}(k=0)\right|=\left|2\Delta+(g_{b}-g_{ab})n_{a}\right|.

(29)

Note that the gap is determined both by the Haldane gap term and interaction terms. In a simple case $\Delta=0$ , the modes at $n_{b}=0$ are stable at $g_{a}^{2}\geq g_{ab}^{2}$ . In addition, while the presence of the gapless modes for 2D materials might naively signify an instability towards long-range power-law correlations, it should not play a crucial role in experimental samples. Indeed, in such a case, the range of correlations is effectively limited by the size of samples and the range of probes.

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}\neq 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 ground-state solutions $\psi_{a,0}$ and $\psi_{b,0}$ . For example, one could replace $k_{y}\to k_{x}$ for $\theta_{b}-\theta_{a}=\pi/2$ . It is worth noting, however, that the state $n_{b}\neq 0$ is not always the ground state of the system (see the discussion in Sec. IV.2). Depending on the values of the interaction constants, the state with $n_{b}=0$ could be also realized. 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 $\psi_{a,0}\to e^{i\theta_{a}}\psi_{a,0}$ and $\psi_{b,0}\to e^{i\theta_{b}}\psi_{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 $\psi_{a,0}\to e^{i\theta_{a}}\psi_{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}\neq 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}\neq 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., $\omega_{\pm}=vk$ , at $2\Delta+n_{a}(g_{a}-g_{ab})=0$ . In this case, $n_{b}=0$ .

Let us discuss the dependence of the frequencies of the collective modes on the Haldane gap parameter $\Delta$ . As expected from the analytical result (27), there is no dependence on $\Delta$ for $\mathbf{k}\parallel\hat{\mathbf{x}}$ . On the other hand, the Haldane gap $\Delta$ can be used to tune the dispersion relation of the collective modes at $\mathbf{k}\parallel\hat{\mathbf{y}}$ . For example, as one can see from Fig. 8, one of the modes ( $\omega_{+}$ ) becomes unstable. The Haldane gap can be also used to stabilize collective modes for $n_{b}=0$ (see Fig. 9). In particular, $\omega_{+}$ and $\omega_{-}$ are stable at large values of $\Delta$ .

Furthermore, we identify three different phases determined by the combination of parameters $2\Delta+(g_{b}-g_{ab})n_{a}$ at $n_{b}=0$ . For $\Delta<-\Delta_{\rm cr}$ , where $\Delta_{\rm 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 $\Delta=\Delta_{\rm cr}$ and have a sound-like dispersion relation $\omega_{\pm}=vk$ . One of the modes remains gapless while the other acquires a gap for $-\Delta_{\rm cr}<\Delta<\Delta_{\rm 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 $\Delta\geq\Delta_{\rm cr}$ . For the parameters used, $\omega_{+}$ is gapped and $\omega_{-}$ 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 the Bose-Einstein condensate of Dirac quasiparticles and outlined possible routes to achieve this state in bosonic systems. The proposed framework is general and can be applied for investigating various physical systems with Dirac bosonic dispersion. We focus on the case of 2D Dirac magnons. To achieve the Bose-Einstein condensation at the Dirac points, a pumping of magnons was proposed. 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-type 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 Nowik-Boltyk-Demokritov:2012 . Among other techniques that might be required to resolve magnons belonging to different Dirac nodes or sublattices, we mention time-resolved magneto-optic Kerr effect measurements Kruglyak-Grundler:2010 and nitrogen-vacancy centers Prananto-An:2020 .

The 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 of the ground state, the dynamics of collective modes could be effectively controlled via the Haldane gap term. This term might be intrinsically present in a system or generated and tuned via the effective pseudo-Zeeman effect in ferrimagnets.

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 Larsson-Balatsky:2019 . As an explicit example, we considered the case of Dirac magnons, whose experimental realization and condensation 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 Dirac magnons. A suitable class of materials for investigating the 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.

The 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 could be observed by using the Brillouin light scattering, just as for conventional magnon BEC Nowik-Boltyk-Demokritov:2012 , as well as time-resolved magneto-optic Kerr effect measurements Kruglyak-Grundler:2010 and nitrogen-vacancy centers Prananto-An:2020 . It would be interesting also to investigate the formation of the condensates in finite samples and the effect of 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 as well as its application to realistic materials is an interesting question that is outside the scope of this study.

Acknowledgements.

We appreciate useful discussions with G. Aeppli, J. Bailey, A. Pertsova, A. Sinner, C. Triola, A. Yakimenko, and V. Zyuzin. This work was supported by the University of Connecticut, VILLUM FONDEN via the Centre of Excellence for Dirac Materials (Grant No. 11744), the European Research Council under the European Unions Seventh Framework Program Synergy HERO, and the Knut and Alice Wallenberg Foundation KAW 2018.0104. PS acknowledges the support through the Yale Prize Postdoctoral Fellowship in Condensed Matter Theory. SB acknowledges support from Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)-TRR 80.

Appendix A Schematic derivation of the energy density for Dirac magnons

In this appendix, we provide a schematic derivation of the free energy density (1) by using magnons on a honeycomb lattice as a characteristic example. The Heisenberg Hamiltonian of the model reads as

\displaystyle H=-J\sum_{\langle ij\rangle}\mathbf{S}_{i}\cdot\mathbf{S}_{j},

(30)

where the summation runs over the nearest-neighbor coupling between spins $\mathbf{S}_{i}$ . Further, $J$ is the coupling constant, which is positive in the ferromagnetic phase.

The position vectors of lattice sites for the sublattice $A$ could be parameterized as follows:

\mathbf{r}_{\mathbf{n}}=n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2},

(31)

where the site label $\mathbf{n}=(n_{1},n_{2})$ is determined by a pair of integers $n_{1}$ and $n_{2}$ , and

\mathbf{a}_{1}=\left(\frac{\sqrt{3}a}{2},\frac{3a}{2}\right),\qquad\mathbf{a}_{2}=\left(\frac{\sqrt{3}a}{2},-\frac{3a}{2}\right)

(32)

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

\mathbf{r}_{\mathbf{n}}^{\prime}=\bm{\delta}_{1}+n_{1}\mathbf{a}_{1}+n_{2}\mathbf{a}_{2},

(33)

where $\bm{\delta}_{1}=(\mathbf{a}_{1}-\mathbf{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 $\bm{\delta}_{2}=\bm{\delta}_{1}+\mathbf{a}_{2}$ and $\bm{\delta}_{3}=\bm{\delta}_{1}-\mathbf{a}_{1}$ .

To bosonise Hamiltonian (30), 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 $\hat{a}_{i}$ and $\hat{b}_{i}$ ( $\hat{a}_{i}^{{\dagger}}$ and $\hat{b}_{i}^{{\dagger}}$ ) corresponding to the two sublattices $A$ and $B$ , respectively, we obtain

$\displaystyle S_{i}^{+}$	$\displaystyle\equiv$	$\displaystyle S_{i}^{x}+iS_{i}^{y}=\sqrt{2S}\left(\hat{a}_{i}-\frac{\hat{a}_{i}^{{\dagger}}\hat{a}_{i}\hat{a}_{i}}{4S}\right)+\mathcal{O}\left(S^{-3/2}\right),$	(34)
$\displaystyle S_{i}^{-}$	$\displaystyle\equiv$	$\displaystyle S_{i}^{x}-iS_{i}^{y}=\sqrt{2S}\left(\hat{a}_{i}^{{\dagger}}-\frac{\hat{a}_{i}^{{\dagger}}\hat{a}_{i}^{{\dagger}}\hat{a}_{i}}{4S}\right)+\mathcal{O}\left(S^{-3/2}\right),$	(35)
$\displaystyle S_{i}^{z}$	$\displaystyle=$	$\displaystyle S-\hat{a}_{i}^{{\dagger}}\hat{a}_{i}.$	(36)

Similar expressions can be written for the sublattice $B$ . By using these relations and performing the Fourier transform $\hat{a}_{i}=N^{-1/2}\sum_{\mathbf{k}}e^{i\mathbf{k}\bm{\delta}_{i}}\hat{a}_{\mathbf{k}}$ , where $N$ is the total number of sites and $\mathbf{k}$ is momentum, we derive the following Hamiltonian in the leading order approximation (see also Ref. Pershoguba-Balatsky:2018 ):

\displaystyle H=\sum_{\mathbf{k}}\hat{\psi}^{{\dagger}}_{\mathbf{k}}H_{0}\hat{\psi}_{\mathbf{k}},

(37)

where

\displaystyle H_{0}=JS\left(\begin{array}[]{cc}3&-\gamma_{\mathbf{k}}\\ -\gamma_{\mathbf{k}}^{*}&3\\ \end{array}\right),

(40)

\displaystyle\hat{\psi}_{\mathbf{k}}=\left(\begin{array}[]{c}\hat{a}_{\mathbf{k}}\\ \hat{b}_{\mathbf{k}}\\ \end{array}\right),

(43)

and

\gamma_{\mathbf{k}}=\sum_{j}e^{i\mathbf{k}\cdot\bm{\delta}_{j}}=\left[e^{iak_{y}}+2e^{-\frac{i}{2}ak_{y}}\cos\left(\frac{\sqrt{3}ak_{x}}{2}\right)\right].

(44)

The spectrum of Hamiltonian (40) reads as

\displaystyle\epsilon_{\mathbf{k}}^{\pm}=JS\left(3\pm|\gamma_{\mathbf{k}}|\right)=3JS\pm JS\sqrt{1+4\cos\left(\frac{\sqrt{3}ak_{x}}{2}\right)\left[\cos\left(\frac{\sqrt{3}ak_{x}}{2}\right)+\cos\left(\frac{3ak_{y}}{2}\right)\right]}.

(45)

As is easy to check, the two branches in Eq. (45) touch at two non-equivalent isolated points in the Brillouin zone: $\mathbf{k}_{D}^{\pm}=\pm 4\pi/(3\sqrt{3}a)\hat{\mathbf{x}}$ , where plus and minus signs correspond to $K$ and $K^{\prime}$ Dirac points or valleys, respectively. Here, $\hat{\mathbf{r}}=\mathbf{r}/|\mathbf{r}|$ is the unit vector. In the vicinity of $K$ and $K^{\prime}$ points, i.e., at $\mathbf{k}=\mathbf{k}_{D}^{\pm}+\delta\mathbf{k}$ , where $\delta\mathbf{k}$ is small, the expression for $\epsilon_{\mathbf{k}}^{\pm}$ takes the following simple form:

\epsilon_{\mathbf{k}_{D}^{\zeta}+\delta\mathbf{k}}^{\pm}\simeq 3JS\pm v|\delta\mathbf{k}|,

(46)

where $\zeta=\pm$ is the valley degree of freedom and $v=3JSa/2$ is an analog of the Fermi velocity in graphene. Then, the Hamiltonian in Eq. (40) reads as

\displaystyle H_{0}^{(\zeta)}=3JS\mathds{1}_{2}+v\zeta\sigma_{x}\delta k_{x}+v\sigma_{y}\delta k_{y}.

(47)

Here, $\bm{\sigma}=\left(\sigma_{x},\sigma_{y}\right)$ is the vector of the Pauli matrices acting in the sublattice space. In terms of the bispinor

\hat{\Psi}_{\delta\mathbf{k}}=\left(\hat{a}_{+,\delta\mathbf{k}},\hat{b}_{+,\delta\mathbf{k}},\hat{b}_{-,\delta\mathbf{k}},\hat{a}_{-,\delta\mathbf{k}}\right)^{T},

(48)

we have the following Hamiltonian in sublattice and valley space:

H_{0}\simeq\hat{\Psi}_{\delta\mathbf{k}}^{\dagger}\left(\begin{array}[]{cc}3JS+v(\delta\mathbf{k}\cdot\bm{\sigma})&0\\ 0&3JS-v(\delta\mathbf{k}\cdot\bm{\sigma})\end{array}\right)\hat{\Psi}_{\delta\mathbf{k}}.

(49)

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 $\mathbf{B}_{i}$

-g\mu_{\rm B}\sum_{i}\left(\mathbf{B}_{i}\cdot\mathbf{S}_{i}\right).

(50)

Note that we assume that $\mathbf{B}_{i}$ could be different at different sublattices. Effectively, this might stem from slightly different gyromagnetic ratios $g$ at the $A$ and $B$ sublattices, which indeed occurs in ferrimagnets Fransson-Balatsky:2016 . It is worth noting also that the pseudo-Zeeman effect was predicted for fermions in graphene in Ref. Manes-Vozmediano:2013 . Its origin, however, is related to strains. For the sake of definiteness, we assume that $\mathbf{B}_{i}\parallel\hat{\mathbf{z}}$ . Performing the Holstein-Primakoff transformation and omitting the constant term $\sim g\mu_{\rm B}BS$ , we obtain the following term in the Hamiltonian:

H_{\rm B}=g\mu_{\rm B}\sum_{\mathbf{{k}}}\hat{\psi}_{\mathbf{{k}}}^{{\dagger}}\left(B\mathds{1}_{2}+\tilde{B}\sigma_{z}\right)\hat{\psi}_{\mathbf{{k}}}.

(51)

As one can see, the usual magnetic field simply shifts the spectrum of magnons as $3JS\to 3JS+g\mu_{\rm B}B$ . Therefore, it can be used to tune the effective chemical potential of magnons. Furthermore, it stabilizes the 2D Dirac ferromagnetic material. The effect of the sublattice-dependent field $\tilde{B}=(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 (45) reads as

\epsilon_{\mathbf{k}}^{\pm}=JS\left(3\pm\sqrt{|\gamma_{\mathbf{k}}|^{2}+\left(g\mu_{\rm B}\tilde{B}\right)^{2}/(JS)^{2}}\right).

(52)

The structure of the effective field $\tilde{B}$ for Dirac magnons allows us to identify it with the Haldane gap term $\Delta$ used in the main text. The energy spectrum of Dirac magnons at $\tilde{B}=0$ and $\tilde{B}\neq 0$ is shown in Figs. 1(a) and 1(b), respectively.

The next-order Holstein-Primakoff transformation gives the following interaction term Pershoguba-Balatsky:2018 :

	$\displaystyle H_{4}$	$\displaystyle=$	$\displaystyle\frac{J}{4N}\sum_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4}}\delta_{\mathbf{k}_{1}+\mathbf{k}_{2},\mathbf{k}_{3}+\mathbf{k}_{4}}\Bigg{\{}\gamma_{\mathbf{k}_{2}}^{}\hat{a}_{\mathbf{k}_{1}}^{{\dagger}}\hat{b}_{\mathbf{k}_{2}}^{{\dagger}}\hat{a}_{\mathbf{k}_{3}}\hat{a}_{\mathbf{k}_{4}}+\gamma_{\mathbf{k}_{4}}\hat{a}_{\mathbf{k}_{1}}^{{\dagger}}\hat{a}_{\mathbf{k}_{2}}^{{\dagger}}\hat{a}_{\mathbf{k}_{3}}\hat{b}_{\mathbf{k}_{4}}+\gamma_{\mathbf{k}_{2}}\hat{b}_{\mathbf{k}_{1}}^{{\dagger}}\hat{a}_{\mathbf{k}_{2}}^{{\dagger}}\hat{b}_{\mathbf{k}_{3}}\hat{b}_{\mathbf{k}_{4}}+\gamma_{\mathbf{k}_{4}}^{}\hat{b}_{\mathbf{k}_{1}}^{{\dagger}}\hat{b}_{\mathbf{k}_{2}}^{{\dagger}}\hat{b}_{\mathbf{k}_{3}}\hat{a}_{\mathbf{k}_{4}}$		(53)
		$\displaystyle-$	$\displaystyle 4\gamma_{\mathbf{k}_{4}-\mathbf{k}_{2}}^{*}\hat{a}_{\mathbf{k}_{1}}^{{\dagger}}\hat{b}_{\mathbf{k}_{2}}^{{\dagger}}\hat{a}_{\mathbf{k}_{3}}\hat{b}_{\mathbf{k}_{4}}\Bigg{\}}.$		(53)

Let us derive the interaction terms in the vicinity of the Dirac points and neglect the gradient terms. This approximation significantly simplifies the final expression. Furthermore, these terms are expected to be weak for small deviations from the Dirac nodes. We use

	$\displaystyle\gamma_{\delta\mathbf{k}}$	$\displaystyle\approx$	$\displaystyle 3+\mathcal{O}(\delta k^{2}),$		(54)
	$\displaystyle\gamma_{\mathbf{k}_{D}^{\zeta}+\delta\mathbf{k}}$	$\displaystyle\approx$	$\displaystyle-\frac{3}{2}a\left(\zeta\delta k_{x}-i\delta k_{y}\right)+\mathcal{O}(\delta k^{2}).$		(55)

Further, we employ the Fourier transform of the operators

	$\displaystyle\hat{a}_{\zeta,\delta\mathbf{k}}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{N}}\int d\mathbf{r}e^{-i\mathbf{r}\cdot\delta\mathbf{k}}\hat{a}_{\zeta,\mathbf{r}},$		(56)
	$\displaystyle\hat{a}_{\zeta,\mathbf{r}}$	$\displaystyle=$	$\displaystyle\frac{1}{\sqrt{N}}\int\frac{d\delta\mathbf{k}}{(2\pi)^{2}}e^{i\mathbf{r}\cdot\delta\mathbf{k}}\hat{a}_{\zeta,\delta\mathbf{k}}.$		(57)

Note also that

\delta_{\mathbf{k}_{1},\mathbf{k}_{2}}=\frac{(2\pi)^{2}}{N}\delta(\mathbf{k}_{1}-\mathbf{k}_{2}).

(58)

Then, the leading order interaction term is

	$\displaystyle H_{4}$	$\displaystyle\simeq$	$\displaystyle-\frac{3J}{N}\sum_{\mathbf{k}_{1},\mathbf{k}_{2},\mathbf{k}_{3},\mathbf{k}_{4}}\sum_{\zeta_{1},\zeta_{2},\zeta_{3},\zeta_{4}}\delta_{\zeta_{2},\zeta_{4}}\hat{a}_{\zeta_{1},\delta\mathbf{k}_{1}}^{{\dagger}}\hat{b}_{\zeta_{2},\delta\mathbf{k}_{2}}^{{\dagger}}\hat{a}_{\zeta_{3},\delta\mathbf{k}_{3}}b_{\zeta_{4},\delta\mathbf{k}_{4}}\frac{(2\pi)^{2}}{N}\int d\mathbf{r}e^{-i\left(\delta\mathbf{k}_{1}+\delta\mathbf{k}_{2}-\delta\mathbf{k}_{4}-\delta\mathbf{k}_{3}\right)\cdot\mathbf{r}}$		(59)
		$\displaystyle=$	$\displaystyle-3J(2\pi)^{2}\int d\mathbf{r}\sum_{\zeta_{1},\zeta_{2}}\hat{a}_{\zeta_{1}}^{{\dagger}}\hat{b}_{\zeta_{2}}^{{\dagger}}\hat{a}_{\zeta_{1}}\hat{b}_{\zeta_{2}}.$		(59)

Therefore, by using the Bogolyubov approximation for the field operators Pitaevskii-Stringari:book , e.g., $\hat{a}\to\psi_{a}$ , in

	$\displaystyle H$	$\displaystyle=$	$\displaystyle\int d\mathbf{r}\Bigg{\{}3JS\sum_{\zeta}\left[\hat{a}_{\zeta}^{{\dagger}}\hat{a}_{\zeta}+\hat{b}_{\zeta}^{{\dagger}}\hat{b}_{\zeta}\right]-iv\sum_{\zeta}\hat{a}_{\zeta}^{{\dagger}}\left[\zeta\partial_{x}-i\partial_{y}\right]\hat{b}_{\zeta}-iv\sum_{\zeta}\hat{b}_{\zeta}^{{\dagger}}\left[\zeta\partial_{x}+i\partial_{y}\right]\hat{a}_{\zeta}-3J(2\pi)^{2}\sum_{\zeta_{1},\zeta_{2}}\hat{a}_{\zeta_{1}}^{{\dagger}}\hat{b}_{\zeta_{2}}^{{\dagger}}\hat{a}_{\zeta_{1}}\hat{b}_{\zeta_{2}}$		(60)
		$\displaystyle+$	$\displaystyle\Delta\sum_{\zeta}\left[\hat{a}_{\zeta}^{{\dagger}}\hat{a}_{\zeta}-\hat{b}_{\zeta}^{{\dagger}}\hat{b}_{\zeta}\right]\Bigg{\}},$		(60)

we obtain the final expression for the energy density

\epsilon=c_{0}n_{\rm tot}-iv\sum_{\zeta}\psi_{a,\zeta}^{*}\left[\zeta\partial_{x}-i\partial_{y}\right]\psi_{b,\zeta}-iv\sum_{\zeta}\psi_{b,\zeta}^{*}\left[\zeta\partial_{x}+i\partial_{y}\right]\psi_{a,\zeta}+g_{ab}n_{a}n_{b}+\Delta(n_{a}-n_{b}).

(61)

Here, $n_{a,\zeta}=\psi_{a,\zeta}=^{*}\psi_{a,\zeta}=$ , $n_{b,\zeta}=\psi_{b,\zeta}^{*}\psi_{b,\zeta}$ , $n_{a}=\sum_{\zeta=\pm}n_{a,\zeta}$ , $n_{b}=\sum_{\zeta=\pm}n_{b,\zeta}$ , and $n_{\rm tot}=n_{a}+n_{b}$ . This energy density serves as a starting point for the free energy in Eq. (1). It is worth noting also that we considered a simple model for Dirac magnons in the derivation above. Other interaction terms such as terms corresponding to the Dzyaloshinskii–Moriya and Kitaev interactions, were omitted in this schematic derivation. However, their effects can be still included phenomenologically, at least in part, by adding additional terms to the free energy (1). The corresponding microscopical analysis will be reported elsewhere.

Appendix B Classification of non-trivial phases of Dirac BEC

In this appendix, we discuss possible non-trivial phases of the Dirac BECs formed at the two valleys of the Brillouin zone. In the absence of interaction, the symmetry class for the BEC corresponds to $U(1)_{a+}\times U(1)_{a-}\times U(1)_{b+}\times U(1)_{b-}$ , where each of the symmetries is related to the independent phase rotation of the wave function. We denote the wave-functions $\psi_{\alpha,\zeta}$ ( $\alpha=a,b,\zeta=\pm$ ) at each valley with the pseudospin components as $\psi_{\alpha,\zeta}=\sqrt{n_{\alpha,\zeta}}e^{i\theta_{\alpha,\zeta}}$ . Interaction could lead to phase-locked condensates, a few of which are provided in Table 1.

Table 1: Classification of Dirac condensates residing in two valleys of the spectrum shown in Fig. 1 based on their phase locking.

Density	Type	Phase		Designation
$n_{a,+}=n_{a,-}$	Symmetric	$\theta_{a,+}=\theta_{a,-}$	$\theta_{b,+}=\theta_{b,-}$	Corotating
		$\theta_{a,+}=-\theta_{a,-}$	$\theta_{b,+}=-\theta_{b,-}$	Anti corotating
$n_{b,+}=n_{b,-}$		$\theta_{a,+}=-\theta_{a,-}$	$\theta_{b,+}=\theta_{b,-}$	Corotating “B” condensate
$n_{a,+}\neq n_{a,-}$	Asymmetric			Anti corotating “A” condensate
		$\theta_{a,+}=\theta_{a,-}$	$\theta_{b,+}=-\theta_{b,-}$	Corotating “A” condensate
$n_{b,+}\neq n_{b,-}$				Anti corotating “B” condensate

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 $\theta_{a/b,+}=\theta_{a/b,-}$ , we call them corotating. In the other case, when the phases are equal in magnitude but differ in sign, the condensates are called anti corotating. However, for suitable interaction channels between the multicomponent condensates, a completely different phase locking can appear as $\theta_{a,+}=-\theta_{a,-}$ and $\theta_{b,+}=\theta_{b,-}$ (or $\theta_{a,+}=\theta_{a,-}$ and $\theta_{b,+}=-\theta_{b,-}$ ). In this case, we define it as anti corotating “A” (“B”) and corotating “B” (“A”).

Appendix C Effect of intercondensate interactions on the condensate dynamics

Let us consider the effect of nonzero interaction constant $g_{ab}$ on the condensate dynamics. By using the model defined in Sec. III.1 in the main text, we numerically calculate the magnon densities $n_{a}$ , $n_{b}$ , and $n_{3}$ for $g_{a}=g_{b}=g_{ab}=t_{0}^{-1}$ (see Fig. 10) as well as $g_{a}=g_{b}=0$ , and $g_{ab}=-t_{0}^{-1}$ (see Fig. 11). Note that the interaction constants $g_{a}$ , $g_{b}$ , $g_{ab}$ are not manifested in the condensate dynamics when the phase-mixing terms are absent (e.g., the Josephson coupling term $u_{ab}$ ). As one can see by comparing Fig. 4 with Figs. 10 and 11, the effect of the intercondensate interaction constant $g_{ab}$ is similar to that of $g_{a}$ and $g_{b}$ at small $u_{ab}$ . On the other hand, the anharmonic oscillations of the condensates are more easily accessible for a nonzero $g_{ab}$ .

Appendix D Four-component model

To describe the magnon BECs at the Dirac nodes and at the $\Gamma$ point, we extend the phenomenological model in Sec. III.1. In addition, to Eqs. (5)–(7), a Gross-Pitaevskii equation for the magnon population at $k=0$ is included. The full system then reads as

	$\displaystyle i\partial_{t}\psi_{a}+i\Gamma_{a}\psi_{a}=\left(c_{0}-\mu+\Delta\right)\psi_{a}-iv\left[\left(\partial_{x}-i\partial_{y}\right)\psi_{b}\right]+g_{a}n_{a}\psi_{a}+g_{ab}n_{b}\psi_{a}+u_{ab}\psi_{a}^{*}\left(\psi_{b}\right)^{2}$
	$\displaystyle+iP_{a3}\Gamma_{a}n_{3}\psi_{a}-iP_{a0}\Gamma_{0}n_{0}\psi_{a},$		(62)
	$\displaystyle i\partial_{t}\psi_{b}+i\Gamma_{b}\psi_{b}=\left(c_{0}-\mu-\Delta\right)\psi_{b}-iv\left[\left(\partial_{x}+i\partial_{y}\right)\psi_{a}\right]+g_{b}n_{b}\psi_{b}+g_{ab}n_{a}\psi_{b}+u_{ab}^{}\psi_{b}^{}\left(\psi_{a}\right)^{2}$
	$\displaystyle+iP_{b3}\Gamma_{b}n_{3}\psi_{b}-iP_{b0}\Gamma_{0}n_{0}\psi_{b},$		(63)
	$\displaystyle i\partial_{t}\psi_{0}+i\Gamma_{0}\psi_{0}=-\frac{1}{2m_{\rm eff}}\left(\partial_{x}^{2}+\partial_{y}^{2}\right)\psi_{0}+g_{0}n_{0}\psi_{0}+iP_{0b}\Gamma_{a}n_{a}\psi_{0}+iP_{0a}\Gamma_{b}n_{b}\psi_{0},$		(64)
	$\displaystyle\partial_{t}n_{3}+\Gamma_{3}n_{3}=\Gamma_{3}\tilde{P}(t)-\sum_{\zeta=\pm}\left(P_{3a}\Gamma_{a}n_{a,\zeta}+P_{3b}\Gamma_{b}n_{b,\zeta}\right)n_{3}.$		(65)

Here, the magnon BEC at the $\Gamma$ point is described by Eq. (64). Due to kinematic constraints, it is unlikely that magnons can scatter to the $\Gamma$ point directly from the pumped population. On the other hand, the lowest branch can acquire magnons via thermalization of magnons at the Dirac points. For example, this could be achieved via four-magnon processes (in the case of YIG, see, e.g., Ref. Serga-Hillebrands:2014 ). This process is described via the terms with $P_{0b}$ and $P_{0a}$ .

By using the values of parameters in Sec. III.2 as well as assuming that $\Gamma_{0}=\Gamma_{a}=\Gamma_{b}=\Gamma_{3}=0.25\,t_{0}^{-1}$ and $P_{0a}=P_{0b}=0.5$ , we present the results for the population densities in Fig. 12. As one can see, the magnon density at the $\Gamma$ point raises at the later stages of the pumping. Since the source of these magnons is the magnon BEC at the Dirac points, the rise of the magnon density at the lowest branch leads to the depletion of the Dirac magnons. Another effect that would be interesting to test experimentally and/or verify microscopically is the rise of the pumped magnon density. In the model at hand, this fact stems from the decrease of the Dirac magnons densities.

Appendix E 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 do 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_{0}^{-1}$ , $g_{ab}=0$ , and the ground state defined in Sec. IV.1. The corresponding results for $u_{ab}=0.2\,t_{0}^{-1}$ are shown in Figs. 13 and 14 at $n_{b}\neq 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.3. Moreover, since we no longer have independent phases for the ground state wave functions at $n_{b}\neq 0$ , one of the modes can become gapped (see Fig. 13(a)).

References

(1) M. I. Katsnelson, Graphene: Carbon in two dimensions, Mater. Today 10, 20 (2007).
(2) A. K. Geim and K. S. Novoselov, The rise of graphene, Nat. Mater. 6, 183 (2007).
(3) A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
(4) A. K. Geim, Graphene: Status and prospects, Science 324, 1530 (2009).
(5) M. I. Katsnelson, Graphene: Carbon in two dimensions (Cambridge University Press, 2012).
(6) M. Z. Hasan and C. L. Kane, Colloquium: Topological insulators, Rev. Mod. Phys. 82, 3045 (2010).
(7) X.-L. Qi and S.-C. Zhang, Topological insulators and superconductors, Rev. Mod. Phys. 83, 1057 (2011).
(8) B. A. Bernevig and T. L. Hughes, Topological insulators and topological superconductors (Princeton University Press, 2013).
(9) A. M. Turner and A. Vishwanath, Beyond band insulators: Topology of semimetals and interacting phases (Elsevier Science, 2013).
(10) T. O. Wehling, A. M. Black-Schaffer, and A. V. Balatsky, Dirac materials, Adv. Phys. 63, 1 (2014).
(11) B. Yan and C. Felser, Topological materials: Weyl semimetals, Annu. Rev. Condens. Matter Phys. 8, 337 (2017).
(12) M. Z. Hasan, S.-Y. Xu, I. Belopolski, and C.-M. Huang, Discovery of Weyl fermion semimetals and topological Fermi arc states, Annu. Rev. Condens. Mattter Phys. 8, 289 (2017).
(13) N. P. Armitage, E. J. Mele, and A. Vishwanath, Weyl and Dirac semimetals in three-dimensional solids, Rev. Mod. Phys. 90, 015001 (2018).
(14) T. Jacqmin, I. Carusotto, I. Sagnes, M. Abbarchi, D. D. Solnyshkov, G. Malpuech, E. Galopin, A. Lemaître, J. Bloch, and A. Amo, Direct observation of Dirac cones and a flatband in a honeycomb lattice for polaritons, Phys. Rev. Lett. 112, 116402 (2014).
(15) L. Lu, Z. Wang, D. Ye, L. Ran, L. Fu, J. D. Joannopoulos, M. Soljačić, Experimental observation of Weyl points, Science 349, 622 (2015).
(16) H.-X. Wang, L. Xu, H.-Y. Chen, and J.-H. Jiang, Three-dimensional photonic Dirac points stabilized by point group symmetry, Phys. Rev. B 93, 235155 (2016).
(17) B. Yang, Q. Guo, B. Tremain, R. Liu, L. E. Barr, Q. Yan, W. Gao, H. Liu, Y. Xiang, J. Chen, C. Fang, A. Hibbins, L. Lu, S. Zhang, Ideal Weyl points and helicoid surface states in artificial photonic crystal structures, Science 359, 1013 (2018).
(18) Q. Guo, O. You, B. Yang, J. B. Sellman, E. Blythe, H. Liu, Y. Xiang, J. Li, D. Fan, J. Chen, C. T. Chan, and S. Zhang, Observation of three-dimensional photonic Dirac points and spin-polarized surface arcs, Phys. Rev. Lett. 122, 203903 (2019).
(19) Y. Yang, Z. Gao, H. Xue, L. Zhang, M. He, Z. Yang, R. Singh, Y. Chong, B. Zhang, and H. Chen, Realization of a three-dimensional photonic topological insulator, Nature 565, 622 (2019).
(20) M. Xiao, W.-J. Chen, W.-Y. He, and C. T. Chan, Synthetic gauge flux and Weyl points in acoustic systems, Nat. Phys. 11, 920 (2015).
(21) M. Serra-Garcia, V. Peri, R. Süsstrunk, O. R. Bilal, T. Larsen, L. G. Villanueva, and S. D. Huber, Observation of a phononic quadrupole topological insulator, Nature 555, 342 (2018).
(22) X. Cai, L. Ye, C. Qiu, M. Xiao, R. Yu, M. Ke, and Z. Liu, Symmetry-enforced three-dimensional Dirac phononic crystals, Light Sci. Appl. 9, 38 (2020).
(23) S. S. Pershoguba, S. Banerjee, J. C. Lashley, J. Park, H. Ågren, G. Aeppli, and A. V. Balatsky, Dirac magnons in honeycomb ferromagnets, Phys. Rev. X 8, 011010 (2018).
(24) L. Chen, J.-H. Chung, B. Gao, T. Chen, M. B. Stone, A. I. Kolesnikov, Q. Huang, and P. Dai, Topological spin excitations in honeycomb ferromagnet CrI₃, Phys. Rev. X 8, 041028 (2018).
(25) I. Lee, F. G. Utermohlen, D. Weber, K. Hwang, C. Zhang, J. van Tol, J. E. Goldberger, N. Trivedi, and P. C. Hammel, Fundamental spin interactions underlying the magnetic anisotropy in the Kitaev ferromagnet CrI₃, Phys. Rev. Lett. 124, 017201 (2020).
(26) S. Banerjee, J. Fransson, A. M. Black-Schaffer, H. Ågren, and A. V. Balatsky, Granular superconductor in a honeycomb lattice as a realization of bosonic Dirac material, Phys. Rev. B 93, 134502 (2016).
(27) J. Fransson, A. M. Black-Schaffer, and A. V. Balatsky, Magnon Dirac materials, Phys. Rev. B 94, 075401 (2016).
(28) C. Gong, L. Li, Z. Li, H. Ji, A. Stern, Y. Xia, T. Cao, W. Bao, C. Wang, Y. Wang, Z. Q. Qiu, R. J. Cava, S. G. Louie, J. Xia, and X. Zhang, Discovery of intrinsic ferromagnetism in two-dimensional van der Waals crystals, Nature 546, 265 (2017).
(29) B. Huang, G. Clark, E. Navarro-Moratalla, D. R. Klein, R. Cheng, K. L. Seyler, D. Zhong, E. Schmidgall, M. A. McGuire, D. H. Cobden, W. Yao, D. Xiao, P. Jarillo-Herrero, and X. Xu, Layer-dependent ferromagnetism in a van der Waals crystal down to the monolayer limit, Nature 546, 270 (2017).
(30) A. C. Gossard, V. Jaccarino, and J. P. Remeika, Experimental test of the spin-wave theory of a ferromagnet, Phys. Rev. Lett. 7, 122 (1961).
(31) H. L. Davis and A. Narath, Spin-wave renormalization applied to ferromagnetic CrBr₃, Phys. Rev. 134, A433 (1964).
(32) F.-Y. Li, Y.-D. Li, Y. B. Kim, L. Balents, Y. Yu, and G. Chen, Weyl magnons in breathing pyrochlore antiferromagnets, Nat. Commun. 7, 12691 (2016).
(33) A. Mook, J. Henk, and I. Mertig, Tunable magnon Weyl points in ferromagnetic pyrochlores, Phys. Rev. Lett. 117, 157204 (2016).
(34) S. A. Owerre, Weyl magnons in noncoplanar stacked kagome antiferromagnets, Phys. Rev. B 97, 094412 (2018).
(35) W. Yao, C. Li, L. Wang, S. Xue, Y. Dan, K. Iida, K. Kamazawa, K. Li, C. Fang, and Y. Li, Topological spin excitations in a three-dimensional antiferromagnet, Nat. Phys. 14, 1011 (2018).
(36) L.-C. Zhang, Y. A. Onykiienko, P. M. Buhl, Y. V. Tymoshenko, P. Čermák, A. Schneidewind, J. R. Stewart, A. Henschel, M. Schmidt, S. Blügel, D. S. Inosov, and Y. Mokrousov, Magnonic Weyl states in Cu₂OSeO₃, Phys. Rev. Research 2, 013063 (2020).
(37) S. Banerjee, D. Abergel, H. Agren, G. Aeppli, and A. V. Balatsky, Interacting Dirac materials, J. Phys. Condens. Matter 32, 405603 (2020).
(38) L. Zhang, J. Ren, J.-S. Wang, and B. Li, Topological magnon insulator in insulating ferromagnet, Phys. Rev. B 87, 144101 (2013).
(39) A. Mook, J. Henk, and I. Mertig, Edge states in topological magnon insulators, Phys. Rev. B 90, 024412 (2014).
(40) S. A. Owerre, A first theoretical realization of honeycomb topological magnon insulator, J. Phys. Condens. Matter 28, 386001 (2016).
(41) K. Nakata, J. Klinovaja, and D. Loss, Magnonic quantum Hall effect and Wiedemann-Franz law, Phys. Rev. B 95, 125429 (2017).
(42) A. A. Kovalev, V. A. Zyuzin, and B. Li, Pumping of magnons in a Dzyaloshinskii-Moriya ferromagnet, Phys. Rev. B 95, 165106 (2017).
(43) S. A. Owerre, Magnonic analogs of topological Dirac semimetals, J. Phys. Commun. 1, 025007 (2017).
(44) X. S. Wang, A. Brataas, and R. E. Troncoso, Bosonic Bott index and disorder-induced topological transitions of magnons, Phys. Rev. Lett. 125, 217202 (2020).
(45) Y. Ferreiros and M. A. H. Vozmediano, Elastic gauge fields and Hall viscosity of Dirac magnons, Phys. Rev. B 97, 054404 (2018).
(46) M. M. Nayga, S. Rachel, and M. Vojta, Magnon Landau levels and emergent supersymmetry in strained antiferromagnets, Phys. Rev. Lett. 123, 207204 (2019).
(47) T. Liu and Z. Shi, Twist-induced magnon Landau levels in honeycomb magnets, arXiv:2004.01371 (2020).
(48) P. A. McClarty and J. G. Rau, Non-Hermitian topology of spontaneous magnon decay, Phys. Rev. B 100, 100405(R) (2019).
(49) A. Pertsova and A. V. Balatsky, Dynamically induced excitonic instability in pumped Dirac materials, Ann. Phys. (Berlin) 532, 1900549 (2020).
(50) F. Dalfovo, S. Giorgini, L. P. Pitaevskii, and S. Stringari, Theory of Bose-Einstein condensation in trapped gases, Rev. Mod. Phys. 71, 463 (1999).
(51) R. Ozeri, N. Katz, J. Steinhauer, and N. Davidson, Colloquium: Bulk Bogoliubov excitations in a Bose-Einstein condensate, Rev. Mod. Phys. 77, 187 (2005).
(52) L. P. Pitaevskii and S. Stringari, Bose-Einstein condensation and superfluidity (Oxford University Press, 2016).
(53) L. H. Haddad and L. D. Carr, The nonlinear Dirac equation in Bose–Einstein condensates: Foundation and symmetries, Physica D 238, 1413 (2009).
(54) L. H. Haddad and L. D. Carr, Relativistic linear stability equations for the nonlinear Dirac equation in Bose-Einstein condensates, EPL 94, 56002 (2011).
(55) L. H. Haddad, K. M. O’Hara, and L. D. Carr, Nonlinear Dirac equation in Bose-Einstein condensates: Preparation and stability of relativistic vortices, Phys. Rev. A 91, 043609 (2015).
(56) L. H. Haddad, C. M. Weaver, L. D. Carr, The nonlinear Dirac equation in Bose–Einstein condensates: I. Relativistic solitons in armchair nanoribbon optical lattices, New J. Phys. 17, 063033 (2015).
(57) L. H. Haddad and L. D. Carr, The nonlinear Dirac equation in Bose–Einstein condensates: II. Relativistic soliton stability analysis, New J. Phys. 17, 063034 (2015).
(58) L. H. Haddad and L. D. Carr, The nonlinear Dirac equation in Bose–Einstein condensates: vortex solutions and spectra in a weak harmonic trap, New J. Phys. 17, 113011 (2015).
(59) L. H. Haddad and L. D. Carr, The nonlinear Dirac equation in Bose–Einstein condensates: Superfluid fluctuations and emergent theories from relativistic linear stability equations, New J. Phys. 17, 093037 (2015).
(60) J. Cuevas-Maraver, P. G. Kevrekidis, A. Saxena, A. Comech, and R. Lan, Stability of solitary waves and vortices in a 2D nonlinear Dirac model, Phys. Rev. Lett. 116, 214101 (2016).
(61) V. Carmona, J. Cuevas-Maraver, F. Fernández-Sánchez, and E. García- Medina (eds.), Nonlinear systems, Vol. 1. Understanding complex systems. (Springer, Cham, 2018).
(62) A. N. Poddubny and D. A. Smirnova, Ring Dirac solitons in nonlinear topological systems, Phys. Rev. A 98, 013827 (2018).
(63) S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, A. A. Serga, B. Hillebrands, and A. N. Slavin, Bose–Einstein condensation of quasi-equilibrium magnons at room temperature under pumping, Nature 443, 430 (2006).
(64) V. E. Demidov, O. Dzyapko, S. O. Demokritov, G. A. Melkov, and A. N. Slavin, Thermalization of a parametrically driven magnon gas leading to Bose-Einstein condensation, Phys. Rev. Lett. 99, 037205 (2007).
(65) O. Dzyapko, V. E. Demidov, S. O. Demokritov, G. A. Melkov, and A. N. Slavin, Direct observation of Bose-Einstein condensation in a parametrically driven gas of magnons, New J. Phys. 9, 64 (2007).
(66) V. E. Demidov, O. Dzyapko, S. O. Demokritov, G. A. Melkov, and A. N. Slavin, Observation of spontaneous coherence in Bose-Einstein condensate of magnons, Phys. Rev. Lett. 100, 047205 (2008).
(67) P. Nowik-Boltyk, O. Dzyapko, V. E. Demidov, N. G. Berloff, and S. O. Demokritov, Spatially non-uniform ground state and quantized vortices in a two-component Bose-Einstein condensate of magnons, Sci. Rep. 2, 482 (2012).
(68) A. A. Serga, V. S. Tiberkevich, C. W. Sandweg, V. I. Vasyuchka, D. A. Bozhko, A. V. Chumak, T. Neumann, B. Obry, G. A. Melkov, A. N. Slavin, and B. Hillebrands, Bose–Einstein condensation in an ultra-hot gas of pumped magnons, Nat. Commun. 5, 3452 (2014).
(69) D. A. Bozhko, A. J. E. Kreil, H. Yu. Musiienko-Shmarova, A. A. Serga, A. Pomyalov, V. S. L’vov, and B. Hillebrands, Bogoliubov waves and distant transport of magnon condensate at room temperature, Nat. Commun. 10, 2460 (2019).
(70) M. Mohseni, A. Qaiumzadeh, A. A. Serga, A. Brataas, B. Hillebrands, and P. Pirro, Bose-Einstein condensation of nonequilibrium magnons in laterally confined systems, New J. Phys. 22 083080 (2020)
(71) M. Mohseni, Q. Wang, B. Heinz, M. Kewenig, M. Schneider, F. Kohl, B. Lägel, C. Dubs, A. V. Chumak, and P. Pirro, Controlling of nonlinear relaxation of quantized magnons in nano-devices, arXiv:2006.03400 (2020).
(72) J. Kasprzak, M. Richard, S. Kundermann, A. Baas, P. Jeambrun, J. M. J. Keeling, F. M. Marchetti, M. H. Szymańska, R. André, J. L. Staehli, V. Savona, P. B. Littlewood, B. Deveaud, and L. S. Dang, Bose–Einstein condensation of exciton polaritons, Nature 443, 409 (2006).
(73) H. Deng, H. Haug, and Y. Yamamoto, Exciton-polariton Bose-Einstein condensation, Rev. Mod. Phys. 82, 1489 (2010).
(74) B. Deveaud, Exciton-polariton Bose-Einstein condensates, Annu. Rev. Condens. Matter Phys. 6 155 (2015)
(75) Y. Sun, P. Wen, Y. Yoon, G. Liu, M. Steger, L. N. Pfeiffer, K. West, D. W. Snoke, and K. A. Nelson, Bose-Einstein condensation of long-lifetime polaritons in thermal equilibrium, Phys. Rev. Lett. 118, 016602 (2017).
(76) N. P. Proukakis, D. W. Snoke, P. B. Littlewood (Eds.), Universal themes of Bose-Einstein condensation (Cambridge University Press, 2017).
(77) V. S. L’vov Wave turbulence under parametric excitation: Applications to magnets (Springer, Berlin Heidelberg, 1994).
(78) E. V. Gorbar, V. A. Miransky, and I. A. Shovkovy, Chiral asymmetry of the Fermi surface in dense relativistic matter in a magnetic field, Phys. Rev. C 80, 032801(R) (2009).
(79) Y. J. Lin, K. Jimenez-Garcia, and I. B. Spielman, Spin–orbit-coupled Bose–Einstein condensates, Nature 471, 83 (2011).
(80) J. L. Mañes, F. de Juan, M. Sturla, and M. A. H. Vozmediano, Generalized effective Hamiltonian for graphene under nonuniform strain, Phys. Rev. B 88, 155405 (2013).
(81) V. Hahn and P. Kopietz, Collisionless kinetic theory for parametrically pumped magnons, Eur. Phys. J. B 93, 132 (2020).
(82) S. O. Demokritov, V. E. Demidov, O. Dzyapko, G. A. Melkov, and A. N. Slavin, Quantum coherence due to Bose-Einstein condensation of parametrically driven magnons, New J. Phys. 10, 045029 (2008).
(83) S. M. Rezende, Theory of coherence in Bose-Einstein condensation phenomena in a microwave-driven interacting magnon gas, Phys. Rev. B 79, 174411 (2009).
(84) S. M. Rezende, Wave function of a microwave-driven Bose-Einstein magnon condensate, Phys. Rev. B 81, 020414(R) (2010).
(85) D. Prananto, Y. Kainuma, K. Hayashi, N. Mizuochi, K. Uchida, and T. An, Probing thermal magnon current via nitrogen-vacancy centers in diamond, arXiv:2007.13433 (2020).
(86) V. V. Kruglyak, S. O. Demokritov, and D. Grundler, Magnonics, J. Phys. D: Appl. Phys. 43 264001 (2010).
(87) See Supplemental Material for the animation of dynamics of the Dirac magnon populations.
(88) E. Fradkin, S. A. Kivelson, M. J. Lawler, J. P. Eisenstein, and A. P. Mackenzie, Nematic fermi fluids in condensed matter physics, Annu. Rev. Condens. Matter Phys. 1, 153 (2010).
(89) C.-M. Jian and H. Zhai, Paired superfluidity and fractionalized vortices in systems of spin-orbit coupled bosons, Phys. Rev. B 84, 060508(R) (2011).
(90) M. Larsson and A. Balatsky, Paul Dirac and the Nobel Prize in physics, Phys. Today 72(11), 46 (2019).