Non-Hermitian Chern number in rotating Rayleigh–Bénard convection

Abstract

We show that rotating Rayleigh–Bénard convection, where a rotating fluid is heated from below, exhibits a non-Hermitian topological invariant. Recently, Favier & Knobloch (J. Fluid Mech., vol. 895, 2020, R1) hypothesized that the robust sidewall modes in rapidly rotating convection are topologically protected. By considering a Berry curvature defined in the complex wavenumber space, we reveal that the bulk states can be characterized by a non-zero integer Chern number, implying a potential topological origin of the edge modes based on the Atiyah–Patodi–Singer index theorem (Fukaya et al., Phys. Rev. D, vol. 96 2017, 125004; Yu et al., Nucl. Phys. B, vol. 916, 2017, pp. 550–566). The linearized eigenvalue problem is intrinsically non-Hermitian, therefore, the definition of Berry curvature generalizes that of the stably stratified problem. Moreover, the three-dimensional set-up naturally regularizes the eigenvector, avoiding the compactification problem in shallow water waves (Tauber et al., J. Fluid Mech., vol. 868, 2019, R2). Under the hydrostatic approximation, it recovers a two-dimensional analogue of the one which explains the topological origin of the equatorial Kelvin and Yanai waves (Delplace et al., Science, vol. 358, issue 6366, 2017, pp. 1075–1077). The non-zero Chern number relies only on rotation when the fluid is stratified, no matter whether it is stable or unstable. However, the neutrally stratified system does not support a topological invariant. In addition, we define a winding number to visualize the topological nature of the fluid. Our results represent a step forward for the topologically protected states in convection, but the bulk-boundary correspondence requires a further direct analysis for proof, and the robustness of the edge states under varying boundary conditions remains a question to be answered.

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 2000; Ahlers, Grossmann & Lohse 2009; King, Stellmach & Aurnou 2012; Guervilly, Hughes & Jones 2014). 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 1991, 1993; Ecke, Zhong & Knobloch 1992; Goldstein et al. 1993; Ning & Ecke 1993), 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 2020; Ecke 2023; Ecke & Shishkina 2023). These states were first discovered through experiments indirectly by Rossby (1969), where a surprising convection occurs with a Rayleigh number ($Ra$) below the critical value for the laterally unbounded system (Chandrasekhar 1961). Zhang & Liao (2009) give the asymptotic expression of the critical $Ra$ for the onset of these sidewall modes, regaining the leading terms previously calculated by Herrmann & Busse (1993). 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. 1993; Liu & Ecke 1999; Horn & Schmid 2017). 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. 2020; Zhang et al. 2020; Ecke, Zhang & Shishkina 2022). 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 2020). 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. 2009; Terrien, Favier & Knobloch 2023).

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 2010; Delplace, Marston & Venaille 2017) or the winding number (Zhu, Li & Marston 2023). When the topological invariant is non-zero, the bulk-boundary correspondence guarantees the existence of robust edge states (Essin & Gurarie 2011; Venaille et al. 2023; Onuki, Venaille & Delplace 2024). The topological physics first appeared in condensed matter (Kosterlitz & Thouless 1973; Thouless et al. 1982; Haldane 1988), such as topological insulators (Hasan & Kane 2010), topological photonics (Ozawa et al. 2019) and topological phononics (Liu, Chen & Xu 2020). 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. 2017; Venaille & Delplace 2021), the topological waves in fluids with odd viscosity (Souslov et al. 2019) and the presence of topological invariants in active matter systems (Shankar et al. 2022). This universality of topological states can be attributed to the similar symmetries embodied in systems of different scales (Senthil 2015).

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 2020). 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 2021; Ding, Fang & Ma 2022). 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 2020)

(2.1a)$$\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}$$

(2.1b)$$\begin{gather}\boldsymbol{\nabla}\boldsymbol{\cdot}\boldsymbol{u}=0, \end{gather}$$

(2.1c)$$\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

(2.2)\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}

(2.3)\begin{gather} \boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{\hat{u}}=0, \end{gather}

(2.4)\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

(2.5)\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

(2.6)\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

(2.7)\begin{equation} \boldsymbol{\mathsf{H}}\boldsymbol{\psi}=\omega\boldsymbol{\psi}, \end{equation}

with

(2.8) \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

(2.9)\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 1989; Xiao et al. 2010). 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 1984; Perez, Delplace & Venaille 2021).

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 2010). The bulk-boundary correspondence principle links the topological properties of the bulk and edge states (Essin & Gurarie 2011; Venaille et al. 2023; Onuki et al. 2024).

In the wavenumber space, the Berry connection (Berry 1984) is defined as

(2.10)\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

(2.11)\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

(2.12)\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 2019; Yang et al. 2020; Wu et al. 2022). As a result, a non-Bloch Chern number is defined as the integral of the Berry curvature over the GBZ

