\section{Introduction}
\label{sec: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~\cite{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~\cite{Jacqmin-Amo:2014,Lu-Soljacic:2015,Wang-Jiang:2016,Yang-Zhang:2018,Guo-Zhang:2018,Yang-Chen:2019} and phononic~\cite{Xiao-Chan:2015,Serra-Garcia-Huber:2018,Cai-Liu:2020} crystals and magnets~\cite{Pershoguba-Balatsky:2018,Chen-Dai:2018,Lee-Hammel:2020}.
Honeycomb arrays of superconducting grains were also predicted to host bosonic Dirac excitations~\cite{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.~\cite{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~\cite{Fransson-Balatsky:2016,Pershoguba-Balatsky:2018}.
The corresponding magnon Dirac materials~\cite{Fransson-Balatsky:2016} such as van der Waals crystals Cr$_2$Ge$_2$Te$_6$ ~\cite{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$_3$~\cite{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~\cite{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~\cite{Li-Chen:2016,Mook-Mertig:2016,Owerre:2018}. Linear crossing in the magnon spectrum was indeed observed in the antiferromagnet Cu$_3$TeO$_6$~\cite{Yao-Li:2018} and the multiferroic ferrimagnet Cu$_2$OSeO$_3$~\cite{Zhang-Mokrousov:2020}.

A wide range of properties was observed within this emergent class of Dirac magnon materials. We mention
interaction effects~\cite{Pershoguba-Balatsky:2018,Banerjee-Balatsky:2020} that renormalize the slope of the dispersion, topological properties~\cite{Zhang-Li:2013,Mook-Mertig:2014,Owerre:2016,Nakata-Loss:2017,Kovalev-Li:2017,Owerre:2017,Zhang-Mokrousov:2020,Wang-Troncoso:2020}, strain effects~\cite{Ferreiros-Vozmediano:2018,Nayga-Vojta:2019,Liu-Shi:2020}, and non-Hermitian dynamics~\cite{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 {\em 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.~\cite{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~\cite{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.~\cite{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.~\cite{Cuevas-Maraver-Lan:2016,Cuevas-Saxena:2018,Poddubny-Smirnova:2018}). Among the most interesting features, we notice various types of vortices and solitons~\cite{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~\cite{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.~\cite{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~\cite{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.~\ref{sec:Model}. 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.~\ref{sec:phenomenology-pumping}. Section~\ref{sec:phenomenology-collective} is devoted to collective modes in the Dirac BECs. The results are discussed and summarized in Sec.~\ref{sec:Summary}. A schematic derivation of the free energy ansatz used in the main text is given in Appendix~\ref{sec:App-magnon-H}. A few non-trivial phases of the Dirac BEC are summarized in Appendix~\ref{sec:App-BEC-Class}. The additional results for the population dynamics at different interaction constants are presented in Appendix~\ref{sec:App-BEC-gab}. Appendix~\ref{sec:App-Josephson-omega} 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~\ref{sec:App-phenomenology-4component}. Through this study, we set $\hbar=k_{B}=1$.

\section{Model}
\label{sec:Model}

In this section, we present a phenomenological approach to Dirac BECs. The approach is similar to that for conventional BECs~\cite{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~\ref{sec:App-magnon-H}):
\begin{eqnarray}
\label{phen-F-def}
F &=& (c_0-\mu) n_{\rm tot}+\Delta  \sum_{\zeta=\pm}(n_{a,\zeta}-n_{b,\zeta}) -iv \Psi^{\dag}(\tau_z\otimes\bm{\sigma}\cdot\bm{\nabla})\Psi \nonumber\\
&+& \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).
\end{eqnarray}
The four-component wave function of Dirac BEC is
\begin{equation}
\label{model-bispinor-def}
\Psi= \left(\psi_{a,+},\psi_{b,+},\psi_{b,-},\psi_{a,-}\right)^{T},
\end{equation}
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~\cite{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~\cite{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~\ref{sec:App-magnon-H} for the discussion of the Dirac spin-wave dispersion relation) and BECs corresponding to different valleys and pseudospins are shown schematically in Fig.~\ref{fig:phen-spectrum-schematic}. 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.~\ref{fig:phen-spectrum-schematic}(b), the presence of the $\Delta$ term opens the gap in the spectrum.

\begin{figure*}[!ht]
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig1a.eps}}
\hspace{0.1\textwidth}
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig1b.eps}}
\caption{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\neq0$ 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.~\ref{sec:phenomenology-pumping} for the description of the pumping mechanism). Finally, $n_0$ denotes the population at the $\Gamma$ point (see also Appendix~\ref{sec:App-phenomenology-4component}).
}
\label{fig:phen-spectrum-schematic}
\end{figure*}

By varying the free energy (\ref{phen-F-def}) 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:
\begin{eqnarray}
\label{phen-GPE-a-def}
i\partial_t \psi_{a,\zeta} &=& \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},\\
\label{phen-GPE-b-def}
i\partial_t \psi_{b,\zeta} &=& \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}.
\end{eqnarray}
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~\ref{sec:App-BEC-Class}.

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.

\section{Generation of Dirac Bose-Einstein condensates via pumping}
\label{sec:phenomenology-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.~\cite{Pertsova-Balatsky:2020}.) As we discussed in the introduction, Dirac magnons can be realized in, e.g., CrI$_3$ (see Refs.~\cite{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.~\cite{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$_3$, see, e.g., Ref.~\cite{Lee-Hammel:2020}.

\subsection{Three-component model}
\label{sec:phenomenology-3component}

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.~\ref{fig:phen-spectrum-schematic}). The extension of this model to the case where the magnon BEC is formed at the $\Gamma$ point is presented in Appendix~\ref{sec:App-phenomenology-4component}. These populations could be achieved by the parametric pumping similarly to the case of conventional magnetic materials~\cite{Lvov:book} or any other suitable pumping scheme. The results below do not depend on the details of the pumping. The system reads as
\begin{eqnarray}
\label{phen-GPE-a-1}
&&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},\\
\label{phen-GPE-b-1}
&&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},\\
\label{phen-GPE-pump-1}
&&\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.
\end{eqnarray}
Compared to Eqs.~(\ref{phen-GPE-a-def}) and (\ref{phen-GPE-b-def}), 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 (\ref{phen-GPE-pump-1}) 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~\cite{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
\begin{equation}
\label{phen-t-tP-def}
\tilde{P}(t) = \sqrt{P_{\rm max}}\, \theta(t-t_{\rm Begin}) \theta(t_{\rm End}-t).
\end{equation}
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~\cite{Demidov-Slavin:2008,Demokritov-Slavin:2008} (see also Refs.~\cite{Rezende:2009,Rezende:2010} for the theoretical discussion). While only a simple step-like profile in Eq.~(\ref{phen-t-tP-def}) is considered in this work, the time profiles of magnon BECs can be modulated by applying more complicated time-dependent pump.

\subsection{Time evolution of Dirac magnons}
\label{sec:phenomenology-time}

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:
\begin{eqnarray}
\label{phen-three-t-num-vars-be}
&&\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,\\
&&P_{3a}=P_{3b}=P_{a3}=P_{b3}=1, \quad t_{\rm Begin} = 10\,t_0, \quad t_{\rm End} = 200\,t_0.
\label{phen-three-t-num-vars-ee}
\end{eqnarray}
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.~\cite{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~\ref{sec:App-BEC-gab}. 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.~(\ref{phen-GPE-a-1}) and (\ref{phen-GPE-b-1})]. This follows from the fact that the density dynamics is determined by the imaginary part of the Gross-Pitaevskii equations (\ref{phen-GPE-a-1}) and (\ref{phen-GPE-b-1}).

The time profiles of magnon densities are shown in Figs.~\ref{fig:phen-three-t-P-Delta-I} and \ref{fig:phen-three-t-P-I-density} 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.~\ref{fig:phen-three-t-P-I-density}). 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.~(\ref{phen-GPE-a-1}) through (\ref{phen-GPE-pump-1}), there is a noticeable depletion of the pumped population. After turning off the pump, both condensates and pumped population decay during the time $\tau\sim1/\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.

\begin{figure*}[!ht]
\centering
\subfigure[]{\includegraphics[height=0.28\textwidth]{Fig2a.eps}}
\hspace{0.02\textwidth}
\subfigure[]{\includegraphics[height=0.28\textwidth]{Fig2b.eps}}
\hspace{0.02\textwidth}
\subfigure[]{\includegraphics[height=0.28\textwidth]{Fig2c.eps}}
\caption{The dependence of the magnon densities on time $t$ at a few values of the Haldane gap $\Delta$ for $P_{\rm max}=0.5\,P_0$. The densities are normalized to the total magnon density at $t=t_{\rm End}$. Here, the values of parameters given in Eqs.~(\ref{phen-three-t-num-vars-be}) and (\ref{phen-three-t-num-vars-ee}) are used. In addition, $P_0$ is the characteristic strength of the pump ($P_0=t_0^{-2}$). The Haldane gap leads to the splitting of the magnon densities corresponding to different pseudospins.}
\label{fig:phen-three-t-P-Delta-I}
\end{figure*}

\begin{figure*}[!ht]
\centering
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig3a.eps}}
\hspace{0.1\textwidth}
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig3b.eps}}
\caption{The dependence of the magnon densities (Dirac BEC in panel (a) and pumped population in panel (b)) on time $t$ and pump power $P_{\rm max}$. Here, the values of parameters given in Eqs.~(\ref{phen-three-t-num-vars-be}) and (\ref{phen-three-t-num-vars-ee}) are used. In addition, $P_0$ is the characteristic strength of the pump ($P_0=t_0^{-2}$). The pumped population demonstrates a box-shaped profile. On the other hand, the Dirac BEC forms much faster at larger powers.}
\label{fig:phen-three-t-P-I-density}
\end{figure*}

Further, we investigate the effect of the Haldane gap $\Delta$. The corresponding results are shown in Fig.~\ref{fig:phen-three-t-P-Delta-I} 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.~\ref{fig:phen-three-t-P-Josephson-uab1-uab}. 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.~\ref{fig:phen-three-t-P-Josephson-period}. 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.~\cite{Prananto-An:2020,Kruglyak-Grundler:2010} for reviews).

\begin{figure*}[!ht]
\centering
\subfigure[]{\includegraphics[height=0.27\textwidth]{Fig4a.eps}}
\hspace{0.02\textwidth}
\subfigure[]{\includegraphics[height=0.27\textwidth]{Fig4b.eps}}
\hspace{0.02\textwidth}
\subfigure[]{\includegraphics[height=0.27\textwidth]{Fig4c.eps}}
\caption{The dependence of the magnon densities on time $t$ at a few values of $u_{ab}$ for $P_{\rm max}=0.5\,P_0$ and $\Delta=0.1\,t_0^{-1}$. The densities are normalized to the total magnon density at $t=t_{\rm End}$. Here, the values of parameters given in Eqs.~(\ref{phen-three-t-num-vars-be}) and (\ref{phen-three-t-num-vars-ee}) are used. In addition, $P_0$ is the characteristic strength of the pump ($P_0=t_0^{-2}$). The Rabi oscillations, whose period is determined by $1/\Delta$ at small $u_{ab}$ and large $\Delta$ ($\Delta>0.1\,t_0^{-1}$ at $u_{ab}=0.2\,t_0^{-1}$), are clearly visible for the magnon BECs.}
\label{fig:phen-three-t-P-Josephson-uab1-uab}
\end{figure*}

\begin{figure*}[!ht]
\centering
\includegraphics[height=0.35\textwidth]{Fig5.eps}
\caption{The period of inter-condensate (Rabi) oscillations is shown by the black dots. The red line represents a linear fit. For large values of $\Delta$, we indeed observe the expected for the Rabi oscillations scaling for the period as inverse gap. On the other hand, for smaller values of the Haldane gap, the dynamic is more complicated and does not follow a simple linear behavior. We see a significant scatter in the period due to nonlinear dynamic at $\Delta\lesssim t_0^{-1}$. In this case, other interaction parameters also become important.}
\label{fig:phen-three-t-P-Josephson-period}
\end{figure*}

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~\cite{SM}. The snapshot after the pump is turned off ($t=300\,t_0$) is shown in Fig.~\ref{fig:phen-three-Bloch}.

\begin{figure*}[!ht]
\centering
\includegraphics[height=0.35\textwidth]{Fig6.eps}
\caption{The trajectory of the system in the density space $\left(n_{a},n_{b},n_3\right)/n_{\rm tot}$ 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 $n_3$ 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.~(\ref{phen-three-t-num-vars-be}) and (\ref{phen-three-t-num-vars-ee}) are used. In addition, we set $u_{ab}=0.1\,t_0^{-1}$ and $\Delta=0.025\,t_0^{-1}$. For the full time dependence, see the Supplemental Material~\cite{SM}.}
\label{fig:phen-three-Bloch}
\end{figure*}

\section{Collective modes}
\label{sec:phenomenology-collective}

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~\cite{Dalfovo-Stringari:1999-rev,Ozeri-Davidson:2005-rev,Pitaevskii-Stringari:book} for the magnon BEC. In this case, one looks for the solution as
\begin{eqnarray}
\label{phen-CM-a-uv}
\psi_{a} &=& \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],\\
\label{phen-CM-b-uv}
\psi_{b} &=& \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].
\end{eqnarray}
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.

\subsection{Ground state}
\label{sec:phenomenology-collective-ground}

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.~(\ref{phen-CM-a-uv}) and (\ref{phen-CM-b-uv}), the Gross-Pitaevskii equations (\ref{phen-GPE-a-def}) and (\ref{phen-GPE-b-def}) lead to
\begin{eqnarray}
\label{PCG-GPE-a}
&&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,\\
\label{PCG-GPE-b}
&&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.
\end{eqnarray}
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.~(\ref{PCG-GPE-a}) and (\ref{PCG-GPE-b}), 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
\begin{equation}
\label{phen-coll-J-phase-uab1}
\theta_{u_{ab}}= 2(\theta_{a}-\theta_b) +\pi l, \quad l=0,1,2,3,\ldots,
\end{equation}
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
\begin{eqnarray}
\label{P-C-ground-mu-1}
\mu &=& c_0 +\Delta +g_{a}n_{a} +g_{ab} n_b +u_{ab}n_b,\\
\label{P-C-ground-mu-2}
\mu &=& c_0 -\Delta +g_{b}n_{b} +g_{ab} n_a +u_{ab}n_a.
\end{eqnarray}
Note that even and odd values of $l$ correspond to positive and negative values  of $u_{ab}$ in the equation above.
The system (\ref{P-C-ground-mu-1}) and (\ref{P-C-ground-mu-2}) has a solution for
\begin{equation}
\label{P-C-ground-nb}
n_{b} = \frac{2\Delta +\left(g_{a} -g_{ab} -u_{ab}\right)n_a}{g_{b}-g_{ab} -u_{ab}}.
\end{equation}
(One can fix $n_b$ and solve for $n_a$, which is physically equivalent). The effective chemical potential for the density (\ref{P-C-ground-nb}) reads as
\begin{eqnarray}
\label{P-C-ground-mu-1-nbneq0}
\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}}.
\end{eqnarray}

