1. Introduction
Rotating Rayleigh–Bénard convection, where a rotating fluid between two horizontal plates is driven by bottom heating, is an important prototype for geophysics and astrophysics (Bodenschatz, Pesch & Ahlers Reference Bodenschatz, Pesch and Ahlers2000; Ahlers, Grossmann & Lohse Reference Ahlers, Grossmann and Lohse2009; King, Stellmach & Aurnou Reference King, Stellmach and Aurnou2012; Guervilly, Hughes & Jones Reference Guervilly, Hughes and Jones2014). It occurs in various natural systems, such as the Earth's atmosphere and oceans, as well as in the interiors of stars and exoplanets.
In this convection system, the rotation leads to the emergence of edge states at the sidewall (Zhong, Ecke & Steinberg Reference Zhong, Ecke and Steinberg1991, Reference Zhong, Ecke and Steinberg1993; Ecke, Zhong & Knobloch Reference Ecke, Zhong and Knobloch1992; Goldstein et al. Reference Goldstein, Knobloch, Mercader and Net1993; Ning & Ecke Reference Ning and Ecke1993), which are characterized by their unidirectional propagation and their ability to persist even in the presence of different types of barriers and turbulence (Favier & Knobloch Reference Favier and Knobloch2020; Ecke Reference Ecke2023; Ecke & Shishkina Reference Ecke and Shishkina2023). These states were first discovered through experiments indirectly by Rossby (Reference Rossby1969), where a surprising convection occurs with a Rayleigh number (
$Ra$) below the critical value for the laterally unbounded system (Chandrasekhar Reference Chandrasekhar1961). Zhang & Liao (Reference Zhang and Liao2009) give the asymptotic expression of the critical 
$Ra$ for the onset of these sidewall modes, regaining the leading terms previously calculated by Herrmann & Busse (Reference Herrmann and Busse1993). With increasing 
$Ra$, more sidewall modes emerge, and they may undergo modulational instabilities and interact with bulk modes, resulting in a complex nonlinear dynamics (Zhong et al. Reference Zhong, Ecke and Steinberg1993; Liu & Ecke Reference Liu and Ecke1999; Horn & Schmid Reference Horn and Schmid2017). In the turbulent flows with high 
$Ra$, the sidewall modes may be transformed to boundary zonal flows, which can also coexist with bulk convection (de Wit et al. Reference de Wit, Aguirre Guzmán, Madonia, Cheng, Clercx and Kunnen2020; Zhang et al. Reference Zhang, van Gils, Horn, Wedi, Zwirner, Ahlers, Ecke, Weiss, Bodenschatz and Shishkina2020; Ecke, Zhang & Shishkina Reference Ecke, Zhang and Shishkina2022). Beyond these detailed studies, a new insight is proposed that these edge states may be related to the topological nature of the system (Favier & Knobloch Reference Favier and Knobloch2020). Thus, a new set of questions arises, such as, can we obtain the topological invariant of the system? How robust in quantitative terms is the topology in the presence of turbulence? The answers to these follow-up questions will deepen our understanding of the overall flow structure and the spatial distribution of heat flux in geophysical and astrophysical contexts (Zhong et al. Reference Zhong, Stevens, Clercx, Verzicco, Lohse and Ahlers2009; Terrien, Favier & Knobloch Reference Terrien, Favier and Knobloch2023).
Edge states under topological protection have been a topic of great interest in recent years. These edge states are robust against disorder and perturbations, and their topological protection originates from the non-trivial topology of the underlying bulk system. This topological nature can be indicated by a corresponding topological invariant, such as the Chern number (Xiao, Chang & Niu Reference Xiao, Chang and Niu2010; Delplace, Marston & Venaille Reference Delplace, Marston and Venaille2017) or the winding number (Zhu, Li & Marston Reference Zhu, Li and Marston2023). When the topological invariant is non-zero, the bulk-boundary correspondence guarantees the existence of robust edge states (Essin & Gurarie Reference Essin and Gurarie2011; Venaille et al. Reference Venaille, Onuki, Perez and Leclerc2023; Onuki, Venaille & Delplace Reference Onuki, Venaille and Delplace2024). The topological physics first appeared in condensed matter (Kosterlitz & Thouless Reference Kosterlitz and Thouless1973; Thouless et al. Reference Thouless, Kohmoto, Nightingale and den Nijs1982; Haldane Reference Haldane1988), such as topological insulators (Hasan & Kane Reference Hasan and Kane2010), topological photonics (Ozawa et al. Reference Ozawa2019) and topological phononics (Liu, Chen & Xu Reference Liu, Chen and Xu2020). Recently, similar topological properties have been found in macroscopic systems, especially in hydrodynamics, such as the topological origin of trapped waves near the equator or coastlines (Delplace et al. Reference Delplace, Marston and Venaille2017; Venaille & Delplace Reference Venaille and Delplace2021), the topological waves in fluids with odd viscosity (Souslov et al. Reference Souslov, Dasbiswas, Fruchart, Vaikuntanathan and Vitelli2019) and the presence of topological invariants in active matter systems (Shankar et al. Reference Shankar, Souslov, Bowick, Marchetti and Vitelli2022). This universality of topological states can be attributed to the similar symmetries embodied in systems of different scales (Senthil Reference Senthil2015).
In this paper, we point out that linearized rotating convection is a non-Hermitian system and that these edge states are manifestations of the non-trivial topological Berry phase in the bulk. Non-Hermitian systems are open systems that do not obey the Hermitian symmetry property (Ashida, Gong & Ueda Reference Ashida, Gong and Ueda2020). Thus, their eigenvalues are not necessarily real, and the corresponding eigenvectors can be unorthogonal. Rayleigh–Bénard convection is non-Hermitian because it exchanges energy with an external heat source. Non-Hermitian systems exhibit unconventional physical properties, including non-reciprocal transmission, exceptional points and topological phase transitions (Bergholtz, Budich & Kunst Reference Bergholtz, Budich and Kunst2021; Ding, Fang & Ma Reference Ding, Fang and Ma2022). These phenomena have attracted significant attention in recent years due to their potential applications in various fields of physics, such as optics, condensed matter physics and quantum information science.
With the non-Hermitian and topological nature, we foresee that rotating convection will become a more abundant system and act as a platform for probing topological states in turbulent flows. This paper is structured as follows. Section 2.1 introduces the non-Hermitian Hamiltonian matrix of the linearized governing equations. Section 2.2 explains the calculation of Chern number in non-Hermitian systems. Section 3 focuses on calculating the topological invariants in the unstably stratified case. Section 4 explores the topological properties of the system in stable and critical cases. Section 5 discusses the calculation of the Chern number under the hydrostatic approximation. Section 6 visualizes the topological nature through winding numbers. Finally, discussions and summaries are presented in § 7. Appendix A discusses the solutions satisfying realistic boundary conditions.
2. The non-Hermitian eigenvalue problem of linearized rotating convection
2.1. Non-Hermitian Hamiltonian matrix
The governing equations for rotating Rayleigh–Bénard convection include the Navier–Stokes, continuity and temperature equations. The dimensionless form can be written as (cf. Favier & Knobloch Reference Favier and Knobloch2020)
$$\begin{gather} \frac{\partial\boldsymbol{u}}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\boldsymbol{u} ={-}\boldsymbol{\nabla}p-\lambda\boldsymbol{e}_{z}\times\boldsymbol{u} +\alpha\theta\boldsymbol{e}_{z}+E\nabla^{2}\boldsymbol{u}, \end{gather}$$
$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$
$$\begin{gather}\frac{\partial\theta}{\partial t}+\left(\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\right)\theta=w+\frac{E}{Pr}\nabla^{2}\theta, \end{gather}$$
where 
$\boldsymbol {u}=(u,v,w)$ is the velocity vector, 
$\lambda =\pm 1$ indicates the direction of rotation, 
$p$ is the pressure, 
$\theta$ is the temperature fluctuation relative to the linear conduction profile, 
$\alpha =RaE^{2}/Pr$ is the square of the convective Rossby number. The Rayleigh number is 
$Ra=\beta g\Delta Th^{3}/\nu \kappa$, where 
$\beta$ is the thermal expansion coefficient, 
$g$ is the acceleration due to gravity, 
$\Delta T$ is the temperature difference between the top and bottom plates, 
$h$ is the height of the fluid layer and 
$\nu$ and 
$\kappa$ are the viscosity and thermal diffusivity, respectively. It is a measure of the strength of the buoyancy-driven flow relative to the viscous forces. The Ekman number, 
$E =\nu /(2\varOmega h^2)$, describes the balance of viscous forces to Coriolis forces, where 
$\varOmega$ is the angular velocity of rotation. The Prandtl number is 
$Pr=\nu /\kappa$ and is a measure of the relative importance of viscous and thermal diffusion. It is worth noting that 
$\alpha$ is independent of both 
$\nu$ and 
$\kappa$, so that we can take the inviscid limit 
$E \to 0$ when necessary at a fixed non-zero 
$\alpha$.
Considering a normal-mode ansatz 
$(u,v,w,p,\theta ) = {\rm Re}\{(\hat {u},\hat {v},\hat {w},\hat {p},\hat {\theta }) \exp ({{\rm i}(\boldsymbol {k}\boldsymbol {\cdot } \boldsymbol {r}\,{-}\,\omega t)})\}$ with the wavenumber 
$\boldsymbol {k}=(k_x,k_y,k_z)$, the linearized equations become
\begin{gather} -{\rm i}\omega\boldsymbol{\hat{u}}={-}{\rm i}\hat{p}\boldsymbol{k}-\lambda\boldsymbol{e}_{z}\times\boldsymbol{\hat{u}}+\alpha\hat{\theta}\boldsymbol{e}_z-Ek^{2}\boldsymbol{\hat{u}},\end{gather}
\begin{gather} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\hat{u}}=0, \end{gather}
\begin{gather} -{\rm i}\omega\hat{\theta}=\hat{w}-\frac{E}{Pr}k^{2}\hat{\theta}, \end{gather}
where 
$k=\sqrt {k_{x}^{2}+k_{y}^{2}+k_{z}^{2}}$. Taking the divergence of both sides of (2.1a) and combining with (2.1b), we get
\begin{equation} 0=\hat{p} k^{2}+{\rm i}\lambda(k_{x}\hat{v}-k_{y}\hat{u})+{\rm i}\alpha k_{z}\hat{\theta}. \end{equation}Then
\begin{equation} \hat{p}={\rm i}\frac{\lambda({-}k_{x}\hat{v}+k_{y}\hat{u})-\alpha k_{z}\hat{\theta}}{k^{2}}. \end{equation} With the wave vector 
$\boldsymbol{\psi} \equiv (\hat {u},\hat {v},\hat {w},\hat {\theta })$, we get the eigenequation from (2.2) and (2.4) as
\begin{equation} \boldsymbol{\mathsf{H}}\boldsymbol{\psi}=\omega\boldsymbol{\psi}, \end{equation}with
\begin{equation} \boldsymbol{\mathsf{H}}={\rm
i}\left[\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}}
E_{0}+\lambda\dfrac{k_{x}k_{y}}{k^{2}} &
\lambda\left(1-\dfrac{k_{x}^{2}}{k^{2}}\right) & 0 &
-\alpha\dfrac{k_{x}k_{z}}{k^{2}}\\
\lambda\left({-}1+\dfrac{k_{y}^{2}}{k^{2}}\right) &
E_{0}-\lambda\dfrac{k_{x}k_{y}}{k^{2}} & 0 &
-\alpha\dfrac{k_{y}k_{z}}{k^{2}}\\
\lambda\dfrac{k_{y}k_{z}}{k^{2}} &
-\lambda\dfrac{k_{x}k_{z}}{k^{2}} & E_{0} &
\alpha-\alpha\dfrac{k_{z}^{2}}{k^{2}}\\ 0 & 0 & 1 & E_{1}
\end{array}\right], \end{equation}
where 
$E_{0}=-Ek^{2}$, 
$E_{1}=E_{0}/Pr$. This is a non-Hermitian Hamiltonian because 
$H\neq H^{{\dagger} }$, where 
${\dagger}$ denotes the conjugate transpose. Because we define the growth rate as 
${\rm i} \omega$ for convenience, the eigenvalue 
$\omega$ can be a complex number, unlike in the Hermitian case where 
$\omega$ is real. The eigenvalue problem of 
$H^{{\dagger} }$ is
\begin{equation} \boldsymbol{\mathsf{H}}^{{{\dagger}}}\boldsymbol{\psi}^{\prime}=\omega^{{\ast}}\boldsymbol{\psi}^{\prime}, \end{equation}
where 
$\omega ^{\ast }$ is the complex conjugate of 
$\omega$, and 
$\boldsymbol{\psi}^{\prime }$ is generally different from 
$\boldsymbol{\psi}$.
In (2.1), we do not include the boundary conditions, and the main text of this paper focuses on demonstrating the idea that, to calculate a Chern number for rotating convection, a complex vertical wavenumber is required. A discussion on realistic boundary conditions can be found in Appendix A.
2.2. Chern number in a non-Hermitian system
The Chern number is a topological invariant that was originally defined for systems with a periodic structure, such as crystals, by integrating the Berry curvature over the Brillouin zone (Zak Reference Zak1989; Xiao et al. Reference Xiao, Chang and Niu2010). It is a quantity with integer values and is numerically equal to the Berry phase divided by 
$2{\rm \pi}$. The Berry curvature characterizes the local geometry of wave polarization and is recognized to manifest in the equations of motion of wave packets with multiple components (Berry Reference Berry1984; Perez, Delplace & Venaille Reference Perez, Delplace and Venaille2021).
As a topological invariant, the Chern number is insensitive to small perturbations of the Hamiltonian that do not change the topology of the system. A non-zero Chern number implies the existence of a non-trivial bulk topology and the presence of robust edge states (Hasan & Kane Reference Hasan and Kane2010). The bulk-boundary correspondence principle links the topological properties of the bulk and edge states (Essin & Gurarie Reference Essin and Gurarie2011; Venaille et al. Reference Venaille, Onuki, Perez and Leclerc2023; Onuki et al. Reference Onuki, Venaille and Delplace2024).
In the wavenumber space, the Berry connection (Berry Reference Berry1984) is defined as
\begin{equation} \boldsymbol{A}_{n}(\boldsymbol{k})={\rm i}\langle \boldsymbol{\psi}_{n}(\boldsymbol{k})|\boldsymbol{\nabla}_{\boldsymbol{k}}|\boldsymbol{\psi}_{n}(\boldsymbol{k})\rangle, \end{equation}
where 
$\boldsymbol{\psi} _{n}$ is the eigenvector, and 
$\boldsymbol {\nabla }_{\boldsymbol {k}}$ is the gradient with respect to the wave vector 
$\boldsymbol {k}$. As with the mathematical symbols of quantum mechanics, the right vector 
$|\boldsymbol{\psi} _{n}\rangle$ represents the general eigenvector, the left vector 
$\langle \boldsymbol{\psi} _{n}|$ represents the conjugate transpose of 
$|\boldsymbol{\psi} _{n}\rangle$ and 
$\langle \cdots |\cdots \rangle$ represents the inner product. The Berry connection is a vector-valued function whose integral along a closed path gives the Berry phase. Then, the Berry curvature is calculated by the curl of the Berry connection
\begin{equation} \boldsymbol{\varOmega}_{n}(\boldsymbol{k})=\boldsymbol{\nabla}_{\boldsymbol{k}}\times\boldsymbol{A}_{n}(\boldsymbol{k}). \end{equation}The Chern number is then defined as the integral of the Berry curvature over the Brillouin zone
\begin{equation} C_{n}=\frac{1}{2{\rm \pi}}\int_{BZ}\boldsymbol{\varOmega}_{n}(\boldsymbol{k})\boldsymbol{\cdot} d\boldsymbol{S}, \end{equation}
where 
$d\boldsymbol {S}$ is the surface element of the Brillouin zone. The Berry curvature 
$\boldsymbol {\varOmega }_{n}$ is a pseudovector in three dimensions. Due to the direction of rotation in this system, we only consider its component in the 
$z$ direction, 
$\varOmega _{n}^{z}$, and 
$d\boldsymbol {S}$ takes the horizontal plane.
In non-Hermitian systems, the classic Brillouin zone is no longer sufficient to describe the band topology due to the significant difference in frequency spectra between open (non-Bloch) and periodic (Bloch and classic) boundaries. Therefore, a generalized Brillouin zone (GBZ) defined on the complex wavenumber space is needed (Yokomizo & Murakami Reference Yokomizo and Murakami2019; Yang et al. Reference Yang, Zhang, Fang and Hu2020; Wu et al. Reference Wu, Xie, Zhou and An2022). As a result, a non-Bloch Chern number is defined as the integral of the Berry curvature over the GBZ
\begin{equation} C_{n}=\frac{1}{2{\rm \pi}}\int_{GBZ}\boldsymbol{\varOmega}_{n}(\tilde{\boldsymbol{k}})\boldsymbol{\cdot} d\tilde{\boldsymbol{S}}, \end{equation}
where 
$\tilde {\boldsymbol {k}}$ is complex and 
$d\tilde {\boldsymbol {S}}$ is the surface element of the GBZ. The Berry curvature can take a complex number, as long as the Chern number obtained by integration remains an integer (Yao, Song & Wang Reference Yao, Song and Wang2018).
Due to the non-Hermitian nature of the system (2.7), obtaining the Berry curvature requires the consideration of a biorthogonal basis set, 
$\{|\boldsymbol{\psi} _{n}\rangle \}$ and 
$\{|\boldsymbol{\psi} ^{\prime }_{n}\rangle \}$, which satisfies
$$\begin{gather} \boldsymbol{\mathsf{H}}|\boldsymbol{\psi}_{n}\rangle=\omega_{n}|\boldsymbol{\psi}_{n}\rangle, \end{gather}$$
$$\begin{gather}\boldsymbol{\mathsf{H}}^{{{\dagger}}}|\psi^{\prime}_{m}\rangle=\omega_{m}^{{\ast}}|\psi^{\prime}_{m}\rangle, \end{gather}$$
$$\begin{gather}\langle\psi^{\prime}_{m}|\boldsymbol{\psi}_{n}\rangle=\delta_{m,n}, \end{gather}$$
where 
$\omega _m$ and 
$\omega _n$ are eigenvalues with band indices 
$m$ and 
$n$. Then we can define the non-Hermitian Berry curvature of band 
$\omega _{n}$ as (Yao et al. Reference Yao, Song and Wang2018)
\begin{equation} \varOmega_{n}^{z}={\rm i}\left[\left\langle \frac{\partial\boldsymbol{\psi}_{n}^{\prime}}{\partial k_{x}}\left|\frac{\partial\boldsymbol{\psi}_{n}}{\partial k_{y}}\right.\right\rangle -\left\langle \frac{\partial\boldsymbol{\psi}_{n}^{\prime}}{\partial k_{y}}\left|\frac{\partial\boldsymbol{\psi}_{n}}{\partial k_{x}}\right.\right\rangle \right]. \end{equation}3. Topology of rotating convection
3.1. Complex wavenumber 
$k_{z}$
When the eigenvalues of the linearized equations (2.2)–(2.4) are not real, which often happens in non-Hermitian systems, the convection is linearly unstable with oscillatory modes that can grow or decay over time, and exploring the topological properties gets complicated. For example, in the inviscid case where 
$E=0$, the eigenvalues and eigenvectors are degenerated at 
$k_{x}^{2}+k_{y}^{2}=k_{z}^{2}/\alpha$ (cf. (3.4)), and a larger 
$k_{x}$ or 
$k_{y}$ would make the eigenvalues purely imaginary. Shen, Zhen & Fu (Reference Shen, Zhen and Fu2018) proposed that the non-Hermitian system loses its well-defined Chern numbers unless 
$k_{z}$ becomes complex: with a non-zero imaginary part in 
$k_{z}$, the frequency bands will always remain gapped in the complex space, and the integration of the Berry curvature will give an integer Chern number. After such a procedure, we have obtained a GBZ, using which we can treat non-Hermitian problems under open boundary conditions as if they were under periodic boundary conditions (Yokomizo & Murakami Reference Yokomizo and Murakami2019; Yang et al. Reference Yang, Zhang, Fang and Hu2020; Wu et al. Reference Wu, Xie, Zhou and An2022).
For the rotating convection system, figure 1 illustrates the eigenfrequency gap opened by introducing a complex vertical wavenumber 
$k_z$ with 
$k_z^2=1+0.1{\rm i}$ and 
$\alpha =1$. We can see that both the real and imaginary parts of 
$\omega$ approach 
$0$, but its absolute value is bounded from below. When 
$k_z^2=1+\epsilon {\rm i}$ with 
$\epsilon \ll 1$, one can asymptotically find that
\begin{equation} |\omega|_{gap} = \sqrt{\frac{\epsilon\alpha}{1+\alpha}} + {\rm h.o.t.}, \end{equation}
where 
${\rm h.o.t.}$ denotes high-order terms.
Figure 1. Illustration of the eigenvalue with complex 
$k_z$. Here, 
$k_z^2=1+0.1{\rm i}$, 
$\alpha =1$ and 
$E=0$.
Complex 
$k_z$ is also physically reasonable. The top and bottom no-slip boundaries introduce an exponential decay form to the asymptotic solutions along the 
$z$-direction (Zhang & Liao Reference Zhang and Liao2009), which agrees with our numerical results in Appendix A. The choice of a complex 
$k_z$ is due to the boundary conditions we have adopted: the top and bottom surfaces are no slip and the sidewalls are smooth (where both 
$k_x$ and 
$k_y$ are real). If we consider another scenario where the sidewalls are no slip but the upper and lower surfaces are smooth, 
$k_x$ and 
$k_y$ need a specific imaginary part. Mathematically speaking, 
$k_z$ can take real values in this case, and we can still obtain a well-defined Chern number. On the other hand, under open boundary conditions, i.e. with external energy injection, the bulk eigenstates of non-Hermitian systems exhibit a localized behaviour towards the boundary, the so-called non-Hermitian skin effect (Okuma et al. Reference Okuma, Kawabata, Shiozaki and Sato2020; Zhang, Yang & Fang Reference Zhang, Yang and Fang2022; Zhang et al. Reference Zhang, Zhang, Lu and Chen2022), which differs from the extended Bloch waves in Hermitian systems. In our system, energy exchanges with the outsiders in the 
$z$-direction due to the temperature difference, so we take a complex 
$k_{z}$.
3.2. Inviscid topological invariant ($E=0$)
With the GBZ where 
$k_{z}$ is complex, we can calculate the Chern number for open boundaries through a normal process, just as we do under periodic boundary conditions. The solutions satisfying realistic boundary conditions can be found in Appendix A. To obtain a single-valued function when dealing with the square root, we redefine the square root operation on the complex domain, making its value domain lie in the right half-plane that does not include the negative imaginary axis, i.e. 
$\sqrt {-1}=i$.
We start to discuss the topological properties of the system with the simplest case 
$E=0$. Since 
$k_{z}$ is complex, the singularity in the eigenvalues and eigenvectors will not exist. The Hamiltonian matrix becomes
\begin{equation} \boldsymbol{\mathsf{H}}_{0}={\rm
i}\left[\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}}
\lambda\dfrac{k_{x}k_{y}}{k^{2}} &
\lambda\left(1-\dfrac{k_{x}^{2}}{k^{2}}\right) & 0 &
-\alpha\dfrac{k_{x}k_{z}}{k^{2}}\\
\lambda\left({-}1+\dfrac{k_{y}^{2}}{k^{2}}\right) &
-\lambda\dfrac{k_{x}k_{y}}{k^{2}} & 0 &
-\alpha\dfrac{k_{y}k_{z}}{k^{2}}\\
\lambda\dfrac{k_{y}k_{z}}{k^{2}} &
-\lambda\dfrac{k_{x}k_{z}}{k^{2}} & 0 &
\alpha-\alpha\dfrac{k_{z}^{2}}{k^{2}}\\ 0 & 0 & 1 & 0
\end{array}\right]. \end{equation}
For the eigenequation
\begin{equation} \boldsymbol{\mathsf{H}}_{0}\boldsymbol{\psi}=\omega\boldsymbol{\psi}, \end{equation}the non-zero eigenvalues and corresponding right eigenvectors are
\begin{gather} \omega_{{\pm}}={\pm}\omega_{0},\quad \omega_{0}=\frac{\sqrt{k_{z}^{2}-\alpha(k_{x}^{2}+k_{y}^{2})}}{k}, \end{gather}
\begin{gather}
\left.\left|\psi_{{\pm}}\right.\right\rangle
=\frac{1}{\sqrt{2\left(\dfrac{k_{z}^{2}}{k_{x}^{2}+k_{y}^{2}}-\alpha\right)}}\left[\begin{array}{@{}c@{}}
\dfrac{k_{z}\left(-{\rm i}\lambda k_{y}\mp
k_{x}\omega_{0}\right)}{k_{x}^{2}+k_{y}^{2}}\\
\dfrac{k_{z}\left({\rm i}\lambda k_{x}\mp
k_{y}\omega_{0}\right)}{k_{x}^{2}+k_{y}^{2}}\\
\pm\omega_{0}\\ {\rm i} \end{array}\right],
\end{gather}
where 
$k=\sqrt {k_x^2+k_y^2+k_z^2}$ and 
$k_x^2+k_y^2\neq 0$. When 
$\rho =\sqrt {k_{x}^{2}+k_{y}^{2}}\xrightarrow {}0$, the eigenvectors become
\begin{equation} \left.\left|\psi_{{\pm}}\right.\right\rangle =\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} \pm \lambda {\rm i}\\ 1\\ 0\\ 0 \end{array}\right]. \end{equation}
When 
$\rho \xrightarrow {}+\infty$, the eigenvectors are still single valued so that
\begin{equation} \left.\left|\psi_{{\pm}}\right.\right\rangle =\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} 0\\ 0\\ \pm 1\\ \dfrac{{\rm i}}{\sqrt{-\alpha}} \end{array}\right]. \end{equation}Thus, different from the shallow water model, the eigenvectors here are regular on a compact manifold, which guarantees a well-defined Chern number without introducing unphysical items (Tauber, Delplace & Venaille Reference Tauber, Delplace and Venaille2019).
For the eigenequation
\begin{equation} \boldsymbol{\mathsf{H}}_{0}^{{{\dagger}}}\boldsymbol{\psi}^{\prime}=\omega^{{\ast}}\boldsymbol{\psi}^{\prime}, \end{equation}the non-zero eigenvalues and corresponding left eigenvectors are
\begin{gather} \omega_{{\pm}}^{{\ast}}={\pm}\omega_{0}^{{\ast}}, \end{gather}
\begin{gather} \left\langle
\left.\psi_{{\pm}}^{\prime}\right|\right.=\frac{1}{\sqrt{2\left(\dfrac{k_{z}^{2}}{k_{x}^{2}+k_{y}^{2}}-\alpha\right)}}\left[\begin{array}{@{}c@{}}
\dfrac{k_{z}\left({\rm i}\lambda
k_{y}\mp\dfrac{k_{x}}{\omega_{0}}\right)}{k_{x}^{2}+k_{y}^{2}}\\
\dfrac{k_{z}\left(-{\rm i}\lambda
k_{x}\mp\dfrac{k_{y}}{\omega_{0}}\right)}{k_{x}^{2}+k_{y}^{2}}\\
\mp\dfrac{\alpha}{\omega_{0}}\\ {\rm i}\alpha
\end{array}\right]^{\mathrm{T}},
\end{gather}
when 
$k_x^2+k_y^2\neq 0$. When 
$\rho \xrightarrow {}0$, the eigenvectors become
\begin{equation} \left\langle \left.\psi_{{\pm}}^{\prime}\right|\right. =\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} \mp \lambda {\rm i}\\ 1\\ 0\\ 0 \end{array}\right]^{\mathrm{T}}. \end{equation}
When 
$\rho \xrightarrow {}+\infty$, the eigenvectors are single valued so that
\begin{equation} \left\langle \left.\psi_{{\pm}}^{\prime}\right|\right. =\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} 0\\ 0\\ \pm 1\\ -{\rm i}\sqrt{-\alpha} \end{array}\right]^{\mathrm{T}}. \end{equation} Substituting the expressions of eigenvectors in (2.15), the Berry curvature of the positive band (associated with 
$\omega _+$ and 
$\omega _+^{\ast }$) becomes
\begin{equation} \varOmega^{z}={-}\lambda\frac{k_{z}^{2}\left[\rho^{2}k_{z}^{2}+2\alpha\left(\rho^{4}+4\rho^{2}k_{z}^{2}+2k_{z}^{4}\right)-\alpha^{2}\rho^{2}\left(2\rho^{2}+k_{z}^{2}\right)\right]}{2k^{3}\left(k_{z}^{2}-\alpha\rho^{2}\right)^{5/2}}, \end{equation}
where 
$\rho =\sqrt {k_{x}^{2}+k_{y}^{2}}$. Then the Chern number is
\begin{align} C&=\frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\varOmega^{z}\,{\rm d}k_{x}\,{\rm d}k_{y} \end{align}
\begin{align} &=\int_{0}^{+\infty}\rho\varOmega^{z}\,{\rm d}\rho \end{align}
\begin{align} &=\lambda\left.\frac{-2k_{z}^{4}+k_{z}^{2}(\alpha-1)\rho^{2}}{2\sqrt{\rho^{2}+k_{z}^{2}}\left(k_{z}^{2}-\alpha\rho^{2}\right)^{3/2}}\right|_0^{+\infty} \end{align}
\begin{align} &=\lambda, \end{align}which implies the existence of topologically protected edge states only related to the direction of fluid rotation. This is a three-dimensional version of the topologically protected equatorial waves (Delplace et al. Reference Delplace, Marston and Venaille2017). Under the hydrostatic approximation, we obtain an exact two-dimensional counterpart, shown in § 5.
3.3. Topological invariant when $E\neq 0$
This section considers 
$E\neq 0$, which does not change the Chern number we obtained in the previous subsection. Firstly, When 
$Pr=1$, the Hamiltonian matrix 
$H$ becomes
\begin{equation} \boldsymbol{\mathsf{H}}=H_{0}+{\rm i}E_{0}\boldsymbol{\hat{1}}, \end{equation}
where 
$H_{0}$ is the simple Hamiltonian where 
$E=0$ and 
$\boldsymbol {\hat {1}}$ is the unit matrix of the same order as 
$H_{0}$. The eigenvalues are
\begin{equation} \omega_{{\pm}}={\pm}\omega_{0}+{\rm i}E_{0}. \end{equation}
It is just a shift of 
$\pm \omega _{0}$, the eigenvalues of 
$H_{0}$ and then the eigenvectors are independent of 
$E$. Therefore, both the Berry curvature and Chern number are unchanged compared with the case where 
$E=0$.
When 
$Pr\neq 1$, there are no simple expressions for the eigenvalues and eigenvectors, and we need to calculate the Chern number with the help of numerical calculations, seen in figure 2. To avoid the derivation of the eigenvectors, we rewrite (2.15) as
\begin{equation} \varOmega_{n}^{z}={\rm
i}\underset{m\neq n}{\sum}\frac{\left\langle\boldsymbol{\psi}_{n}^{\prime}\left|\dfrac{\partial
H}{\partial
k_{x}}\right|\psi_{m}\right\rangle\left\langle\psi_{m}^{\prime}\left|\dfrac{\partial
H}{\partial k_{y}}\right|\boldsymbol{\psi}_{n}\right\rangle-\{x\leftrightarrow
y\}}{(\omega_{n}-\omega_{m})^{2}}.
\end{equation}
It is similar to the Hermitian case (Xiao et al. Reference Xiao, Chang and Niu2010) except that the left vector takes the biorthogonal partner. This form is very useful for numerical calculations because it can be evaluated under an unsmooth phase choice of the eigenstates, which often occurs in the standard diagonalization algorithms.
Figure 2. Numerical results of the non-Hermitian Chern number. (a) Chern number with different Prandtl numbers 
$Pr$ and rotation directions 
$\lambda =\pm 1$. The upper limit of integral for 
$\rho$ is 
$k_{cut}=6$, the 
$k$-point accuracy per unit length is 200 (radial) * 200 (axial) using the Gauss–Legendre integration scheme. (b) Demonstrating the numerical convergence degree with different 
$k$-point accuracies as 
$k_{cut}$ increases. Other parameters are 
$\alpha =1$, 
$E=0.01$, 
$k_{z}=1+0.1 {\rm i}$.
4. The stably stratified and critical cases
Even though this paper focuses on the unstably stratified convection situation, this section shows that, when the stratification is stable, following the above-mentioned procedure, we can obtain the same Chern number as in classic stably stratified situations (e.g. Delplace et al. Reference Delplace, Marston and Venaille2017). For the simple case that 
$E=0$, when 
$\alpha <0$, the eigenvalues of the linearized equations (2.2)–(2.4) are real, and the system is stable without oscillatory modes that can grow or decay over time. At this time the discussion in the previous section where 
$\alpha >0$ still holds, and the system still has a non-zero Chern number 
$C=\lambda$. Thus, the existence of the tenacious sidewall modes in the stratified fluid relies only on rotation, regardless of the driving mechanism.
When 
$\alpha =0$, the system is in a critical neutral state, and the eigenvalues of the Hamiltonian are
$$\begin{gather} \omega_{1}\equiv\omega_+ = \frac{k_{z}}{k}-{\rm i}Ek^{2}, \end{gather}$$
$$\begin{gather}\omega_{2}\equiv\omega_{-}={-}\frac{k_{z}}{k}-{\rm i}Ek^{2}, \end{gather}$$
$$\begin{gather}\omega_{3}={-}{\rm i}\frac{E}{Pr}k^{2}, \end{gather}$$
$$\begin{gather}\quad\omega_{4}={-}{\rm i}Ek^{2}. \end{gather}$$
Considering the symmetry between 
$\omega _+$ and 
$\omega _{-}$, we have assumed that 
${\rm Re}\{k_{z}\}>0$. The Berry curvature of band 
$\omega _{3}$ (or 
$\omega _{4}$) is 0. The Berry curvature of band 
$\omega _{\pm }$ reads
\begin{equation} \varOmega^{z}_{{\pm}}={\mp}\lambda\frac{k_{x}^{2}+k_{y}^{2}}{2k_{z}k^{3}}, \end{equation}
which is independent of 
$E$ and 
$Pr$. The integral does not converge when we evaluate the Chern number
\begin{equation} C_{{\pm}}=\frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\varOmega^{z}_{{\pm}}\,{\rm d}k_{x}\,{\rm d}k_{y} =\int_{0}^{+\infty}\frac{\mp\lambda\rho^{3}}{2k_{z}(\rho^2+k_{z}^2)^{3/2}}\,{\rm d}\rho. \end{equation}
Mathematically, this is due to the degeneration of 
$\omega _+$ and 
$\omega _{-}$ at 
$\rho \xrightarrow {}+\infty$, which leads to the closure of the band gap of the system. The critical case is important to show that stratification is crucial, whether stable or unstable.
5. Hydrostatic approximation
Under the hydrostatic approximation, the 
$z$-direction momentum equation of (2.2) reduces to
\begin{equation} \hat{\theta}=\frac{{\rm i}\hat{p} k_{z}}{\alpha}. \end{equation}From (2.3) we get
\begin{equation} \hat{w}={-}\frac{k_{x}\hat{u}+k_{y}\hat{v}}{k_{z}}. \end{equation} With the wave vector 
$\boldsymbol{\psi} \equiv (\hat {u},\hat {v},\hat {p})$, we get the eigenequations from (2.2) and (2.4) as
\begin{equation} \boldsymbol{\mathsf{H}}=\left[\begin{array}{@{}c@{\quad}c@{\quad}c@{}} {\rm i}E_{0} & {\rm i}\lambda & k_{x}\\ -{\rm i}\lambda & {\rm i}E_{0} & k_{y}\\ -\alpha\dfrac{k_{x}}{k_{z}^{2}} & -\alpha\dfrac{k_{y}}{k_{z}^{2}} & {\rm i}E_{1} \end{array}\right], \end{equation}
where 
$E_{0}=-Ek^{2}$, 
$E_{1}=E_{0}/Pr$. This is a non-Hermitian Hamiltonian because 
$H\neq H^{{\dagger} }$, but for the stable case where 
$\alpha <0$, it can be changed into a Hermitian one by replacing the variables. If we further assume that 
$Pr=1$, the Hamiltonian is an analogue of the one in topological shallow water waves in Delplace et al. (Reference Delplace, Marston and Venaille2017), except for a frequency shift. However, for the unstable case with 
$\alpha >0$, the non-Hermitian nature of the system is intrinsic and cannot be turned into a Hermitian one by variable substitutions.
For the eigenequation
\begin{equation} \boldsymbol{\mathsf{H}}\boldsymbol{\psi}=\omega\boldsymbol{\psi}, \end{equation}
the eigenvalues and eigenvectors contributing to non-zero Chen numbers are (
$Pr=1$)
\begin{gather} \omega_{{\pm}}={\pm}\omega_{0}-{\rm i}Ek^{2},\quad \omega_{0}=\sqrt{1-\frac{\alpha}{k_{z}^{2}}(k_{x}^{2}+k_{y}^{2})}, \end{gather}
\begin{gather} \left.\left|\psi_{{\pm}}\right.\right\rangle =\frac{1}{\sqrt{2\left(\dfrac{k_{z}^{2}}{k_{x}^{2}+k_{y}^{2}}-\alpha\right)}}\left[\begin{array}{@{}c@{}} \dfrac{k_{z}^{2}\left({\rm i}\lambda k_{y}\pm k_{x}\omega_{0}\right)}{k_{x}^{2}+k_{y}^{2}}\\ \dfrac{k_{z}^{2}\left(-{\rm i}\lambda k_{x}\pm k_{y}\omega_{0}\right)}{k_{x}^{2}+k_{y}^{2}}\\ -\alpha \end{array}\right].\end{gather}
When 
$\rho \xrightarrow {}0$, the eigenvectors become
\begin{equation} \left.\left|\psi_{{\pm}}\right.\right\rangle =\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} \pm \lambda {\rm i}\\ 1\\ 0 \end{array}\right]. \end{equation}
When 
$\rho \xrightarrow {}+\infty$, the eigenvectors become
\begin{equation}
\left.\left|\psi_{{\pm}}\right.\right\rangle
=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} \pm
\sqrt{k_z^2}\dfrac{k_x}{\rho}\\ \pm
\sqrt{k_z^2}\dfrac{k_y}{\rho}\\ \sqrt{-\alpha}
\end{array}\right], \end{equation}
which are not single valued in different directions. In order to get a compact manifold where the Chern number is well defined, we may need to introduce an odd viscosity term as in the shallow water model (Tauber et al. Reference Tauber, Delplace and Venaille2019).
For the eigenequation
\begin{equation} \boldsymbol{\mathsf{H}}^{{{\dagger}}}\boldsymbol{\psi}^{\prime}=\omega^{\prime}\boldsymbol{\psi}^{\prime}, \end{equation}we have
\begin{gather} \omega_{{\pm}}^{\prime}=\omega_{{\pm}}^{{\ast}}, \end{gather}
\begin{gather} \left\langle
\left.\psi_{{\pm}}^{\prime}\right|\right.=\frac{1}{\sqrt{2\left(\dfrac{k_{z}^{2}}{k_{x}^{2}+k_{y}^{2}}-\alpha\right)}}\left[\begin{array}{@{}c@{}}
\dfrac{-{\rm i}\lambda k_{y}\pm
k_{x}\omega_{0}}{k_{x}^{2}+k_{y}^{2}}\\ \dfrac{{\rm
i}\lambda k_{x}\pm k_{y}\omega_{0}}{k_{x}^{2}+k_{y}^{2}}\\
1
\end{array}\right]^{\mathrm{T}}.\end{gather}
When 
$\rho \xrightarrow {}0$, the eigenvectors become
\begin{equation} \left\langle \left.\psi_{{\pm}}^{\prime}\right|\right. =\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} \mp \lambda {\rm i}\\ 1\\ 0 \end{array}\right]^{\mathrm{T}}. \end{equation}
When 
$\rho \xrightarrow {}+\infty$, the eigenvectors become
\begin{equation} \left\langle
\left.\psi_{{\pm}}^{\prime}\right|\right.
=\frac{1}{\sqrt{2}}\left[\begin{array}{@{}c@{}} \pm
\dfrac{k_x}{\rho\sqrt{k_z^2}}\\ \pm
\dfrac{k_y}{\rho\sqrt{k_z^2}}\\ \dfrac{1}{\sqrt{-\alpha}}
\end{array}\right]. \end{equation}
The Berry curvature of the positive band (associated with 
$\omega _+$ and 
$\omega _+^{\prime }$) is
\begin{equation} \varOmega^{z}={-}\frac{\lambda\alpha\sqrt{k_{z}^{2}}}{\left(k_{z}^{2}-\alpha\rho^{2}\right)^{3/2}}, \end{equation}
where 
$\rho =\sqrt {k_{x}^{2}+k_{y}^{2}}$. Then the Chern number is
\begin{align} C&=\frac{1}{2{\rm \pi}}\int_{-\infty}^{+\infty}\int_{-\infty}^{+\infty}\varOmega^{z}\,{\rm d}k_{x}\,{\rm d}k_{y} \end{align}
\begin{align} &=\int_{0}^{+\infty}\rho\varOmega^{z}\,{\rm d}\rho \end{align}
\begin{align} &=\left.-\frac{\lambda\sqrt{k_{z}^{2}}}{\sqrt{k_{z}^{2}-\alpha\rho^{2}}}\right|_0^{+\infty} \end{align}
\begin{align} &=\lambda. \end{align}
In the stably stratified case, 
$\alpha <0$, the above calculation holds when 
${\rm Re}\{k_{z}\}\neq 0$; while, in the unstably stratified case, 
$\alpha >0$, we need 
${\rm Im}\{k_{z}\}\neq 0$. If the system is neutral, 
$\alpha =0$, we arrive at two gapped flat bands 
$\omega _{\pm }=\pm 1-{\rm i}Ek^{2}$, which differ from the degenerate scenario in (4.1), and obtain 
$C=0$ since 
$\varOmega ^{z}=0$.
6. Winding number
To put the topological nature of the system into perspective, we define a complex function in the wavenumber space 
$(k_x,k_y)$ (cf. Zhu et al. Reference Zhu, Li and Marston2023)
\begin{equation} \varXi(\boldsymbol{k})=\hat{v}(\boldsymbol{k})\hat{w}^{{\ast}}(\boldsymbol{k}), \end{equation}
where 
$\hat {v}$ and 
$\hat {w}$ are the second and third components of 
$\psi (\boldsymbol {k})$. Using other component combinations to define 
$\varXi (\boldsymbol {k})$ should also be fine. As long as one factor is taken as its complex conjugate, we can remove the gauge redundancy of the eigenfunctions, leaving only the internal phase difference between the two components.
Figure 3 depicts the argument of 
$\varXi (\boldsymbol {k})$, 
$\tan ^{-1}({\rm Re}\{\varXi \}/{\rm Im}\{\varXi \})$, whose 
$x$ and 
$y$ components represent the real and imaginary parts of 
$\varXi$ with a rescaled equal length. When 
$\lambda =1$, there is a vortex formed by the vector arrows going in the anticlockwise direction, and the arrows along a closed circle smoothly wind by a phase of 2
${\rm \pi}$, suggesting a winding number of 
$1$. Conversely, when 
$\lambda =-1$, there is a strain flow with the winding number of 
$-1$.
Figure 3. Arrows representing argument of 
$\varXi (\boldsymbol {k})=\hat {v}(\boldsymbol {k})*\hat {w}^{\ast }(\boldsymbol {k})$ with 
$\lambda =1$ (a) and 
$\lambda =-1$ (b). The 
$x$ and 
$y$ components of the arrow represent the real and imaginary parts of 
$\varXi (\boldsymbol {k})$. The length of the arrows is rescaled to be equal. Colours represent normalized magnitude 
$|\varXi (\boldsymbol {k})|$ in arbitrary units. Other parameters are 
$Pr=1$, 
$\alpha =0.5$, 
$k_{z}=1+0.5 {\rm i}$.
7. Summary and discussion
In summary, we show that the linearized rotating Rayleigh–Bénard convection may support non-Hermitian topological states characterized by a non-zero integer Chern number. This finding is consistent with the numerical prediction by Favier & Knobloch (Reference Favier and Knobloch2020) that the robust sidewall modes in rapidly rotating convection are topologically protected. Due to the unstable stratification, the linear eigenvalue problem is intrinsically non-Hermitian, so the Berry curvature is defined on biorthogonal eigenstates, and the corresponding eigenvalues may be complex, which is very different from the Hermitian system. Beyond the shallow water model, the eigenvectors here are regular on a compact manifold, which guarantees a well-defined Chern number without introducing unphysical items (Tauber et al. Reference Tauber, Delplace and Venaille2019). The emergence of these topological edge states is fundamentally due to rotation breaking the system's time-reversal symmetry, but stratification also plays a crucial role: without the stratification, i.e. for the critical case, the bulk Chern number is either not well defined or equals zero. Under the hydrostatic approximation, the problem transforms into a two-dimensional counterpart of the one that explains the topological origin of equatorial waves (Delplace et al. Reference Delplace, Marston and Venaille2017). Finally, an eigenvector-dependent winding number is introduced to visualize the topological nature of the fluid. Our conclusions hold within the specific constraints of our model, which include no-slip boundary conditions at the top/bottom surface, and smooth walls on the sides. However, it should be noted that we do not provide a direct proof of the explicit connection between the non-zero Chern number and robust edge modes, although we can invoke the Atiyah–Patodi–Singer index theorem (Fukaya, Onogi & Yamaguchi Reference Fukaya, Onogi and Yamaguchi2017; Yu, Wu & Xie Reference Yu, Wu and Xie2017) to support our claims, as is common in typical condensed matter physics articles. We acknowledge deeper exploration is needed to clarify this bulk-boundary correspondence, which will extend our conclusions to more general boundary conditions.
Topological invariants are mathematical quantities that remain unchanged under continuous deformations of a system. They can both qualitatively describe the topological properties of a system (determining whether the robust sidewall states exist or not), and quantitatively describe the system's proximity to topological phase transitions, in other words, the robustness level of the sidewall states. As the system parameters are tuned, the topological invariants may change abruptly at certain critical points, signalling the occurrence of a topological phase transition. The size of the energy gap, which is closely related to the topological invariants, can provide a concrete metric for this proximity (Xia et al. Reference Xia2009). A larger energy gap indicates that the system is further away from the critical point, while a smaller gap suggests the system is closer to the transition and more susceptible to the emergence of interesting topological phenomena.
With the conclusions drawn from topology, we can avoid the secondary dynamical details and better understand the fundamental influencing factors of the problem. In our example, we suggested that the key factors for the existence of robust sidewall modes are the presence of rotation and stratification, rather than the heat distribution, whether the stratification is stable, or whether the boundaries are flat, etc. This is what previous numerical calculations and analytical derivations could not achieve.
On the precession direction of the edge states, they are usually (not all) retrograde in rotating convection (Zhong et al. Reference Zhong, Ecke and Steinberg1991; Goldstein et al. Reference Goldstein, Knobloch, Mercader and Net1993; Herrmann & Busse Reference Herrmann and Busse1993; Kuo & Cross Reference Kuo and Cross1993), but prograde both at the Earth's equator and in the two-dimensional fluids with odd viscosity (Delplace et al. Reference Delplace, Marston and Venaille2017; Souslov et al. Reference Souslov, Dasbiswas, Fruchart, Vaikuntanathan and Vitelli2019). This is a result of the combined contribution of the boundary conditions and the system parameters, especially the effect of the Prandtl number 
$Pr$ (Goldstein et al. Reference Goldstein, Knobloch, Mercader and Net1994; Horn & Schmid Reference Horn and Schmid2017). It is quite different from the conventional Hermitian system with periodic boundary conditions, where the positive or negative sign of the Chen number determines the direction of the edge states (Hasan & Kane Reference Hasan and Kane2010; Essin & Gurarie Reference Essin and Gurarie2011).
The bulk-interface correspondence for equatorial waves was established by Tauber et al. (Reference Tauber, Delplace and Venaille2019) and Venaille & Delplace (Reference Venaille and Delplace2021). There, the topological edge waves are regarded as the interface between two adjacent bulks. It is generally applicable for the rotating convection case, when we consider the sidewall (edge) modes as interfaces between the internal fluid (one bulk) and the sidewalls (another bulk). This reflects the standard procedure when dealing with topologically protected edge states. According to the principles of topology, when two bulk materials with different topological invariants come into contact, there will always be an interface between them, which is gapless and lacks a well-defined topological invariant. However, one point that needs to be noted is that we are currently dealing with the bulk-edge correspondence of a non-Hermitian system, rather than simply systems with periodic boundaries. This is an important generalization of bulk-edge correspondence (Imura & Takane Reference Imura and Takane2019; Xiao et al. Reference Xiao, Deng, Wang, Zhu, Wang, Yi and Xue2019; Yokomizo & Murakami Reference Yokomizo and Murakami2020; Zirnstein, Refael & Rosenow Reference Zirnstein, Refael and Rosenow2021; Rapoport & Goldstein Reference Rapoport and Goldstein2023).
For simplicity, we assume both 
$k_{x}$ and 
$k_{y}$ to be real numbers. Thus, according to the regular formalities (Souslov et al. Reference Souslov, Dasbiswas, Fruchart, Vaikuntanathan and Vitelli2019), the non-zero Chern number we obtain describes the topological properties of the system without boundary in the 
$x-y$ plane. When there is an impenetrable and slippery sidewall, the bulk-boundary correspondence ensures the existence of topologically protected edge states (Essin & Gurarie Reference Essin and Gurarie2011). When the sidewall is no slip, we can still get a well-defined non-zero Chern number, but in this case, 
$k_{x}$ and 
$k_{y}$ need to be taken as complex numbers (Yao et al. Reference Yao, Song and Wang2018). Despite the significant influence of boundary conditions on the topological properties of non-Hermitian systems, we can continue to handle the issue using the procedures under periodic boundary conditions, thanks to the generalized Brillouin zone with complex wavenumbers (Yokomizo & Murakami Reference Yokomizo and Murakami2019; Yang et al. Reference Yang, Zhang, Fang and Hu2020; Wu et al. Reference Wu, Xie, Zhou and An2022).
Our definition of the non-Hermitian Berry curvature is quite traditional, but is not necessarily the only way. There may be different ways of defining it, considering only left eigenvalues or only right eigenvalues or switching their roles, but that does not affect the values of their integrals, i.e. the Chern numbers (Shen et al. Reference Shen, Zhen and Fu2018).
Rayleigh–Bénard convection is essentially a nonlinear system, so as a complement to the conclusions of this paper, it is necessary to continue to investigate the nonlinear effects on the topological properties of the bulk states. We expect that weak nonlinear effects do not destroy the topological invariance of the system, partly because the topological properties are robust to perturbations and partly because previous numerical simulations imply this (Favier & Knobloch Reference Favier and Knobloch2020). For a more detailed quantitative analysis, we need to evaluate the generalized geometric phase in the case of non-eigenstates (Liu, Wu & Niu Reference Liu, Wu and Niu2003; Wu, Liu & Niu Reference Wu, Liu and Niu2005).
Funding
This work has received financial support from the National Natural Science Foundation of China (NSFC) under grant nos. 92052102 and 12272006, the Laoshan Laboratory under grant no. 2022QNLM010201 and from the Fundamental Research Funds for the Central Universities under grant no. BLX202348.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Solutions satisfying realistic boundary conditions
Unlike in Hermitian systems, the eigenfrequencies and solutions of non-Hermitian systems are very sensitive to the boundary conditions. For simplicity, we assume that the horizontal direction is unbounded and the 
$z$-direction takes the no-slip boundary condition such that
\begin{equation} u=v=w=\theta=0, \quad \text{on} \ z=0,1. \end{equation}To ensure that our analyses above match a realistic situation, we construct solutions in the form that
\begin{equation} \boldsymbol{\phi}=\exp({{\rm i}(k_x x+k_y y-\omega t)})\sum_{i=1}^4 \lambda_i {\rm e}^{{\rm i}k_{zi} z}\left.\left|\psi_{i}\right.\right\rangle, \end{equation}
where 
$|\psi _{i}\rangle$ is the eigenvector of Hamiltonian (3.15) corresponding to 
$k_{zi}$. According to our calculations above, when the system evolves along a closed path in its parameter space, 
$|\psi _{1-4}\rangle$ will all get the Berry phase of 
$2{\rm \pi} \lambda$. As a result, the solution 
$\phi$ also gets the same Berry phase. For a particular frequency 
$\omega (k_x,k_y)$, 
$k_{z1-4}$ come from the dispersion relation that (
$Pr=1$)
\begin{equation} \omega=\frac{\sqrt{k_{z}^{2}-\alpha(k_{x}^{2}+k_{y}^{2})}}{k}-{\rm i}E k^2. \end{equation}
Generally, the solved 
$k_{zi}$ is complex, and 
$k_{z3}=-k_{z1}$, 
$k_{z4}=-k_{z2}$.
For the no-slip boundary condition, we have
\begin{equation} \boldsymbol{\mathsf{M}}\left[\begin{array}{@{}c@{}} \lambda_1\\ \lambda_2\\ \lambda_3\\ \lambda_4 \end{array}\right]=0, \end{equation}where
\begin{equation} \boldsymbol{\mathsf{M}}=\left[\begin{array}{@{}c@{\quad}c@{\quad}c@{\quad}c@{}} 1 & 1 & 1 & 1\\ k_{z1} & k_{z2} & k_{z3} & k_{z4}\\ {\rm e}^{{\rm i}k_{z1}} & {\rm e}^{{\rm i}k_{z2}} & {\rm e}^{{\rm i}k_{z3}} & {\rm e}^{{\rm i}k_{z4}}\\ {\rm e}^{{\rm i}k_{z1}}k_{z1} & {\rm e}^{{\rm i}k_{z2}}k_{z2} & {\rm e}^{{\rm i}k_{z3}}k_{z3} & {\rm e}^{{\rm i}k_{z4}}k_{z4} \end{array}\right]. \end{equation}
To get a non-trivial solution of 
$\lambda _i(k_x,k_y)$, we have
\begin{equation} {\rm det}(\boldsymbol{\mathsf{M}})=0, \end{equation}and obtain that
\begin{equation} 2 k_{z1} k_{z2} ({-}1 + \cos{k_{z1}} \cos{k_{z2}}) + (k_{z1}^2 + k_{z2}^2) \sin{k_{z1}} \sin{k_{z2}}=0. \end{equation}
Obviously, 
$k_{z1}=k_{z2}$ satisfies the above equation, but this will give a trivial solution that 
$\lambda _i(k_x,k_y)=0$ and we do not choose it. Combining with (A3), one can identify the non-trivial values of 
$k_{zi}$, 
$\omega (k_x,k_y)$ and 
$\lambda _i(k_x,k_y)$. As an example, for parameters 
$k_x=k_y=0.5$, 
$\alpha =1$ and 
$E=0.01$, numerical calculation gives that 
$\omega \approx 0.9903-0.8108{\rm i}$, 
$k_{z1}\approx 8.9763-0.0195{\rm i}$, 
$k_{z2}\approx 0.3998-0.6643{\rm i}$ and
\begin{equation} \left[\begin{array}{@{}c@{}} \lambda_1\\ \lambda_2\\ \lambda_3\\ \lambda_4 \end{array}\right] =\mathcal{C} \left[\begin{array}{@{}c@{}} 1\\ 0.9877-1.5056{\rm i}\\ 0.9188-0.4422{\rm i}\\ -2.9066+1.9478{\rm i} \end{array}\right], \end{equation}
where 
$\mathcal {C}$ is an arbitrary constant.