(2.13)\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 2018).

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

(2.14a)$$\begin{gather} \boldsymbol{\mathsf{H}}|\boldsymbol{\psi}_{n}\rangle=\omega_{n}|\boldsymbol{\psi}_{n}\rangle, \end{gather}$$

(2.14b)$$\begin{gather}\boldsymbol{\mathsf{H}}^{{{\dagger}}}|\psi^{\prime}_{m}\rangle=\omega_{m}^{{\ast}}|\psi^{\prime}_{m}\rangle, \end{gather}$$

(2.14c)$$\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. 2018)

(2.15)\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 (2018) 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 2019; Yang et al. 2020; Wu et al. 2022).

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

(3.1)\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 2009), 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. 2020; Zhang, Yang & Fang 2022; Zhang et al. 2022), 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

(3.2) \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

(3.3)\begin{equation} \boldsymbol{\mathsf{H}}_{0}\boldsymbol{\psi}=\omega\boldsymbol{\psi}, \end{equation}

the non-zero eigenvalues and corresponding right eigenvectors are

(3.4)\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}

(3.5) \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

(3.6)\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

(3.7)\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 2019).

For the eigenequation

(3.8)\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

(3.9)\begin{gather} \omega_{{\pm}}^{{\ast}}={\pm}\omega_{0}^{{\ast}}, \end{gather}

(3.10) \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

(3.11)\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

(3.12)\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

(3.13)\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

(3.14a)\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}

(3.14b)\begin{align} &=\int_{0}^{+\infty}\rho\varOmega^{z}\,{\rm d}\rho \end{align}

(3.14c)\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}

(3.14d)\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. 2017). 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

(3.15)\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

(3.16)\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

(3.17) \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. 2010) 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. 2017). 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

(4.1a)$$\begin{gather} \omega_{1}\equiv\omega_+ = \frac{k_{z}}{k}-{\rm i}Ek^{2}, \end{gather}$$

(4.1b)$$\begin{gather}\omega_{2}\equiv\omega_{-}={-}\frac{k_{z}}{k}-{\rm i}Ek^{2}, \end{gather}$$

(4.1c)$$\begin{gather}\omega_{3}={-}{\rm i}\frac{E}{Pr}k^{2}, \end{gather}$$

(4.1d)$$\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

(4.2)\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

(4.3)\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

(5.1)\begin{equation} \hat{\theta}=\frac{{\rm i}\hat{p} k_{z}}{\alpha}. \end{equation}

From (2.3) we get

(5.2)\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

(5.3) \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. (2017), 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

(5.4)\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$)

(5.5)\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}

(5.6)\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

(5.7)\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

(5.8) \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. 2019).

For the eigenequation

(5.9)\begin{equation} \boldsymbol{\mathsf{H}}^{{{\dagger}}}\boldsymbol{\psi}^{\prime}=\omega^{\prime}\boldsymbol{\psi}^{\prime}, \end{equation}

we have

(5.10)\begin{gather} \omega_{{\pm}}^{\prime}=\omega_{{\pm}}^{{\ast}}, \end{gather}

(5.11) \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

(5.12)\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

(5.13) \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

(5.14)\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

(5.15a)\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}

(5.15b)\begin{align} &=\int_{0}^{+\infty}\rho\varOmega^{z}\,{\rm d}\rho \end{align}

(5.15c)\begin{align} &=\left.-\frac{\lambda\sqrt{k_{z}^{2}}}{\sqrt{k_{z}^{2}-\alpha\rho^{2}}}\right|_0^{+\infty} \end{align}

(5.15d)\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. 2023)

(6.1)\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 (2020) 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. 2019). 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. 2017). 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 2017; Yu, Wu & Xie 2017) 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. 2009). 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. 1991; Goldstein et al. 1993; Herrmann & Busse 1993; Kuo & Cross 1993), but prograde both at the Earth's equator and in the two-dimensional fluids with odd viscosity (Delplace et al. 2017; Souslov et al. 2019). 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. 1994; Horn & Schmid 2017). 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 2010; Essin & Gurarie 2011).

The bulk-interface correspondence for equatorial waves was established by Tauber et al. (2019) and Venaille & Delplace (2021). 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 2019; Xiao et al. 2019; Yokomizo & Murakami 2020; Zirnstein, Refael & Rosenow 2021; Rapoport & Goldstein 2023).

For simplicity, we assume both $k_{x}$ and $k_{y}$ to be real numbers. Thus, according to the regular formalities (Souslov et al. 2019), 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 2011). 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. 2018). 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 2019; Yang et al. 2020; Wu et al. 2022).

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. 2018).

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 2020). For a more detailed quantitative analysis, we need to evaluate the generalized geometric phase in the case of non-eigenstates (Liu, Wu & Niu 2003; Wu, Liu & Niu 2005).