There is also an additional solution for Eqs.~(\ref{PCG-GPE-a}) and (\ref{PCG-GPE-b}), where $n_b=0$ and
\begin{equation}
\label{P-C-ground-mu-1-n=0}
\mu = c_0 +\Delta +g_{a}n_{a}.
\end{equation}
We note that solution (\ref{P-C-ground-nb}) exists as long as $n_b>0$. Otherwise, only the ground state with $n_b=0$ and $\mu$ given in Eq.~(\ref{P-C-ground-mu-1-n=0})
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\neq0$ and $n_b\neq0$ 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~\cite{Fradkin-Mackenzie:2010} and in Bose nematics (see, e.g., Ref.~\cite{Jian-Zhai:2011}).

\subsection{Free energy}
\label{sec:phenomenology-collective-F}

In order to clarify which of the two possible solutions for the steady-state homogenous Gross-Pitaevskii equations, discussed in Sec.~\ref{sec:phenomenology-collective-ground}, is the true ground state, we calculate the free energy of the system. By using Eq.~(\ref{phen-F-def}), we derive the following expression:
\begin{eqnarray}
\label{P-C-F-nbneq0}
F_{n_b\neq0} &=& \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_an_a\Delta -4 \Delta^2\right\} \nonumber\\
&-&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\},
\end{eqnarray}
where $n_{b}$ given in Eq.~(\ref{P-C-ground-nb}) was used. In the case $n_b=0$, we have a simpler expression
\begin{eqnarray}
\label{P-C-F-n=0}
F_{n_b=0} =-3g_{a} n_a^2.
\end{eqnarray}
In addition to comparing the free energies, one should check whether $n_b$ remains positive for the state with $n_b\neq0$.

We present the phase diagram of the system for a few combinations of parameters in Fig.~\ref{fig:P-C-diff-phase-diagram}. Unless otherwise stated, we use the values of parameters given in Eqs.~(\ref{phen-three-t-num-vars-be}) and (\ref{phen-three-t-num-vars-ee}). 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\neq0$ and $n_b=0$ are possible. In general, however, the ground state with $n_b\neq0$ 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.~\ref{fig:P-C-diff-phase-diagram}(c) and \ref{fig:P-C-diff-phase-diagram}(d)].

\begin{figure*}[!ht]
\centering
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig7a.eps}}
\hspace{0.1\textwidth}
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig7b.eps}}\\
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig7c.eps}}
\hspace{0.1\textwidth}
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig7d.eps}}
\caption{The phase diagram of the system for a few combinations of parameters: $g_{a}$ and $g_{b}$ (panel (a)), $g_a$ and $g_{ab}$ (panel (b)), $\Delta$ and $g_{a}$ (panel (c)), as well as $\Delta$ and $n_{a}$ (panel (d)). Red and blue colors correspond to the ground states with $n_{b}\neq0$ and $n_b=0$, respectively. By default, the values of parameters given in Eqs.~(\ref{phen-three-t-num-vars-be}) and (\ref{phen-three-t-num-vars-ee}) are used.
}
\label{fig:P-C-diff-phase-diagram}
\end{figure*}

\subsection{Dispersion relations}
\label{sec:phenomenology-collective-lin}

In this section, we use the solutions found in Sec.~\ref{sec:phenomenology-collective-ground} and linearize the Gross-Pitaevskii equations (\ref{phen-GPE-a-def}) and (\ref{phen-GPE-b-def}). By using relations (\ref{phen-CM-a-uv}) and (\ref{phen-CM-b-uv}), one obtains the following equations for $u_{a,b}$ and $v_{a,b}$:
\begin{eqnarray}
\label{PCG-GPE-ua}
\omega u_a &=& \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\nonumber\\
&+& u_{ab} \left[2\psi_{b,0}\psi_{a,0}^{*} u_b +(\psi_{b,0})^2 v_a\right],\\
\label{PCG-GPE-ub}
\omega u_b &=& \left(c_0 -\Delta -\mu\right)u_b +v (k_x+ik_y) u_a +2 g_{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\nonumber\\
&+& u_{ab}^{*} \left[2\psi_{a,0}\psi_{b,0}^{*} u_a +(\psi_{a,0})^2 v_b\right],\\
\label{PCG-GPE-va}
-\omega v_a &=& \left(c_0 +\Delta -\mu\right)v_a -v (k_x+ik_y) v_b +2 g_{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\nonumber\\
&+& u_{ab}^{*} \left[2\psi_{b,0}^{*}\psi_{a,0} v_b +(\psi_{b,0}^{*})^2 u_a\right],\\
\label{PCG-GPE-vb}
-\omega v_b &=& \left(c_0 -\Delta -\mu\right)v_b -v (k_x-ik_y) v_a +2 g_{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\nonumber\\
&+& u_{ab} \left[2\psi_{a,0}^{*}\psi_{b,0} v_a +(\psi_{a,0}^{*})^2 u_b\right].
\end{eqnarray}
As before, we omit the valley index $\zeta$ assuming that the valleys have the same dynamics. Note that complex conjugation was performed in Eqs.~(\ref{PCG-GPE-va}) and (\ref{PCG-GPE-vb}).

Equating the determinant of the system (\ref{PCG-GPE-ua}) through (\ref{PCG-GPE-vb}) 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}\neq0$, see Appendix~\ref{sec:App-Josephson-omega}.

In the case of the nontrivial ground state solution given in Eq.~(\ref{P-C-ground-nb}), the frequencies are
\begin{equation}
\label{P-C-lin-omega-u3a=u3b=0}
\omega_{\pm}=\sqrt{v^2k^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]}}.
\end{equation}
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 (\ref{P-C-lin-omega-u3a=u3b=0}) 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~\cite{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
\begin{eqnarray}
\label{P-C-lin-omega-nb=0}
\omega^2_{\pm} &=& v^2k^2 +\frac{\left[2\Delta +n_a(g_{a}-g_{ab})\right]^2}{2} \nonumber\\
&\pm& \frac{1}{2}\sqrt{\left\{2v^2k^2 +\left[2\Delta +n_a(g_{a}-g_{ab})\right]^2\right\}^2 -4\left\{v^4k^4 +v^2k^2g_{a}n_{a} \left[2\Delta +(g_{a}-g_{ab})n_a\right]\right\}}.
\end{eqnarray}
The gap between $\omega_{+}$ and $\omega_{-}$ branches reads as
\begin{equation}
\label{P-C-lin-omega-gap-nb=0}
\left|\omega_{+}(k=0)-\omega_{-}(k=0)\right| = \left|2\Delta +(g_{b}-g_{ab})n_{a}\right|.
\end{equation}
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 (\ref{P-C-lin-omega-u3a=u3b=0}) and (\ref{P-C-lin-omega-nb=0}) in Figs.~\ref{fig:P-C-lin-sol-1-Delta-2D} and \ref{fig:P-C-lin-sol-2-Delta-2D}, respectively. It is notable that the rotational symmetry is spontaneously broken in the interacting BEC of magnons at $n_b\neq0$. Indeed, this is evident from the expression in Eq.~(\ref{P-C-lin-omega-u3a=u3b=0}), 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\neq0$ is not always the ground state of the system (see the discussion in Sec.~\ref{sec:phenomenology-collective-F}). 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\neq0$ 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\neq0$. 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$.

\begin{figure*}[!ht]
\centering
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig8a.eps}}
\hspace{0.1\textwidth}
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig8b.eps}}
\caption{Frequencies of the collective modes $\omega_{+}$ (panel (a)) and $\omega_{-}$ (panel (b)) for $n_b\neq0$ at a few values of the Haldane gap $\Delta$. Here, $k_0=1/(vt_0)$ and $\omega_0=vk_0$. The dispersion relation for $\mathbf{k}\parallel\hat{\mathbf{x}}$ is trivial, i.e., $\omega_{\pm}=vk$. Shaded regions denote the imaginary part of the frequencies. Dispersion relations are gapless for all cases under consideration.}
\label{fig:P-C-lin-sol-1-Delta-2D}
\end{figure*}

\begin{figure*}[!ht]
\centering
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig9a.eps}}
\hspace{0.1\textwidth}
\subfigure[]{\includegraphics[height=0.35\textwidth]{Fig9b.eps}}
\caption{Frequencies of the collective modes $\omega_{+}$ (panel (a)) and $\omega_{-}$ (panel (b)) for $n_b=0$ at a few values of the Haldane gap $\Delta$. Here, $k_0=1/(vt_0)$, $\omega_0=vk_0$, and $\Delta_{\rm cr}=|(g_{b}-g_{ab})n_{a}|/2$. The dependence on $k_y$ is the same as for $k_x$. Shaded regions denote the imaginary part of the frequencies. Except a single point $\Delta =\Delta_{\rm cr}$, one of the modes is gapped and the other is gapless.}
\label{fig:P-C-lin-sol-2-Delta-2D}
\end{figure*}

Let us discuss the dependence of the frequencies of the collective modes on the Haldane gap parameter $\Delta$. As expected from the analytical result (\ref{P-C-lin-omega-u3a=u3b=0}), 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.~\ref{fig:P-C-lin-sol-1-Delta-2D}, one of the modes ($\omega_{+}$) becomes unstable.
The Haldane gap can be also used to stabilize collective modes for $n_b=0$ (see Fig.~\ref{fig:P-C-lin-sol-2-Delta-2D}). 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.

\section{Summary and outlook}
\label{sec:Summary}

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~\cite{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~\cite{Kruglyak-Grundler:2010} and nitrogen-vacancy centers~\cite{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~\cite{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~\cite{Nowik-Boltyk-Demokritov:2012}, as well as time-resolved magneto-optic Kerr effect measurements~\cite{Kruglyak-Grundler:2010} and nitrogen-vacancy centers~\cite{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.

\begin{acknowledgments}
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.
\end{acknowledgments}


\section*{Acknowledgment}

\bibliographystyle{IEEEtran}

\section{Abstract}
New advances in technology have greatly improved the monitoring and controlling of power networks, but these advances leave the system open to cyber attacks. One common attack is known as a False Data Injection Attacks (FDIAs), which poses serious threats to the operation and control of power grids. Hence, recent literature has proposed various detection and identification methods for FDIAs, but few studies have focused on a solution that would prevent such attacks from occurring. However, great strides have been made using deep learning to detect attacks. Inspired by these advancements, we have developed a new methodology for not only identifying AC FDIAs but, more importantly, for correction as well. Our methodology utilizes a Long-Short Term Memory Denoising Autoencoder (LSTM-DAE) to correct attacked-estimated states based on the attacked measurements. The method was evaluated using the IEEE 30 system, and the experiments demonstrated that the proposed method was successfully able to identify the corrupted states and correct them with high accuracy.

\section{Brief summary of State Estimation and FDIAs}

\subsection{State Estimation}
State estimation is a power system application used to estimate the state of power transmission networks through a set of sensor measurements and network topology information \cite{Bobba2010}. Network operators use the estimations of the current state to ensure that the power grid is running in the desired states, take corrective control actions, if needed, and plan for any emergencies \cite{Bobba2010} \cite{Qingyu2014}.

There are two types of state estimation based on the power flow models of the network, namely DC and AC state estimations. A DC state estimation is a linear state estimation that is considered as a simplified version of the more complex AC or nonlinear state estimation. The major differences between DC and AC state estimations are: 1) in the DC state estimation, the solution is obtained in a closed form, while in the AC state estimation, it is obtained through iteration; 2) DC state estimation is based on active power flow analysis while the AC state estimation is based on reactive power flow analysis; and 3) in the DC state estimation, the only variables are the voltage phase angles, however, in the AC state estimation both the voltage phase angles and magnitudes are system variables \cite{Rahman2013}.

For  DC state estimation, the relationship between the network's measurement and its state can be expressed by the following equation \cite{Bobba2010}:

\begin{equation}
z=Hx + e
\end{equation}

where \(z=( z_1 , z_2, z_3, ...,z_m)^T\)is the sensor measurements, \( x = (x_1,x_2,...,x_n)^T\) is the true states of the system, ${H}$ is the Jacobian matrix where ${Hx}$ is a vector of m linear functions linking measurement to states, and \(e = (e_1,e_2,...,e_m)^T\) is the random errors in the measurement.

Since the DC model is a linear estimation of the network state, then the estimated system state ${x_{est}}$ can be expressed as:

\begin{equation}
{x_{est}}  =  (H^T WH)^{-1}H^TWz
\end{equation}

Where W is a diagonal matrix that consists of the measurement weights.
As for the AC state estimation, since the power flows are nonlinearly dependent, the relationship between the network’s measurement and its state can be expressed by the following equation \cite{Gabriela2012}:
\begin{equation}
z=h(x) + e
\end{equation}

where ${h(x)}$ is a function vector that establishes the relationships between measured values and state variables.
Consequently, the estimated state variables are determined from the following optimization function:
\begin{equation}
min F(x) = (z - h(x))^T.W.(z-h(x))
\end{equation}

\subsection{FDIAs}
False data injection attacks can be defined as those in which an attacker changes the sensor's measurements, which in turn changes the estimated value of the state variables without being flagged by the state estimator's bad measurement detection algorithms. FDIAs can be either random or targeted \cite{Bobba2010}. In random attacks, the attacker injects arbitrary errors into the estimates of state variables while in targeted attacks, specific errors are injected into the estimates of specific state variables.  A more significant categorization, however, is based on the extent of the knowledge that the attacker has about the network prior to launching the attack. As such, FDIAs can be either complete or incomplete.

As the name suggests, with complete FDIAs, the attacker has complete knowledge of the power network's topology and transmission-line admittance values, allowing him or her to design a false data injection attack vector \cite{Rahman2012}. This type of attack is rare since the attacker is operating with limited resources and restricted physical access to the power grid. Incomplete attacks, the most common type, occurs when the attacker does not possess complete knowledge about the network's parameters and design \cite{Xuan2015}. For more information on complete and incomplete attacks see\cite{James2018}.

In this research, we deal with complete attacks, assuming that the attacker has the full knowledge required to carry out the attack. The attacker’s full knowledge makes detection of the attack more difficult.

\section{Introduction}\label{sec:introduction}}
\else
\section{Introduction}
\label{sec:introduction}
\fi

\IEEEPARstart{T}{his} demo file is intended to serve as a ``starter file''
for IEEE Computer Society journal papers produced under \LaTeX\ using
IEEEtran.cls version 1.8b and later.
I wish you the best of success.

\hfill mds

\hfill August 26, 2015

\subsection{Subsection Heading Here}
Subsection text here.

\subsubsection{Subsubsection Heading Here}
Subsubsection text here.

\section{Conclusion}
The conclusion goes here.

\appendices
\section{Proof of the First Zonklar Equation}
Appendix one text goes here.

\section{}
Appendix two text goes here.

\ifCLASSOPTIONcompsoc

  \section*{Acknowledgments}
\else

  \section*{Acknowledgment}
\fi

The authors would like to thank...

\ifCLASSOPTIONcaptionsoff
  \newpage
\fi

\section{Case Study}
In this research, we validated our results on the IEEE-30 bus system using real-world load data found on the Independent Electricity System Operator (IESO) website\cite{data}. Specifically, we trained our system on the data from 2002 to 2007 and tested it on the data from 2013 to 2015 to ensure that our system is not load dependent.

For FDIAs construction, we implemented the same methodology described in \cite{James2018}. To train and evaluate our system, ${w}$ timesteps samples were sequentially sampled from the normal dataset, and the last sample in the window has been attacked such that it is not detected by the state estimator. Thus, ${w-1}$ samples are normal ${x}$, and the last sample is the attacked state vector $\hat{x}$. Using this setup, 200,000 of ${w}$ samples have been constructed. We trained, evaluated, and tested our system with a ratio of 60/20/20, respectively.

For our experiments, we found that the best window size setting is 5. For the evaluation metric, we used the Root Mean Square Error (RMSE) to calculate the difference between the normal and the attacked values \cite{Chai2014}. Other metrics were also tested, such as Mean Absolute Error (MAE) and Mean Square Error (MSE), but RMSE gave the most easily interpretable results.

All the results and simulations were conducted on Nvidia K-80 GPU along with 16 GB of RAM. The LSTM-DAE was constructed using Keras, which is a top-level library built on top of TensorFlow for a computational speed boost.

 \begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth,keepaspectratio]{figure/ModelLoss.png}
    \caption{Model Loss (RMSE)}
    \label{fig:Results2}
\end{figure}
The results showed that our model is capable of capturing the temporal-spatial behavior of the data. As shown in Figure \ref{fig:Results2} , the loss of both the training and validation sets was monitored while training the model, and the results revealed that the RMSE was gradually decreasing until it reached less than 0.0013 on the validation data after 15 epochs. The RMSE reported on the testing set was 0.001 based on the best validation model.
Then, to validate that the model was capable of reconstructing any attacked state, a random sample was taken from the test set and the normal, attacked, and the corrected states’ values were measured. This is visualized in Figure \ref{fig:Reconstruction}, as the attacked state number 42 was successfully reconstructed with high accuracy.
Moreover, the RMSE values were calculated for all the 60 states to ensure once again that the model was capable of reconstructing any attacked state. The results in Figure \ref{fig:RMSE_60} show that this was achieved as the maximum RMSE was less than 0.0025, which is very low.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth,keepaspectratio]{figure/Reconstruction.png}
    \caption{Reconstruction of a random sample for 5 random states}
    \label{fig:Reconstruction}
\end{figure}

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth,keepaspectratio]{figure/RMSE_60.png}
    \caption{RMSE values for all states}
    \label{fig:RMSE_60}
\end{figure}
Finally, the difference between the actual states and the corrected one was calculated for all the states and timesteps in the test dataset. It was found that the maximum difference is less than 0.004, and most of the differences are between 0 and 0.002. In order to visualize this, a histogram graph was plotted in Figure \ref{fig:Histogram}.

Surprisingly, our model is capable of almost reconstructing the attacked states to their normal values with a very low reconstruction error rate. Based on the reconstruction high accuracy, we were able to identify the source of the attack with a simple threshold mechanism,  which gave us 100 percent identification accuracy. We believe that our model has enough capacity to capture the behavior of the system in its steady-state given a real-world load on unseen load data. These very promising results have encouraged us to pursue research to model the state estimation itself using our approach.

\begin{figure}[t]
    \centering
    \includegraphics[width=\columnwidth,keepaspectratio]{figure/Histogram.png}
    \caption{Histogram of the difference between actual and corrected statse.}
    \label{fig:Histogram}
\end{figure}

\section{Conclusion}

In this paper, we proposed a new FDIAs identification and correction for AC state estimation. While some papers tackled the detection of FDIAs, little work has been done on identifying and correcting these attacks. The proposed method consists of an LSTM-DAE, which can capture the temporal-spatial dependencies. We validated our model on the IEEE 30-bus system using real-world load data. The results showed the model was able to identify the attacked states and correct them with high accuracy.
\section{Introduction}

Power systems and power grid networks are critical components in our modern-day electrical infrastructure \cite{Kush20018}. As such, ensuring both the physical and cyber security of these systems is of utmost importance. The use of information technology in power systems offers various benefits such as improved monitoring and controlling capability as well as allowed for the integration of renewable energy sources. But with these advantages come new risks. In recent years, cyber threats have been recognized as the most significant and dangerous threats to our power systems as a result of the increased dependency on information technology \cite{Gaoqi2013}.

State estimation is an important power system application used to estimate the state of power transmission networks through a set of sensor measurements and network topology information \cite{Bobba2010}. The state estimator collects data about sensor measurements from the Supervisory Control and Data Acquisition (SCADA) system, including information about bus voltages, active and reactive power injection at each bus, and power flow in branches \cite{Ayad2018}.  It also removes all noise and error \cite{Kush20018}, as well as estimates the state of the power network\cite{Bobba2010}. These estimations are then used by network operators to ensure that the power grid is running in the desired states, allowing them to take corrective control actions if needed, and plan for any emergencies \cite{Bobba2010} \cite{Qingyu2014}.

One of the most serious types of attacks on state estimation is a new one: False Data Injection Attacks (FDIAs). These attacks target the state estimation of the power grid by tampering with the data that are transmitted through communication lines \cite{Ashrafuzzaman2018} \cite{Nawaz2018}. FDIAs occur when an attacker introduces incorrect data or measurements into the system, which will reduce the control the system's operator has over the system, thus altering the state of the system \cite{Ayad2018}. It poses serious threats to the operations and control of power systems as it can evade the existing bad data detection (BDD) systems \cite{Sagnik2017}. In addition, FDIAs can affect deregulated energy markets and result in severe economic loses \cite{Kush20018}.

To address FDIAs, protection-based strategies, and detection-based strategies for defending the power systems have been proposed in the literature \cite{Ayad2018}. Multiple detection-based methods have been proposed. The most advanced methods use deep learning approaches such as Convolutional Deep Belief Network (CDBN) in \cite{Youbiao2017}, Recurrent Neural Network (RNN) in \cite{Ayad2018}, Deep Neural Network (DNN) in \cite{Ashrafuzzaman2018}, and Convolutional Neural Network (CNN) and Long-Short Term Memory (LSTM) in \cite{Xiangyu2018}\cite{Sagnik2017} to address the problem.  While substantial progress has been made in FDIAs detection, the complexity of the problem continues to present challenges to those in the power systems community. \cite{Kush20018} uses a statistical approach to address FDIAs, correcting measurements before feeding them to the state estimator.  One limitation of this method is the challenge of being scalable as it involves a complicated threshold scheme that would need to be applied to each scenario to detect the attack. Nevertheless, the scarcity of studies that have attempted to go a step further beyond attack detection and correct the errors resulting from FDIAs makes this field a promising one for future research \cite{Yanpeng2018}.

Inspired by the significant success of deep learning (DL) methods in detecting attacks, in this paper, we present a new methodology for AC FDIAs identification and correction using DL to both identify and correct attacked measurements. With this method, we introduce the LSTM Denoising Autoencoder (LSTM-DAE) to reconstruct the signal such that it detects the anomaly and can then return it to its normal state. The new correction method was evaluated using the IEEE 30  system, and the experiments showed that the proposed method was successfully able to identify the corrupted states and correct them with high accuracy.

The remainder of this paper is organized as follows: Section II provides an overview of the state estimation process used in control centers and the false data injection attacks. Section III introduces the methodology of attack identification and correction using the LSTM-DAE. Section IV presents the case study, while Section V shows the simulation results of the case study with a brief discussion of the efficiency of the method used. Conclusions and suggestions for future work are presented in Section VI.

\section{The Proposed LSTM-DAE Method}

 \begin{figure*}[t!]
    \centering
    \includegraphics[width=\linewidth,keepaspectratio]{figure/LSTMNew.jpeg}
    \caption{LSTM Denoising Autoencoder Architecture.}
    \label{fig:ARCH2}
\end{figure*}

 \begin{figure}[t]
    \centering

    \includegraphics[width=80mm,keepaspectratio]{figure/PropMM.jpeg}
    \caption{The Proposed Identification and Correction scheme for FDIAs.}
    \label{fig:ProposedM}
\end{figure}

Autoencoder, in general, is a neural network that is trained in order to learn a high-level representation of data by reducing data dimensionality. This is accomplished by having two main parts which are an encoder part \textbf{h = f(x)} and a decoder part \textbf{r = g(h)}. While vanilla Autoencoder is only useful to compress and decompress data without losing useful information \cite{Pierre2012}, Denoising Autoencoders (DAE) \cite{Vincent2008}, which is a special type of Autoencoder, can also learn to remove the noise from the corrupted input by recovering the original undistorted input. The DAE showed better representation more than the other Autoencoder types as the network learns a higher level of representation that is relatively robust to noise.

In this research, we aim to correct data affected by FDIAs. Inspired by the DAE and due to the sequential nature of the data, we are proposing a special type of RNN Autoencoder, namely, LSTM-DAE, to correct data after FDIA.

To do so,  instead of feeding the original input ${x}$, we feed a noisy corrupted version of it $\hat{x}$, which is, in our case, the attacked states. Then, the loss function in that case is \textbf{L(x,g(f($\hat{x}$)))} as the reconstructed version of the original input ${x'}$ is computed from the corrupted input $\hat{x}$.

LSTM, on the other hand, is a modified version of vanilla RNN, called gated RNN, which takes advantage of the sequential information to estimate the output \cite{Alex2018}, while having the ability to preserve long-term memory relations \cite{Hochreiter1997}. This is done by having three different types of gates, which are forget gate, input gate, and output gate. When they are combined together, the LSTM can be trained without worrying about losing information from the long-term memory or not capturing all the available information from the short-term memory.
Adding all this together, the proposed model consists of DAE stacking multiple layers of LSTM, where the input is a corrupted version of the data (Attack), and the output is the correct version of it (Normal).

The architecture is composed of three main parts, which are the encoder and the decoder, with a bridge between them. The first part consists of three parts, which are the input sample along with two LSTM layers. The input sample is a matrix whose dimensions are the timestep window \textit{w} and the number of states. The first LSTM layer of the encoder part consists of 128 units, and the Rectified Linear Unit (ReLU) is used as the activation function because it is the state-of-the-art activation function that can capture complex relations between the input and the output \cite{Vinod2010}. The second LSTM layer of the encoder part consists of 64 units with also a ReLU activation function.

Following that, the bridge between the encoder and the decoder is represented by a Repeat Vector layer to reshape the output of the second LSTM layer in the encoder part so it can be fed into the first LSTM layer in the decoder part. This is achieved by repeating the output of the second LSTM layer in the encoder a number of times equal to the window size, which is 5 in this case.

For the decoder part, it starts with an LSTM layer with 128 units and a ReLU function as the activation function. This is followed by a Dense layer, which is a normal Neural Network, but it is wrapped in a Time Distributed layer to apply the Dense layer to every temporal slice of the input of the layer. Finally, the output of the Dense layer is the corrected states after training. The flow of the architecture can be seen in Figure \ref{fig:ARCH2}.

Our system works by correcting the samples online, such that the corrected sample is fed back to the states' queue in order to correct the next attacked state. Thus, the step window \textit{w} which represents the input of the model is divided into normal samples (\(x_{t-n},...,x_{t-3},x_{t-2},x_{t-1}\)) and an attacked sample ($\hat{x_t}$). To ensure that the model is immune to the noise coming from the corrected sample, Additive White Gaussian Noise (AWGN) is added to the states in the normal samples while training. This will guarantee that our system can handle the noise generated from the correction process.

Our system assumes a perfect FDIAs detection; thus, our proposed model focuses on identification and correction. Our system identifies the attacked states by processing the corrected states from the proposed system. A direct comparison is applied between the corrected and attacked states after applying thresholds for theta and voltage to identify the attacked states. If any state value exceeds these thresholds, this indicates that state as an attacked state. The scheme of the process is illustrated in Figure \ref{fig:ProposedM}.