text:
\section*{Abstract}
{\footnotesize
This lecture note adresses the correspondence between spectral flows, often associated to unidirectional modes, and Chern numbers associated to degeneracy points. The notions of topological indices (Chern numbers, analytical indices) are introduced for non specialists with a wave physics or condensed matter background. The correspondence is detailed with several examples, including the Dirac equations in two dimensions, Weyl fermions in three dimensions, the shallow water model and other generalizations.  }

\tableofcontents

\section{Introduction}

In his successful popular book  \underline{La valeur de la science} \cite{Poincare},  the mathematician physicist Henri Poincar\'e pointed out in 1905 that : ''\textit{we can see mathematical analogies between phenomena which have no physical relation neither apparent nor real, so that the laws of one of these phenomena help us to guess those of the other. [...] The goal of mathematical physics is not only to facilitate the physicist in numerical computation [...]. It is still, it is above all to make the physicists know the hidden harmony of things by showing them from a new angle.}'' As one reads these words today, it is striking how much they make one think about the role of topology in physics.

Topology spread in virtually all branches of physics as it helps unveiling the hidden connections between apparently different phenomena.
In condensed matter physics, its use has been greatly accelerated since the rise of topological insulators \cite{Hasan_RMP, Qi_RMP}, one century after Poincar\'e's words.  In that context, topology  deals with the unavoidable singular properties of electronic wavefunctions associated to the energy bands of  solids. Those singularities are classified by integer numbers, called \textit{topological invariants}, that allow physicists to classify band insulators, superconductors and semimetals, beyond the realm of phase transitions of statistical physics \cite{Schnyder2009,Kitaev2009,Chiu2014,Yang2014,Chiu2015}. Although those topological properties are somehow \textit{hidden}, as they appear in a quite abstract space, their physical consequences are observable and manifest for instance by a robust uni-directional transport.

One the examples chosen by Poincar\'e to illustrate his statement is the wave equation, as it establishes a deep connection between sound, light and radio communications. It is thus not so surprising that the topological approaches that were developed in the realm of single particles quantum physics,  which is essentially wave physics, naturally spread to various classical waves systems, in optics, acoustics, mechanics, plasmas and fluids,  even at  geo/astrophysical scales \cite{Goldman2014,TopoPhot,Zhang2018,Xin2020,Shankar2020,Parker2021,Delplace2017}. All those systems share the same \textit{hidden harmony} (to quote Poincar\'e), that topology reveals by making us see those very different systems from a common angle. Topology thus basically provided a new perspective in the description of waves and a powerful tool for their manipulation.

This manuscript is dedicated to a key concept in wave topology, that is common to various quantum and classical systems.  It was actually already put forward by Volovik to draw formal analogies between $^3$He-A superfluids and the standard model in high energy physics \cite{VolovikBook}. Here we want to present it in a more wave physics / quantum condensed matter oriented mind. We shall refer to that concept as the \textit{monopole - spectral flow correspondence}. In short, this correspondence relates the topological properties of point defects in three dimensional (3D) space in which eigenstates of the system are parametrized, to the existence of peculiar modes that transit from an energy (or frequency) branch to another one; a phenomenon which is referred to as a \textit{spectral flow}. Those point defects correspond to degeneracy points (or band crossing points), that act as the source of a \textit{Berry curvature}, and whose flux through a close surface surrounding them is quantized and expressed by a topological invariant, the \textit{first Chern number}.   For that reason, we shall refer to these special points as  \textit{Berry-Chern} monopoles. The spectral flow, on the other hand, consists by construction to modes that propagate in a fixed direction. The monopole-spectral flow correspondence is a direct link between topology and unidirectional wave transport.

The goal of this manuscript is to introduce this correspondence in a pedagogical way, and make it accessible and useful to non specialists. In a first part, we focus on canonical examples encountered in condensed matter physics that we treat in a unified way. We use a  minimal two-band model  to illustrate this ubiquitous correspondence and describe simultaneously interface states in 2D systems such as Chern insulators, quantum valley Hall insulators and other  classical waves systems, together with the dispersive Landau  levels that trigger the chiral anomaly in 3D Weyl semimetals  \cite{Zyuzin2012, Burkov2015}.
This first part is mainly dedicated to a broad community of quantum and classical physicists that are not familiar with the topological concepts. Basic notions such as the Berry curvature, the first Chern number and the Weyl quantization are introduced, detailed and illustrated. The second part of this manuscript is dedicated to a far richer  model, first proposed by Ezawa  in the context of semimetals in a magnetic field \cite{Ezawa2017}, that allows for an investigation of the correspondence beyond the standard linear two-band crossing framework discussed in the first part. We discuss in details the structure of the spectral flow and its associated eigenmodes. In particular, this model is based on a spin algebra rather than a Clifford algebra, and thus goes beyond the usual framework of Dirac operators where this correspondence has been extensively studied (see e.g. \cite{APS1975, Niemi1986, Fukui2012, Yu2017, Bal2021, Zhao21} and the textbook \cite{NakaharaBook}). We thus introduce an  index for this ''non-Dirac'' operator, that accounts for the spectral flow of the model, and we compute it. We also  derive a general expression for the Chern numbers, and show by an explicit calculation that it is equal to the index we have introduced. This makes clear the intervention of the Atiyah-Singer index theorem to understand topological interface modes and 3D current of Landau levels associated to anomalies beyond two-band models. The unusual Landau levels behavior of effective spin$-J$ quasiparticles protected by crystal symmetries  \cite{Bradlyn2016, Ezawa2017} and the two eastward topological equatorial oceanic waves \cite{Delplace2017} are both direct examples accounted by this model.

\section{ Spectral flows in two-band models from archetypal examples in condensed matter physics }

Let us first revisit two canonical condensed matter examples where spectral flows emerge : the 2D Dirac Hamiltonian with an anisotropic mass term, and the 3D Weyl Hamiltonian under a constant axial magnetic field.

\subsection{2D Dirac fermions with an anisotropic mass}

The simplest and certainly the most ubiquitous physical model that manifests a spectral flow is that of a 2D massive Dirac Hamiltonian
\begin{align}
H_D = c \begin{pmatrix}
m(x)c & {\hat{p}}_x - {\rm i} {\hat{p}}_y \\
{\hat{p}}_x +{\rm i} {\hat{p}}_y  & -m(x)c
\end{pmatrix}
\label{eq:Dirac2D}
\end{align}
where $c$ is a celerity (e.g. the Fermi velocity) and $\hat{\mathbf{p}}=({\hat{p}}_x,{\hat{p}}_y)$ is the momentum operator. Importantly here, the masse term $m(x)$ varies in space along the $x$ direction, and is assumed to change sign at a given $x_0$ that we choose to be $x_0=0$. For simplicity, we shall  linearize $m$ around $m(0)=0$ so that $m(x) = \beta x$ with $\beta\equiv \partial_xm$. This is for convenience, as the properties we aim at describing do not depend on this linearization.
The Hamiltonian \eqref{eq:Dirac2D} describes e.g. the band inversion that occurs at the interface between two topologically insulators when tuning a parameter that regulates the amplitude of the gap, that is the mass term \cite{Haldane1988,BernevigBook}.

 Then, because of the $x$ dependence of the mass, $p_x$ is not a good quantum number, unlike $p_y$. One can then rewrite the Dirac Hamiltonian in position representation  for plane waves of wavenumber $k_y$ in the $y$ direction by substituting ${\hat{p}}_x=-{\rm i} \hbar \partial_x$, as
\begin{align}
H_D = c \begin{pmatrix}
c\beta x  & -{\rm i} \hbar \partial_x - {\rm i} \hbar k_y \\
-{\rm i}\hbar \partial_x  + {\rm i} \hbar k_y  & - c \beta x
\end{pmatrix} \ .
\end{align}

We aim at deriving the spectrum of this Hamiltonian as a function of the wavenumber $k_y$. For that purpose, it is convenient to first apply a unitary transformation $\tilde{H}_D = R H_D R^{-1}$ so that $k_y$ becomes proportional to the diagonal Pauli matrix $\sigma_z$. Writing $H_D=-{\rm i} \hbar c \partial_x\sigma_x +\hbar c k_y \sigma_y + \beta c^2 x \sigma_z $ where the $\sigma_i$'s are the Pauli matrices, such a unitary transformation $R$ amounts to a spin rotation where each axis are cyclically interchanged, and thus
\begin{align}
\tilde{H}_D =  \begin{pmatrix}
\hbar c k_y & \beta c^2 x - \hbar c \partial_x \\
\beta c^2 x  + \hbar c \partial_x & - \hbar c k_y
\end{pmatrix} \ .
\end{align}
Let us then introduce the characteristic length $\ell \equiv \sqrt{\hbar/\beta c}$ that allows us to rewrite the Dirac Hamiltonian in terms of annihilation and creation operators
\begin{align}
\label{eq:ann1}
{\hat{a}}  =\frac{1}{\sqrt{2}}\left( \frac{x}{\ell} + \ell \partial_x \right)\, , \qquad
{\hat{a}^{\dagger}} =\frac{1}{\sqrt{2}}\left( \frac{x}{\ell} - \ell \partial_x \right)
\end{align}
that satisfy $[a,a^\dagger ]=1$ and act on \textit{number states} $\ket{n}$ as
  \begin{align}
  {\hat{a}} \ket{n}  =\sqrt{n} \ket{n-1},\qquad {\hat{a}^{\dagger}} \ket{n}  =\sqrt{n+1} \ket{n+1}, \qquad  {\hat{a}^{\dagger}} {\hat{a}} \ket{n} = n \ket{n} \ .
  \end{align}
 We recall that number states are related to Hermite functions $\varphi_n(x)$ following
 \begin{align}
 \braket{x|n} = \varphi_n(x) = \frac{1}{(2^n n! \sqrt{\pi})^{1/2}} {\rm e}^{-\frac{x^2}{2}} H_n(x)
 \end{align}
such that the mode $n \in \mathbb{N}$ corresponds to the number of zeros of  $\varphi_n(x)$. The formalism of number states $\ket{n}$ will be used here for convenience, but note that it is not restrictive to quantum physics, as everything can equivalently be written in terms of Hermite functions $\varphi_n(x)$. As recalled in figure \ref{fig:Hermite}, those functions  are centered around $x=0$ (i.e. $\braket{x}_n=0$) that is where the mass term $m$ changes sign, but their spreading increases with $n$ as $\braket{x^2}_n = \frac{1}{2}+n$. Therefore, the lower the mode, the better the confinement.
\begin{figure}[h!]
  \centering
\includegraphics[width=6cm]{fig_Hermite}\qquad
\includegraphics[width=6cm]{figure_flow_Dirac}
\caption{\footnotesize{(Left) Spatial dependence of a few Hermite functions. (Right) Spectrum of $\tilde{H}_D$ and of $\tilde{H}_W$. The spectral flow (red) between the two branches correspond to the eigenstate \eqref{eq:anzatz_dirac}.  }}
\label{fig:Hermite}
\end{figure}

The Dirac Hamiltonian finally reads
\begin{align}
\tilde{H}_D =  \mathcal{E}
\begin{pmatrix}
\lambda & {\hat{a}^{\dagger}} \\
{\hat{a}}         & -\lambda
\end{pmatrix}
\label{eq:H_Dop}
\end{align}
where we have introduced the characteristic energy $\mathcal{E}\equiv \sqrt{2}\hbar c/\ell$ and the dimensionless parameter $\lambda \equiv \ell k_y / \sqrt{2}$.
The expression of of $\tilde{H}_D$, with the operators ${\hat{a}}$ and ${\hat{a}^{\dagger}}$ located on the off diagonal, suggests to look for solutions $\ket{\Psi_n}$ that decompose over number states $\ket{n}$ as
\begin{align}
\begin{pmatrix}
\lambda & {\hat{a}^{\dagger}} \\
{\hat{a}}         & -\lambda
\end{pmatrix}
\begin{pmatrix}
\alpha_{n+1} \ket{n+1}   \\  \alpha_{n} \ket{n}
\end{pmatrix}
= \tilde{E}_n
\begin{pmatrix}
\alpha_{n+1} \ket{n+1}   \\  \alpha_{n} \ket{n}
\end{pmatrix}
\label{eq:anzatz_dirac}
\end{align}
where the $\alpha_i$ are normalization coefficients and $\tilde{E}$ is the dimensionless eigenenergy normalized by $\mathcal{E}$. This anzatz yields the discrete spectrum
\begin{align}
\tilde{E}_n^{\pm} = \pm \sqrt{\lambda^2 + n +1} \qquad n\in \mathbb{N}
\end{align}
 shown in black in figure \ref{fig:Hermite}. It
 consists of two \textit{branches} $\pm$, separated by a gap around $\tilde{E}=0$. Each branch is constituted of an infinite number of discrete energy levels $n\in \mathbb{N}$, that correspond to localized modes around $x=0$.
The lower energy mode, $n=0$, decomposes into two branches $\tilde{E}_0^{\pm}=\pm \sqrt{k_y^2+1}$ that are associated to the eigenstate
\begin{align}
\Psi_0(x)=\braket{x|\Psi_0}=\begin{pmatrix}
\braket{x|1}\\
\braket{x|0}
\end{pmatrix}=
\begin{pmatrix}
\sqrt{2}x\\
1
\end{pmatrix} {\rm e}^{-\frac{x^2}{2}} \ .
\end{align}

Actually, the anzatz above does not account for all the solutions of \eqref{eq:H_Dop}. Indeed, an additional solution is suggested by the remarkable structure of $\ket{\Psi_n}$, whose components consist in shifted number states $\ket{n}$ and $\ket{n+1}$ with $n\in\mathbb{N}$.
Then, by extrapolating this structure, one can easily check that there exists another lower energy solution that one may somehow abusively call ''$n=-1$''.\footnote{This mode is actually called the ''zero-mode'' in the literature.} This solution reads
\begin{align}
\Psi_{-1}(x)=
\begin{pmatrix}
\braket{x|0}\\
0
\end{pmatrix}=
\begin{pmatrix}
1\\
0
\end{pmatrix} {\rm e}^{-\frac{x^2}{2}}
\label{eq:psi-1}
\end{align}
and satisfies
\begin{align}
\begin{pmatrix}
\lambda  &   {\hat{a}^{\dagger}}  \\
{\hat{a}}   & -\lambda
\end{pmatrix}
\begin{pmatrix}
 \ket{0}  \\  0
\end{pmatrix}
= \tilde{E}_{-1} \begin{pmatrix}
 \ket{0} \\ 0
\end{pmatrix}  \ .
\end{align}
 Unlike the previous solutions for  $n\in\mathbb{N}$, this solution yields a \textit{single} branch
\begin{align}
\tilde{E}_{-1} = \lambda \
\end{align}
that remarkably transits from the negative energy branch to the positive one when increasing $\lambda$. This is the simplest example of a \textit{spectral flow}, and in the following, $\lambda$ will be referred to as the spectral flow parameter.
Put differently, when $\lambda$ is swept from negative to positive values, the branch $-$ looses one state, while the branch $+$ gains one state. Writing $\Delta \mathcal{N}_\pm$ the net number of states which are gained by the branch $\pm$ when varying $\lambda$, one has $\Delta \mathcal{N}_\pm = \pm1$.

Let us finally notice that the mode that ''flows'' has several remarkable properties: It is (1) non-dispersive, (2) localized around $x=0$, that is where $m(x)$ changes sign, (3) fully polarized onto one component only, (4) the only accessible mode in the gap between the two branches $\pm$, (5) it propagates with a positive group velocity for any $k_y$. We will however see later that several of those properties are specific to this simple model. But for now, let us introduce the second example where the same spectral flow appears.

\subsection{Spectral flow of Landau levels in 3D Weyl semi-metals}

Another well-known condensed matter example where a spectral flow arises is that of 3D Weyl semi-metals subjected to a magnetic field \cite{VolovikBook,Yang2011,Spivak2013}. Such materials are characterized by pairs of 2-fold  band crossing points, around which the Bloch Hamiltonian can be expanded as massless Weyl Hamiltonians of the typical form \cite{Balents2011,Yang2011,Burkov2011}
\begin{align}
H_W = c
\begin{pmatrix}
{\hat{p}}_z  &  {\hat{p}}_x -  {\rm i} {\hat{p}}_y \\
{\hat{p}}_x + {\rm i} {\hat{p}}_y   &  -{\hat{p}}_z
\end{pmatrix} \ .
\end{align}
A magnetic field $\mathbf{B}=-B \mathbf{e}_z$,  is then applied in the $-z$ direction, and is accounted   in this Hamiltonian through the minimal coupling $\mathbf{{\hat{p}}}\rightarrow \hat{\boldsymbol{\pi}}=\mathbf{{\hat{p}}}+e \mathbf{A}(\mathbf{r})$ with $e$  the elementary electric charge and  where the electromagnetic vector potential reads $\mathbf{A}(\mathbf{r})=B/2(y, -x ,0)$ in the circular gauge. The generalized momentum $\mathbf{\pi}$ is gauge-invariant, unlike $\mathbf{p}$.

Since the system remains invariant by translation in the $z$ direction, one looks for plane wave solutions of wave number $k_z$ along that direction, which amounts to substitute  ${\hat{p}}_z \rightarrow \hbar k_z$ in the Hamiltonian. The situation is different in the $(x,y)$ plane: There, the commutator of the generalized momenta  $\hat{\pi}_x$ and $\hat{\pi}_y$ does not vanish, but verifies $[\hat{\pi}_x,\hat{\pi}_y]={\rm i} \hbar e B = {\rm i} \left( \hbar / \ell_B \right)^2$ where the characteristic length $\ell_B \equiv \sqrt{\hbar/eB}$ is called the magnetic length. This means that, unlike ${\hat{p}}_x$ and ${\hat{p}}_y$ that do commute, $\hat{\pi}_x$ and $\hat{\pi}_y$ are canonical conjugate variables. This allows us to introduce the annihilation and creation operators
\begin{align}
{\hat{a}}= \left(\frac{\ell_B}{\hbar}\right)\frac{\hat{\pi}_x+{\rm i} \hat{\pi}_y}{\sqrt{2}} \qquad {\hat{a}^{\dagger}} = \left(\frac{\ell_B}{\hbar}\right) \frac{\hat{\pi}_x-{\rm i} \hat{\pi}_y}{\sqrt{2}}
\end{align}
that we abusively also designate by ${\hat{a}}$ and ${\hat{a}^{\dagger}}$ because they similarly act on the number states as ${\hat{a}}\ket{n} = \sqrt{n}\ket{n}$, ${\hat{a}^{\dagger}}\ket{n} = \sqrt{n+1}\ket{n+1}$ and $[{\hat{a}},{\hat{a}^{\dagger}}]=1$.\footnote{Those ladder operators are actually related to the previous ones given in \eqref{eq:ann1} of the quantum harmonic oscillators in both directions as ${\hat{a}}={\hat{a}}_y-{\rm i} {\hat{a}}_x$, where ${\hat{a}}_x$ is given by \eqref{eq:ann1},  ${\hat{a}}_y$ is ${\hat{a}}_x$ after the substitution $x\rightarrow y$, and with $\ell = \ell_B$. This is thus essentially a 2D quantum harmonic oscillator in the $x$ $y$ plane \cite{CohenBook}. } The Weyl Hamiltonian under an axial constant magnetic field reads
\begin{align}
H_W =
\mathcal{E}
 \begin{pmatrix}
\lambda & {\hat{a}^{\dagger}} \\
{\hat{a}} & -\lambda
\end{pmatrix}
\label{eq:H_Wop}
\end{align}
where $\mathcal{E}\equiv \sqrt{2} c \hbar / \ell_B$ and $\lambda = k_z \ell_B / \sqrt{2}$. This Hamiltonian formally exactly the same as that of a 2D Dirac fermion with a mass that varies linearly along the $x$ direction \eqref{eq:H_Dop}, although the characteristic lengths of those two problems have different origin.
The energy spectrum of \eqref{eq:H_Wop} is therefore already shown in figure \ref{fig:Hermite}, and thus displays the same spectral flow as the  2D massive Dirac fermions.  Such a spectral flow of Landau levels is currently associated with the condensed matter signature of the chiral anomaly in 3D Weyl semimetals:  Assuming the Fermi energy lies around $E=0$, applying an electric field  along the $z$ direction generates a bulk current of electrons. Because Weyl points appear by pairs with opposite Chern numbers, they display opposite spectral flows. The chiral anomaly in that context refers to the fact that, although the number of particles is globally conserved when taking into accounts all the Weyl points, it is not conserved around each Weyl points (also called valleys) separetelly \cite{Zyuzin2012, Burkov2015, Landsteiner2016}.  \\

It is worth noticing that the ladder operators ${\hat{a}}$ and ${\hat{a}^{\dagger}}$  popped up in the 2D massive Dirac problem and in the 3D magnetic Weyl problem for very different reasons.  In the 3D Weyl problem, the key point is the coupling of the fermion charge to the magnetic field through the potential vector.  In contrast, in the 2D Dirac case, what matters is the existence of an anisotropic mass term that changes sign. This second mechanism can be found in a large variety of classical waves systems whose dynamics is governs by a set of equations of the form ${\rm i} \partial_t \psi = \hat{H}  \psi$ where $\hat{H}$ is an  Hermitian and linear operator, and  $\psi$ is the vector of the Fourier components of the classical  fields involved, such as pressure and velocity fields in fluids or electric and magnetic fields in optics. The  very spectral flow shown in figure \ref{fig:Hermite} was actually found in a classical analog of a quantum Chern insulator, namely in 2D gyromagnetic photonic crystals, where the role of the mass term was  played by a Faraday coupling \cite{Raghu2008}.

A last remark can be made about the discretization of the spectrum. This discretization can be expected in the 3D Weyl case with a magnetic field, as it corresponds to the emergence of Landau levels. For the 2D massive Dirac fermions, the discretization originates from the assumption of the linear spatial dependence of the mass, and would disappear if the shape of $m(x)$ becomes sharp (step function). But even in that case, the spectral flow between the two branches persists. In other words, it is robust to smooth changes of $m(x)$ provided it changes sign. It is thus natural to account for its robustness with topological tools.

\subsection{Topology of the symbol for spin $S=1/2$ fermions}

\subsubsection{Monopole - spectral flow correspondence}

The spectral flow $\Delta \mathcal{N}_\pm$ of  \textit{operators} such as $H_W$ and $\tilde{H}_D$, has indeed a topological counterpart in the associated \textit{symbol} of such operators, through topological indices called \textit{first Chern numbers} $\mathcal{C}_\pm$. The notions of symbol and of Chern numbers will be introduced in the next sections.  In the meantime, let us introduce the simple key relation
\begin{align}
\Delta \mathcal{N}_\pm = - \mathcal{C}_\pm
\label{eq:spectral_flow_simple}
\end{align}
that constitutes the monopole-spectral flow correspondence.
Such a relation is not obvious, but is certainly the simplest illustration of a very abstract theorem, called the Atiyah-Singer index theorem which, under some assumptions, essentially relates an \textit{analytical index} of an operator to a \textit{topological index} of its symbol \cite{AS1963, GiampieroBook, NakaharaBook}. This analytical index can then be related to a spectral flow of the operator, implying an equality of the type of equation \eqref{eq:spectral_flow_simple} \cite{APS1976}. In the next sections, we shall give an explicit and computable expression of these  indices for a quite general physical model beyond the spin-$1/2$ case discussed so far.

For now, let us keep focusing on the canonical $S=1/2$ example as described by the Hamiltonian \eqref{eq:H_Dop} (or \eqref{eq:H_Wop}), and remark that in the limit $|\lambda| \gg 1$, the diagonal term, that is proportional to the Pauli matrix $\sigma_z$, becomes dominant. The eigenmodes of the operator $H_W$ or $\tilde{H}_D$ then resemble that of  $\sigma_z$, and the typical energy between two consecutive number states $\ket{n}$ and $\ket{n+1}$  in the same branch, $+$ or $-$, tends to zero. There is then a clear separation of energy (or time) scales between the ''spin'' contribution, that one can associate to a fast motion, and the "oscillator" terms that imply  the operators ${\hat{a}}$ and ${\hat{a}^{\dagger}}$,  that one can thus associate to a slow motion.
With this in mind, it is natural to approximate the system by an ''adiabatic'' model, where the operators ${\hat{a}}$ and ${\hat{a}^{\dagger}}$ are treated classically. The Hamiltonian deduced from this procedure is precisely the symbol.

As we shall see, the symbol is, in our case, a $2\times 2$ matrix of parameters (the slow variables); the matrix structure being inherited from the two spin components $\ket{\pm}$ (fast variables) that are left untouched. The eigenmodes $\ket{\psi_\pm}$ of the symbol are thus vectors with two complex components, parametrized by the slow variables. These parametrized eigenstates are determined up to a phase, which gives rise to a structure called \textit{fiber-bundle} (see section \ref{sec:fiber}).
So, in the framework of the symbol, the space of fast motion is a fiber bundle over the phase space of slow variables (see e.g. \cite{Faure2001}). Fiber bundles are mathematical structures that generically possess a topology which accounts for how the parametrized eigenstates arrange and possibly ''twist''. In our case, those twists are encoded into the first Chern numbers $\mathcal{C}_\pm \in \mathbb{Z}$.

\subsubsection{Symbol Hamiltonian ${H_S}$}

The symbol-operator correspondence is routinely used by physicists when replacing by hand a wavenumber $k_x$ by the differential operator $-{\rm i} \partial_x$. Such a substitution is performed when the invariance by translation is lost in a the corresponding direction, for instance by adding a boundary condition.
Conversely, the replacement $-{\rm i} \partial_x \rightarrow k_x $ is performed when invariance by translation is assumed in a linear system governs by differential equations, so that a Fourier transform can be applied.
Such substitutions are routinely used by quantum physicists and were formalized by mathematicians  in order to make rigorous  the mapping between classical and quantum formalisms, i.e. between classical observables defined over phase space (e.g. the momentum $p$) and quantum observables defined by operators acting on a Hilbert space (e.g. $p\rightarrow \hat{p}=-{\rm i} \hbar \partial_x$). In this language, the classical observable $p$ is called the \textit{symbol} and the quantum observable $\hat{p}$ is called the \textit{operator} \cite{FollandBook, ZworskiBook}.

We will use this terminology and denote by ${\hat{\mathcal{H}}}_\text{op}$ the \textit{operator Hamiltonian} and by ${H_S}$ the \textit{symbol Hamiltonian}. Concretely, the previously discussed $H_W$ and $\tilde{H}_D$ in \eqref{eq:H_Wop} and \eqref{eq:H_Dop} constitute our operator Hamiltonian ${\hat{\mathcal{H}}}_\text{op}$.
Of course, the symbol Hamiltonian is itself also an (Hermitian) operator, since, being a matrix, it acts on vectors in $\mathbb{C}^2$. But its components are entities that commute. In contrast, ${\hat{\mathcal{H}}}_\text{op}$ acts on a Hilbert space that is a product $\mathscr{F}\otimes \mathscr{S}_{1/2}$ of an infinite dimensional Fock space $\mathscr{F}$ spanned by number states $\ket{n}$ with $n\in \mathbb{N}$, with  the two-dimensional spin $1/2$ Hilbert space $\mathscr{S}_{1/2}$. In other words, ${H_S}$ is a matrix of numbers, while ${\hat{\mathcal{H}}}_\text{op}$ is a matrix of operators, which makes its study much more involved.  Atiyah-Singer theorem provides a way to capture some properties of ${\hat{\mathcal{H}}}_\text{op}$ -- namely its spectral flow -- thanks to the topological property of a much simpler object, its symbol ${H_S}$.

Actually the mapping between the symbol and the operator is not unique. In particular, there are many ways to quantize a symbol. However, the symbol carries a topological information that does not depend on the choice of the quantization scheme (this part is called the \textit{principal} or \text{lead} symbol). We give a brief introduction to this formalism in the appendix \ref{sec:operator-symbol} (more advanced aspects can be found in various textbooks e.g. \cite{FollandBook, ZworskiBook}), and write the quite natural correspondence
\begin{align}
& \lambda_1 + {\rm i} \lambda_2   \quad  \leftrightarrow  \quad   {\hat{a}}    \\
&\lambda_1 - {\rm i} \lambda_2   \quad \leftrightarrow   \quad   {\hat{a}^{\dagger}}
\end{align}
where $\lambda_1$ and $\lambda_2$ are symbolic notations that designate two canonical conjugate observables in phase space.
In the 2D massive Dirac fermion problem \eqref{eq:H_Dop}, they stand for the classical position $\lambda_1=x$ and momentum $\lambda_2=p_x$ observables, while for the 3D Weyl fermion problem \eqref{eq:H_Wop}, they stand for the generalized momenta perpendicular to the applied magnetic field $\lambda_1=p_x+eA(y)$ and $\lambda_2=p_y+eA(x)$.

One can then come back to the  operator Hamiltonian of our Dirac/Weyl problems and write the correspondence between the symbol Hamiltonian ${H_S}$ and the operator Hamiltonian ${\hat{\mathcal{H}}}_\text{op}$
\begin{align}
{H_S} = \begin{pmatrix}
\lambda     & \lambda_1-{\rm i} \lambda_2 \\
\lambda_1 +{\rm i} \lambda_2  &    -\lambda
\end{pmatrix}
\quad
 \leftrightarrow
 \quad
{\hat{\mathcal{H}}}_\text{op} =
\begin{pmatrix}
\lambda     & {\hat{a}^{\dagger}} \\
{\hat{a}}  &    -\lambda
\end{pmatrix} \ .
\label{eq:symbol_op_correspondence}
\end{align}

Now that the symbol ${H_S}$ of ${\hat{\mathcal{H}}}_\text{op}$ has been introduced, one can discuss its topological properties.

\subsubsection{Degeneracy point of the symbol as a topological defect in 3D parameter space}
\label{sec:topo_defect}

The purpose of this section is to make concrete what is meant by the \textit{topology} of the symbol ${H_S}$.
Let us denote by $\Lambda=\mathbb{R}^3 \backslash \{ {\boldsymbol{\lambda}}^\star \}$ the \textit{parameter space}, where  ${\boldsymbol{\lambda}}^\star$ designates the degeneracy point(s) of ${H_S}$, and ${\boldsymbol{\lambda}}=(\lambda_x,\lambda_y,\lambda_z)=(\lambda_1,\lambda_2,\lambda)$ a point in $\Lambda$.
The symbol Hamiltonian \eqref{eq:symbol_op_correspondence} we are interested in takes the form
\begin{align}
{H_S} = {\boldsymbol{\lambda}} . \boldsymbol{\sigma}
\label{eq:2x2simple}
\end{align}
with $ \boldsymbol{\sigma} = (\sigma_x,\sigma_y,\sigma_z)$ the vector of Pauli matrices.
Such a $2\times 2$ Hermitian matrix constitutes the ubiquitous minimal model in condensed matter physics that owns a topological property. As we shall detail in the following, this topological property can be apprehended in two different ways.

Let us start by pointing out that the eigenvalues $E_\pm=\pm |{\boldsymbol{\lambda}}|$  of $H$ are degenerated at a point ${\boldsymbol{\lambda}}^\star=0$.  This is actually a generic property of Hermitian matrices, as  found by Wigner and von-Neumann: Indeed, the dimension of the space of Hermitian matrices with a two-fold degeneracy is $n^2-3$ (with $n$ the dimension of the matrix), while it is $n^2$ for generic Hermitian matrices. In other words, the codimension of two-fold degenerated Hermitian matrices is $3$. This implies that, generically, i.e. in the absence of other constraints imposed by symmetries, $3$ parameters must be varied to find a degeneracy between two eigen-energies of a Hermitian matrix. This apparently trivial property actually turns out to play a key role in what follows.

In fact, a degeneracy point can be see as a "defect" in parameter space. Similarly to actual topological defects of vector fields in real space (e.g. vortices, dislocations...) one can characterize the degeneracy point by the homotopy property of the map ${\boldsymbol{\lambda}} : \mathcal{S}_{\boldsymbol{\lambda}} \rightarrow \mathcal{S} \subset \mathbb{R}^3$ where $\mathcal{S}_{\boldsymbol{\lambda}}$ is any surface that encloses the degeneracy point. Such maps, between two closed orientable manifolds, (here two surfaces), are classified by an integer-valued number, called the \textit{degree}. The degree of a map is an homotopy invariant, that is, it cannot change value under a continuous deformation of the map. In particular, the two surfaces  $\mathcal{S}_{\boldsymbol{\lambda}} $ and  $\mathcal{S} $ can conveniently be continuously deformed into spheres $S^2$ in $\mathbb{R}^3$. More generally, the homotopy properties of maps from a sphere $S^n$ to a sphere $S^m$ are classified by homotopy groups noted $\pi_n(S^m)$. The search for homotopy groups of spheres is in general a difficult open problem in mathematics. However, there is a particular situation where the result is known, which is $\pi_n(S^n) = \mathbb{Z}$, and thus is particular $\pi_2(S^2)= \mathbb{Z}$. This means that point-like topological defects of 3D vector fields are classified by integers. These integers are precisely the degree, which, geometrically, counts how many times the image of $\mathcal{S}_{\boldsymbol{\lambda}}$ by the map ${\boldsymbol{\lambda}} \in \mathcal{S}_{\boldsymbol{\lambda}} \rightarrow h({\boldsymbol{\lambda}}) \in \mathcal{S} $ covers  $\mathcal{S}$, or equivalently how many times the normalized map $\mathcal{S}_{\boldsymbol{\lambda}} \rightarrow \mathbf{n} \equiv h/\mathbf{h} \in S^2$ wraps the unit sphere. In our simple example, where $h$ is identity, the degree is $1$. In the following, we shall discuss much richer situations where the degree can take arbitrary values.
The vector field around topological defects characterized with different degrees cannot be smoothly deformed into each other. This characterizes the topological robustness of the defect.

\subsubsection{Phase singularity of the eigenstates}
\label{sec:fiber}
Characterizing the degeneracy point of the symbol Hamiltonian with the degree of a map amounts to ignoring the spinorial structure of the problem, carried by the Pauli matrices in \eqref{eq:2x2simple}. Indeed, $H_S$ is a matrix and a subtle topological property is hidden in its normalized eigenvectors
 $\psi_\pm$ %
\begin{align}
\ket{\psi_+} &=
\frac{ {\rm e}^{{\rm i} \chi_+}}{\sqrt{2 |{\boldsymbol{\lambda}}|(|{\boldsymbol{\lambda}}|-\lambda_z)}}
 \begin{pmatrix}
\lambda_x-{\rm i} \lambda_y \\
|{\boldsymbol{\lambda}}|-\lambda_z
 \end{pmatrix}
\\
\ket{ \psi_-} &=
\frac{ {\rm e}^{{\rm i} \chi_-}}{\sqrt{2 |{\boldsymbol{\lambda}}|(|{\boldsymbol{\lambda}}|+\lambda_z)}}
 \begin{pmatrix}
\lambda_x+{\rm i} \lambda_y \\
 -|{\boldsymbol{\lambda}}|-\lambda_z
 \end{pmatrix}
 \label{eq:state2x2_1}
\end{align}
where  $\chi_\pm$ is an arbitrary phase that depends on ${\boldsymbol{\lambda}}$. For any value of ${\boldsymbol{\lambda}}$, there is a gauge freedom for the choice of $\chi_\pm$, as usual in quantum mechanics. This means that for a given ${\boldsymbol{\lambda}}$, there is not a single state that describes the system, but a continuous family of them, whose members that all differ from each other by a phase factor. This is, in a sense, a multi-valuation of the solution. In contrast, the projectors  $\Pi_\pm \equiv \ket{\psi_\pm}\bra{\psi_\pm}$ do not depend on the gauge and are thus always single-valued. Once the gauge is chosen, the eigenstate should be single-valued as well. This is however not always the case, as one can suspect from \eqref{eq:state2x2_1} when ${\boldsymbol{\lambda}}$ is aligned along the $z$ direction, i.e. for $\pm \lambda_z = \pm |{\boldsymbol{\lambda}}|$.

To see it more clearly, let us use the spherical coordinates, and parametrize the spinor on the Bloch sphere, where it becomes explicit that the eigenstates  only depend on $\mathbf{n}$ (up to the gauge choice)~:
\begin{align}
\ket{\psi_+} =
{\rm e}^{{\rm i} \chi_+}
 \begin{pmatrix}
\cos \frac{\theta}{2} {\rm e}^{-{\rm i} \phi}\\
\sin \frac{\theta}{2}
 \end{pmatrix}
 \qquad \text{and} \qquad
 \ket{\psi_-} =
{\rm e}^{{\rm i} \chi_-}
 \begin{pmatrix}
\sin \frac{\theta}{2} {\rm e}^{-{\rm i} \phi}\\
-\cos \frac{\theta}{2}
 \end{pmatrix} \ .
 \label{eq:state2x2_2}
\end{align}
There, setting $\mathbf{n} = n_z >0$ or equivalently $\theta=0$ (north pole), the eigenstate $\Psi_+$ becomes
\begin{align}
\ket{\psi_+} \overset{\theta=0}{=}  {\rm e}^{{\rm i} \chi_+}
\begin{pmatrix}
{\rm e}^{-{\rm i} \phi}\\
0
 \end{pmatrix}
 \qquad \text{(at the north pole of the Bloch sphere)}
\end{align}
 which is multi-valued even after the gauge has been fixed, unless we make the particular gauge choice $ {\rm e}^{{\rm i} \chi_+} = {\rm e}^{{\rm i} \phi}$.  This fixes the multi-valued of the eigenstate, but only at the north pole. Indeed, in this gauge, the eigenstate becomes multivalued at the south pole $\theta=\pi$
\begin{align}
\ket{\psi_+} \overset{\theta=\pi}{=}
\begin{pmatrix}
0\\
{\rm e}^{{\rm i} \phi}
\end{pmatrix} \qquad (\text{with}\  {\rm e}^{{\rm i} \chi_+} = {\rm e}^{{\rm i} \phi})  \qquad \text{(at the south pole of the Bloch sphere)} \ .
\end{align}
Actually, the multi-valuedness of the eigenstates cannot be fixed by a \textit{global} gauge choice. Or say otherwise, there does not exist a smooth function $\chi$ that makes the eigenstate single-valued  everywhere on the Bloch sphere. One can only choose \textit{locally} a gauge that makes the eigenstate single-valued, namely
\begin{align}
\label{eq:gauge_S}
&\text{south gauge:}\quad  \chi=0 \qquad  \ket{\psi_+^{S}} =
\begin{pmatrix}
\cos \frac{\theta}{2} \, {\rm e}^{-{\rm i} \phi} \\
\sin \frac{\theta}{2}
\end{pmatrix}
&\qquad \text{for}\  \theta \neq 0 \\
&
\text{north gauge:}\quad  \chi=\phi \qquad  \ket{\psi_+^{N}} =
\begin{pmatrix}
\cos \frac{\theta}{2}  \\
\sin \frac{\theta}{2} \, {\rm e}^{{\rm i} \phi}
\end{pmatrix}
&\qquad \text{for}\  \theta \neq \pi
\label{eq:gauge_N}
\end{align}
 Thus,  the eigenstates are only piece-wise single-valued, and related by a gauge transformation
 \begin{align}
 \ket{\psi_+^N} = {\rm e}^{{\rm i} \phi} \ket{\psi_+^S} \ .
 \label{eq:gaugeSN}
 \end{align}
 This reflects a topological property of the eigenstates, as their phase cannot be defined smoothly globally: it has to have a singularity somewhere, similarly to a vortex. This topological property is inherently related to the gauge freedom, and is mathematically described by the theory of fiber-bundles \cite{NakaharaBook, GeometricPhaseBook}.

Fiber bundles constitute a rigorous construction of a product of spaces. Locally, the product of two spaces, say a segment $[-1,1]$ with a circle $S^1$, is given by the usual cartesian product. But globally, this does not necessarily hold since a twist may occur. In that case, the fiber-bundle is said to be topological. This is the celebrated example of the Moebius strip that differs by a twist from the cylinder. These two structures both result from the product of the segment with the circle, but the Mobius strip can only be seen as the cartesian product $[-1,1]\times S^1$  locally. Fiber bundles are in particular useful to make meaningful the notion of parametrized vector space. The vector space (here the segment) is called the fiber, while the space that parametrizes the fiber (here the circle) is called the base space.

Another intuitive example of a fiber-bundle is the family of tangent vector spaces to the sphere $S^2$. In that case, the vector space is the tangent plane to the sphere and the base space is the sphere itself. The tangent vector bundle is thus the collection of all tangent vectors at any point of the sphere. This vector bundle is also topological, and its twist can be seen in the tangent vector field to the sphere, that necessarily has two vortices, which are precisely points where the phase (i.e. the angle of the tangent vectors with respect to an arbitrary direction) is multi-valued. The position of these singularities depends on the gauge choice, in the language of physicists, or of the section of the fiber bundle as mathematicians would say. But the total vorticity, equals to $2$ in that example, cannot change. This is a topological invariant of this fiber bundle, and this result is known as the hairy ball theorem.

The fiber bundles we are interested in resemble this hairy ball~: the base space is also the sphere $S^2$ but the fibers (the hairs) are complex vector spaces made of all the eigenstates $\psi$ that only differ by a phase (or gauge choice), that is
\begin{align}
F_\pm({\boldsymbol{\lambda}}) = \{ {\rm e}^{{\rm i} \chi_\pm({\boldsymbol{\lambda}})}  \psi_\pm({\boldsymbol{\lambda}}), {\rm e}^{{\rm i} \chi_\pm({\boldsymbol{\lambda}})}  \in U(1) \} \ .
\label{eq:fiber}
\end{align}
It is said that the fibers $F_\pm({\boldsymbol{\lambda}})$  define the \textit{equivalent class} of the states $\psi_\pm({\boldsymbol{\lambda}})$
as  they have the same projector and thus describe the same physical state. For each mode $\pm$, the fiber bundle (of eigenvector bundle) of interest consists in the continuous collection of such $U(1)$-fibers $F_+$ or $F_-$, where the base space is a sphere that encloses the degeneracy point. A two-fold degeneracy point is thus associated with two eigenvector bundles.
There are many other fiber bundles, of different dimensions, of real or complex nature, with group structures different from $U(1)$, but this simple example already englobes many interesting physical situations.

  \begin{figure}[h!]
  \centering
\includegraphics[width=14cm]{fig_bundle}
\caption{\footnotesize{(a) $U(1)$ fiber-bundle as a continuous collection of fibers $F$ over the sphere defined in parameter space ${\boldsymbol{\lambda}}$. Two states in the same fiber only differ by a phase; they thus have the same projector onto the base space. The impossibility to define smoothly this phase over the base space is a manifestation of the topological property of the fiber bundle, encoded into the Chern number. (b) This Chern number is the quantized flux through the sphere of the Berry curvature generated by the degeneracy point, called Berry monopole. The existence of a non-zero Chern number imposes the 1-form Berry connection to be  defined smoothly only piecewise. }}
\label{fig:bundle}
\end{figure}

\subsubsection{Chern number from the Berry curvature as an obstruction to Stokes theorem}

To express the topological invariants of fiber bundles, one needs now to introduce some definitions of differential geometry.
The language of differential forms is particularly convenient to reveal the topological and geometrical structures of physical systems. Of main importance to characterize the eigenstates $\ket{\psi_n({\boldsymbol{\lambda}})}$ of a quantum Hamiltonian that is smoothly parametrized in $\Lambda$, or actually any eigenstates of a parametrized linear Hermitian eigenvalues problem of the form
\begin{align}
H({\boldsymbol{\lambda}}) \ket{\psi_n({\boldsymbol{\lambda}})} = E_n({\boldsymbol{\lambda}}) \ket{\psi_n({\boldsymbol{\lambda}})}
\label{eq:schro_param}
\end{align}
 that one can encounter in particular in wave physics, are geometrical tools introduced by Berry and Simon, called \textit{Berry connection} and \textit{Berry curvature} \cite{Berry1984, Simon}.  The Berry connection is a 1-form defined from each parametrized eigenstate as\footnote{This definition fits the one of Berry and Simon but differs by a minus sign from the one mostly encountered in solid states physics when dealing with Bloch bands.}
 %
 %
\begin{equation}
\begin{aligned}
\mathcal{A}^{(n)}({\boldsymbol{\lambda}}) \equiv&\,  {\rm i}  \bra{\psi_n} {\rm d} \ket{\psi_n} \\
=&\,  {\rm i}  \bra{\psi_n} \frac{\partial}{\partial \lambda_{j}} \ket{\psi_n} {\rm d} \lambda_j\\
\equiv& \, A_j^{(n)}({\boldsymbol{\lambda}}) {\rm d} \lambda_j  \ .
\label{eq:connection_def}
\end{aligned}
\end{equation}
where the sum over the $j$ index is implicit.\footnote{We keep the upper index to denote the band index.}
 The one-form Berry connection is not gauge invariant~: when a $U(1)$-gauge transformation is performed on a state as
 \begin{equation}
 \begin{aligned}
 \label{eq:gauge_transform}
 \ket{\psi} \rightarrow   {\rm e}^{{\rm i}\chi({\boldsymbol{\lambda}})} \ket{\psi}  = \ket{\tilde{\psi}}
 \end{aligned}
 \end{equation}
 then the Berry connection transforms as
\begin{align}
\mathcal{A}^{(n)} \rightarrow \tilde{\mathcal{A}}^{(n)} &=\, {\rm i} \bra{\tilde{\psi}_n} {\rm d} \ket{\tilde{\psi}_n} \\
& =\, -{\rm d} \chi + \mathcal{A}^{(n)} \ .
\label{eq:connection_gauge}
\end{align}
This is very analogous to a vector potential in electromagnetism, where moreover the language of differential forms automatically allows us to work in arbitrary dimension, where the notion of e.g. the \textit{curl} becomes obscure. Similarly to electromagnetism, one can introduce the equivalent of a magnetic field, called the Berry curvature, as
\begin{align}
 \mathcal{F}^{(n)}({\boldsymbol{\lambda}}) & \equiv \, {\rm d} \mathcal{A}^{(n)}({\boldsymbol{\lambda}}) \\
 &= \frac{\partial A^{(n)}_j}{\partial \lambda_i} {\rm d} \lambda_i \wedge {\rm d} \lambda_j
 \end{align}
 where the implicit sum now runs over $i$ and $j$. Because of the anti-symmetry ${\rm d} \lambda_i \wedge {\rm d} \lambda_j = - {\rm d} \lambda_j \wedge {\rm d} \lambda_i$ of the wedge product $\wedge$, one has
 \begin{align}
 \label{eq:coeff_Berry}
 \mathcal{F}^{(n)}({\boldsymbol{\lambda}}) &= \frac{1}{2} \left(
 \frac{\partial }{\partial \lambda_i} A^{(n)}_j - \frac{\partial }{\partial \lambda_j} A^{(n)}_i
    \right) {\rm d} \lambda_i \wedge {\rm d} \lambda_j \\
    &\equiv \frac{1}{2} F^{(n)}_{ij}\,   {\rm d} \lambda_i \wedge \lambda_j
 \end{align}
 where the second equality defines the coefficients $F^{(n)}_{ij} ({\boldsymbol{\lambda}})$ of the 2-form Berry curvature $ \mathcal{F}^{(n)}$. Those  are left invariant under the $U(1)$-gauge transformation \eqref{eq:gauge_transform}.
Consistently, they  are observable quantities that manifest themselves for instance in the semi-classical equations of motion of wavepackets travelling in spatially slowly varying media, and in the adiabatic evolution of quantum states \cite{Sundaram1999, Xiao_Niu2010, Jotzu2014, Wimmer2017, Luu2018, Cho2018}.

Mathematically, this gauge invariance of the Berry curvature is straightforwardly inferred from the formalism of differential forms. Indeed, at this stage, we have introduced the Berry curvature as a 2-form that derives from of the 1-form Berry connection by taking the exterior derivative ${\rm d}$. Such a construction is called an \textit{exact} form (see appendix B), and it is known as a theorem that the exterior derivative of an exact form vanishes. Therefore, under the gauge transformation \eqref{eq:gauge_transform}, the Berry curvature transforms as $\mathcal{F} \rightarrow \mathcal{F} + {\rm d} \mathcal{F} = \mathcal{F}$.

 A very subtle and crucial point is that the Berry curvature may actually not be  an exact form for every ${\boldsymbol{\lambda}} \in \Lambda$.
The topology arises precisely in such a situation, as the obstruction to define  \textit{globally}  the Berry curvature as $\mathcal{F}({\boldsymbol{\lambda}})= {\rm d} \mathcal{A}({\boldsymbol{\lambda}})$ for every ${\boldsymbol{\lambda}}$. In that case, this equality is only valid  \textit{locally}, i.e. on sub-domains of $\Lambda$, and    $\mathcal{A}^{(n)}({\boldsymbol{\lambda}})$, (and consequently ${\rm d} \mathcal{A}^{(n)}({\boldsymbol{\lambda}})$), is only piecewise well-defined.
This is precisely the case when the phase of the eigenstates, such as  $\psi_\pm$  (see \eqref{eq:state2x2_2}) of the symbol Hamiltonian ${H_S}$ for the spin $S=1/2$, cannot be smoothly defined over $\Lambda$.

As we shall see, the obstruction to define globally the phase of the eigenstates, or equivalently their Berry curvature, is captured quantitatively by an integer-valued topological number, called the  \textit{first Chern number}, expressed as
\begin{align}
\mathcal{C}_n \equiv \frac{1}{2\pi} \int_{S^2} \mathcal{F}^{(n)}  \in \mathbb{Z} \ .
\label{eq:cc_def}
\end{align}
Again, there is a strong analogy with electromagnetism, where the expression \eqref{eq:cc_def} of the Chern number is given by the flux of the Berry curvature. This analogy can be pushed further by expressing the coefficients $F^{(n)}_{ij}$ of the Berry curvature as
\begin{align}
\label{eq:practical_Berry}
F^{(n)}_{ij} ({\boldsymbol{\lambda}})= {\rm i} \sum_{m\neq n} \frac{ \bra{\psi_n} \partial_{\lambda_i} H\ket{\psi_m} \bra{\psi_m} \partial_{\lambda_j} H \ket{\psi_m} -  \bra{\psi_n} \partial_{\lambda_j} H\ket{\psi_m} \bra{\psi_m} \partial_{\lambda_i} H\ket{\psi_m}}{(E_n({\boldsymbol{\lambda}}) - E_m({\boldsymbol{\lambda}}))^2}
\end{align}
which is obtained by differentiating \eqref{eq:schro_param} and inserting it into \eqref{eq:coeff_Berry}. This formula has a practical interest, and we shall use it in a next section to compute analytically the Berry curvature for a general spin-$S$ model. For now, let us stress that the denominator of \eqref{eq:practical_Berry} reveals that the amplitude of the Berry curvature of a band $n$ increases as this band gets closer to the other bands $m$ in parameter space $\Lambda$. In particular, the Berry curvature diverges at a degeneracy point (provided the numerator is regular). For that reason, degeneracy points (or band crossing points) are play the role of point-like source of the Berry curvature. The Chern number number \eqref{eq:cc_def} can then be understood as the (topological) charge of this monopole, with the (close) surface of integration $S^2$  enclosing the degeneracy point $\lambda^\star$  in $\Lambda$ space. Analogies with the  magnetic monopole are naturally often made, since the quantization of the electric charge was found by Dirac to be a consequence of the obstruction to smoothly define the vector potential that couples to the quantum wavefunction around such monopoles. Note however that a $N$-fold band crossing points -- that we dub here Berry-Chern monopoles -- imply $N$ bands, and thus $N$ fiber bundles, each of them being characterized with a Chern number $\mathcal{C}_n$. Such degeneracy points are thus associated with $N$ topological charges rather than a single one. In the simple case of two-band crossing points, which are the most often encountered in the literature,  such as Weyl points, this terminology of topological charge is currently used without too much ambiguity, because the two Chern numbers are necessarily opposite. This follows from the important property of the Berry curvature that satisfies
\begin{align}
\sum_n F_{ij}^{(n)}({\boldsymbol{\lambda}}) =0
\end{align}
where the sum runs over all the bands $n$. This property, that can be derived from the expression \eqref{eq:practical_Berry}, yields the important relation on the Chern numbers
\begin{align}
\label{eq:chern_constraint}
\sum_n \mathcal{C}_n =0 \ .
\end{align}
Thus, an $N$-fold degeneracy point ${\boldsymbol{\lambda}}^\star \in \Lambda$ is characterized by $N$ Chern numbers $\mathcal{C}_n$ satisfying this constraint. We refer here to such points as Berry-Chern monopoles.\\

To see how the Chern number \eqref{eq:cc_def} captures the obstruction to smoothly define the eigenstates (or to have the Berry curvature as an exact $2$-form) globally, let us apply Stokes theorem \eqref{eq:stokes}
\begin{align}
\int_{S^2}   \mathcal{F} =\int_{S^2} {\rm d}  \mathcal{A}     =  \int_{\partial S^2}    \mathcal{A}
\label{eq:paradox}
\end{align}
(where we drop the band index when unecessary).
Since the sphere $S^2$ is a closed surface, it does not have a boundary (i.e. $\partial S^2=0$), so that the right hand side member of this equality always vanishes. This result is only valid if $ \mathcal{F} =  {\rm d} \mathcal{A}$ everywhere on the domain of integration. If, on the other hand, the equality  $\mathcal{F} = {\rm d} \mathcal{A}$  only holds piecewise,  then Stokes theorem cannot be used over $S^2$. For that reason, it is said that the Chern number is an obstruction to Stokes theorem. Let us illustrate this point with the eigenstate $\psi_+$ of the spin $1/2$ problem. We saw in Eqs. \eqref{eq:gauge_S} and \eqref{eq:gauge_N} that such a state cannot can only be  smoothly defined piecewise, in different gauges, and therefore so is the Berry connections
\begin{align}
\mathcal{A} = {\rm i} \bra{\psi} \partial_\theta \ket{\psi}{\rm d} \theta  + {\rm i} \bra{\psi} \partial_\phi \ket{\psi} {\rm d} \phi
\end{align}
leading to
\begin{align}
&\text{south gauge:}\quad  \chi=0 \qquad \mathcal{A}^{(+)}_S = \cos^2\frac{\theta}{2}{\rm d} \phi &\qquad \text{for}\  \theta \neq 0 \\
&\text{north gauge:}\quad  \chi=\pi \qquad \mathcal{A}^{(+)}_N = -\sin^2\frac{\theta}{2}{\rm d} \phi &\qquad \text{for}\  \ \theta \neq \pi \, \
\end{align}
The introduction of the Berry curvature as the derivative of the Berry connections $\mathcal{F} = {\rm d} \mathcal{A}$ must accordingly be performed over the corresponding domains, associated with the two different gauge choices.
The Chern number can then be computed via Stokes theorem by cutting the sphere on different patches where this equality makes sense, for instance over north hemisphere (NH) and south hemisphere (SH) respectively, as illustrated in figure \ref{fig:bundle}~:
\begin{align}
\mathcal{C}_+ &= \frac{1}{2\pi} \int_{S^2} \mathcal{F}_+\\
&=  \frac{1}{2\pi} \int_{\text{NH}}\mathcal{F}_+  +\frac{1}{2\pi} \int_{\text{SH}}\mathcal{F}_+ \\
&=  \frac{1}{2\pi} \int_{\text{NH}}{\rm d} \mathcal{A}_+^N  +\frac{1}{2\pi} \int_{\text{SH}}{\rm d} \mathcal{A}_+^S
\end{align}
so that Stokes theorem can be safely used to get
\begin{align}
\mathcal{C}_+ =   \frac{1}{2\pi} \int_{\partial \text{NH}}\mathcal{A}_+^N  +\frac{1}{2\pi} \int_{\partial \text{SH}}\mathcal{A}_+^S  \ .
\end{align}
where the boundaries $\partial \text{NH}$ and $\partial \text{SH}$ are any common path $\gamma$ that does not cross a pole, and are oriented in opposite directions for the two domains, so that
\begin{align}
\mathcal{C}_+ =   \frac{1}{2\pi} \int_{\gamma}\mathcal{A}_+^N  - \mathcal{A}_+^S  \ .
\label{eq:chern_gauge}
\end{align}
The two connections involved are defined on domains where the eigenstates $\ket{\psi_+^N}$ and $\ket{\psi_+^S}$ are related by the gauge transformation \eqref{eq:gaugeSN}. According to \eqref{eq:connection_gauge}, the Berry connections $\mathcal{A}_+^N$  and $\mathcal{A}_+^S$ inherit this relation as
\begin{align}
\mathcal{A}_+^N  = \mathcal{A}_+^S - {\rm d} \phi
\end{align}
 the Chern number is directly inferred as
\begin{align}
\mathcal{C}_+= \frac{1}{2\pi}\int_0^{2\pi} -{\rm d} \phi
=-1 \ .
\end{align}
Conversely, one finds $\mathcal{C}_-=+1$ if $\psi_-$ is considered instead of $\psi_+$, consistently with the constraint \eqref{eq:chern_constraint}.
The Chern numbers $\mathcal{C}_\pm$ of the symbol Hamiltonians $H_S$ thus satisfy the monopole-spectral flow correspondence \eqref{eq:spectral_flow_simple} for the spin $1/2$ band crossing problem.\\

Let us conclude this section by a few remarks on the computation of the Chern numbers. First one can check the value of the Chern number by inserting the expressions of the Berry connections $\mathcal{A}^{(+)}_S$ and $\mathcal{A}^{(+)}_N$ in \eqref{eq:chern_gauge}  and then taking  the path $\gamma$ as the equator of $S^2$, that is  $\theta=\pi/2$. In passing, one can also check that these two different expressions for the Berry connection indeed yield the same Berry curvature $\mathcal{F}^{+}=-\frac{1}{2}\sin{\theta} {\rm d} \theta \wedge {\rm d} \phi$ whose expression is valid for any $\theta$. A direct integration of this quantity over $S^2$ in spherical coordinates also consistently yields $\mathcal{C}_+ = -1$.
However,  let us stress that the calculation  of the Chern number presented above actually did not require the explicit computation of the Berry connection, neither that of the Berry curvature. It only lies on the gauge transformation between the different patches where the eigenstates are well-defined. This can be slightly rewritten with
\begin{align}
&\ket{\psi} \rightarrow \ket{\tilde{\psi}}  = {\rm e}^{{\rm i}\chi({\boldsymbol{\lambda}})} \ket{\psi} \\
&\mathcal{A} \rightarrow \mathcal{A} + {\rm i} {\rm e}^{-{\rm i} \chi} {\rm d} {\rm e}^{{\rm i} \chi}
\end{align}
where ${\rm e}^{{\rm i} \chi} \in U(1)$ is called the \textit{transition function}, and the Chern numbers then correspond to the winding numbers of this function
\begin{align}
\mathcal{C}  = \frac{1}{2\pi} \int_\gamma \mathcal{A} -  \tilde{\mathcal{A}}
            = -\frac{1}{2\pi {\rm i}}  \int_\gamma {\rm e}^{-{\rm i} \chi} {\rm d} {\rm e}^{{\rm i} \chi} \ .
            \label{eq:transition_function}
\end{align}

\subsection{Spectral flow from the Berry-Chern monopole of the symbol : a simple methodology}

\subsubsection{Predicting a spectral flow from 2D and 3D dispersion relations}

At this stage, we have verified the monopole - spectral flow correspondence  \eqref{eq:spectral_flow_simple} for spin $1/2$ models with a linear band crossing point. It is worth stressing that the correspondence \eqref{eq:symbol_op_correspondence}  between the symbol and the operator  allows one to predict a spectral flow, that is the relevant physical information, from the symbol Hamiltonian which is a much simpler object to manipulate. This provides us with a powerful tool, especially because in the cases we address here, the symbol's spectrum is nothing but the dispersion relation of the homogeneous system. A generic method to predict a spectral flow therefore consists in looking for a band crossing in the dispersion relation that appears as point in a 3D parameter space. This is actually the standard situation for two-band crossings in the absence of symmetry, according to Wigner-von-Neumann theorem.
The computation of the Chern number, detailed in the previous section, may seem technical to non familiar readers. But, since generic band crossings only involve two bands, the effective Hamiltonian describing locally such crossings is the spin $1/2$ Hamiltonian we have discussed, up to unitary transformations, so that the Chern numbers can only be $+1$ or $-1$.  Therefore, it is sufficient to find out a band crossing point with a linear dispersion relation (which is an easy task) in order to predict a spectral flow of one mode (which is much more involved to find in full generality). If the band-crossing point involves more bands and/or a nonlinear dispersion relation, then a careful treatment of the Chern numbers must be done. This is the purpose of the section \ref{sec:general}.

To guarantee the emergence of a spectral flow, one needs to turn the symbol problem, that describes an homogeneous system,  into an inhomogeneous one described by the operator problem. The examples discussed in the first section provides two physical recipes to perform such a transformation, depending on the dimension of the system. Those two methods remain valid even beyond two-fold band crossing points with a linear dispersion relation.

If, one the one hand, the physical problem is three-dimensional, (i.e. the dispersion relation $E({\mathbf{k}})$ is such that ${\mathbf{k}}$ is three-dimensional), then one must look for a physical mechanism that turns two of the three components of the wave number ${\mathbf{k}}$ into canonical conjugate variables. This is what happens when one uses the minimal coupling to account for the application of a magnetic field in Weyl semi-metals. The spectral flow then consists in the transfer of the lowest Landau level leading to a bulk propagating mode in the direction parallel to the magnetic field. This propagating mode leads to a current when an electric field is applied (provided the mode crosses the Fermi energy), which is interpreted as a signature of the chiral anomaly.

If, on the other hand, the physical system is two-dimensional, (i.e. the dispersion relation $E({\mathbf{k}})$ is such that ${\mathbf{k}}$ is two-dimensional), then one must look for a third parameter to be tuned  such that it changes sign when a gap between wavebands closes. This is the role of the mass term in the 2D massive Dirac problem. Note that the symmetries are a good guide to reveal such a parameter: Indeed, degeneracy points generically occur when a symmetry is exceptionally restored in the problem, such as inversion or time reversal. The parameter playing the role of the mass term must thus break such a symmetry. The operator Hamiltonian is then obtained when considering that the mass term becomes a smooth function of space $m(x)$, and a spectral flow is expected along the perpendicular direction $k_y$ if $m(x_0)$ vanishes at some $x_0$. This necessary vanishing of $m(x)$ in the inhomogeneous problem, is the dual (i.e. operator) counter-part of the vanishing of the parameter $m$ in the homogeneous  problem, which is at the source of the Berry-Chern monopole.
The physical interpretation of that spectral flow is that of an uni-directional interface mode, i.e. a wave confined in the perpendicular direction ($x$) to the line $m(x_0)=0$ but that propagates along the interface ($y$) in one direction only.

\subsubsection{Bulk-interface correspondence in 2D}

The trapped unidirectional modes of the spectral flow in 2D inhomogeneous systems look very much like the topological chiral edge states propagating along the boundary of two-dimensional quantum Hall like systems, i.e. Chern insulators.

Let us recall that Chern insulators are band insulators whose energy (or frequency)  Bloch bands  $n$ are characterized by a Chern number $C_n$. Those \textit{band} Chern numbers are given by the integral of the Berry curvature of each Bloch bundle (the fiber bundle of Bloch eigenstates) over the 2D Brillouin zone. The periodicity in both directions of the 2D Brillouin zone (BZ) making it equivalent to a 2-torus $T^2$ torus, the band Chern number reads
 \begin{align}
 \label{eq:band_Chern}
 C_n=\frac{1}{2\pi}\int_{BZ\sim T^2}\mathcal{F}^{(n)}({\mathbf{k}}) \ .
 \end{align}
 In contrast, the surface of integration of the Berry curvature we considered for the monopole Chern numbers $\mathcal{C}_n$ encloses a degeneracy point in $\mathbb{R}^3\backslash \{{\boldsymbol{\lambda}}^\star\}$ space.

The band Chern number \eqref{eq:band_Chern} characterizes the obstruction to smoothly define the phase of the Bloch eigenstates over the Brillouin zone. Time-reversal symmetry must be broken for the band Chern number to be non-zero. A non-zero band Chern number of a (homogeneous) periodic system,  guarantees the existence of a unidirectional (chiral) edge mode that manifests itself as a spectral flow in the non-periodic system with boundaries. This is the celebrated bulk-edge correspondence \cite{HatsugaiPRL1993, HatsugaiPRB1993, GrafPorta2013}.
There is a close connection between the bulk-edge correspondence and the monopole - spectral flow correspondence for 2D systems, since  an interface between two distinct Chern insulators can be described by a single Hamiltonian operator with a varying mass term whose sign controls the topological transitions between the two band insulators. Let us illustrate this point with a canonical model for Chern insulators: the Haldane model \cite{Haldane1988}. It consists of a honeycomb lattice model, with real nearest neighbours couplings (essentially the tight-binding model of graphene), plus complex hopping ${\rm e}^{{\rm i} \phi}$ terms to second nearest neighbours, thus breaking time-reversal symmetry, and a on-site staggered potential $\Delta$ that breaks inversion symmetry (see figure \ref{fig:Haldane} (a)).  This two-band model displays three topologically distinct phases with $C_\pm=0,\pm1,\mp1$, depending on the competition between the two gap opening mechanisms controlled by $\phi$ and $\Delta$, and are separated by  gap closing points (in orange and blue in figure \ref{fig:Haldane} (b)).

\begin{figure}[h!]
\centering
\includegraphics[width=14cm]{figure_Haldane}
\caption{\footnotesize{(a) Haldane model on the honeycomb lattice. (b) Phase diagram of the Haldane model. The topologically distinct insulating phases are characterized by the bulk Chern number  $C_\pm$ of each Bloch band over the Brillouin zone. They are separated by a gap closing point at either ${\mathbf{k}}=K$ or ${\mathbf{k}}=K'$. The gap amplitude is controlled by valley-dependent  mass terms $m_K$ and $m_{K'}$. Along the path $1$, the gap closes at a single valley}}.
\label{fig:Haldane}
\end{figure}

Those gap closing points occur around one of the two valleys $K$ and $K'$ of the graphene dispersion relation (figure \ref{fig:Haldane} (b)). At low energy (i.e. around $E=0$), this gap-opening management is controlled by two valley dependent mass terms  $m_K$ and $m_{K'}$. The topological transitions between two distinct Chern phases are thus  encoded into a low energy description around each valley $K$ or $K'$. In that limit, the Haldane model reduces to two 2D valley dependent Dirac Hamiltonians with a distinct masses  $m_K$ or $m_{K'}$, that depend on the two parameters of the lattice model, $\Delta$ and $\phi$. These Dirac Hamiltonians can be interpreted as symbol Hamiltonians parametrized respectively in $(\delta k_x, \delta k_y, m_{K^{(')}})$ parameter spaces with $\delta {\mathbf{k}}$ a small expansion around the valley $K^{(')}$.  Those symbol Hamiltonians have a Berry-Chern monopole that coincides with the transition between the two insulators, i.e. for $m_{K^{(')}}=0$. The interface problem of a Chern insulator with a trivial one is therefore captured by Hamiltonian operators ${\hat{\mathcal{H}}}_\text{op}$ when taking a mass varying term $m_K$ or $m_{K'}$  (path 1 in figure \ref{fig:Haldane}), say in the $x$ direction. The values $\mathcal{C}_\pm= \mp1$ of those Berry-Chern monopoles agree with the difference of bulk Chern numbers of the Bloch bands between the two insulators, and the associated spectral flow corresponds to the interface chiral mode. \\

Interestingly, the monopole-spectral flow correspondence can also be used to describe a transition between two trivial insulators (path 2 in figure \ref{fig:Haldane}). In that case, the bulk-edge correspondence predicts no interface state, because the band Chern numbers $C_\pm$ of the Bloch bands of both insulators are zero. However, the parameter $\Delta$ ,that breaks inversion symmetry, is varied when going from one insulator to the other one, and the gap closes when this parameter changes sign. At the gap closing point, the mass terms  associated to each valley $m_K$ and $m_{K'}$ vanish and one recovers  the usual graphene model with two Dirac points. The monopole-spectral flow correspondence is \textit{local} in parameter space, and thus applies for each band crossing point, and not to the bulk Bloch bands at fixed $\Delta$. Chern numbers $\mathcal{C}_\pm^{K}$ and $\mathcal{C}_\pm^{K'}$ of the Berry-Chern monopoles can thus be computed for each valley, like previously, and are actually opposite to each other. It follows that the interface problem, that is when $\Delta$ is varied in space and changes sign, manifests a double spectral flow which is opposite in each valley, as illustrated by the numerical spectrum of the Haldane model for a sharp interface in figure \ref{fig:QVHE}. This is the \textit{(quantum) valley Hall effect}, and the Berry-Chern monopole is reminiscent of the so-called \textit{valley} Chern number \cite{YaoNiu09, ZhangValley13}. This phenomenon triggered particular  interest in classical analogs of topological materials due to the experimental ease to break inversion symmetry \cite{Zhang2018, Ma2019_review, TopoPhot}. Note that because the monopole-spectral flow correspondence is local in parameter space, it actually does not require the system to be an insulator: gap inversion can be performed in a given valley, while the bands always touch somewhere else in the Brillouin zone, giving rise to topological spectral flows in (semi-)metallic like systems, where the band Chern numbers are ill defined \cite{Upreti20}.

\begin{figure}[h!]
\centering
\includegraphics[width=11cm]{figure_QVHE}
\caption{\footnotesize{Numerical spectrum of a sharp interface between Haldane models with opposite staggered potential $\Delta$, in a strip geometry, finite in the $x$ direction and infinite in the $y$ direction. the color represent the mean position of the eigenstate in the $x$ direction. The two valleys host opposite spectral flows, corresponding to valley polarized counter propagating modes localized in the middle of the strip where $\Delta$ changes sign. In contrast, the two flat bands edge states do not contribute to the spectral flow. }}
\label{fig:QVHE}
\end{figure}

Finally, we should also stress that the monopole-spectral flow correspondence is particularly suitable to investigate the existence of unidirectional  waves in continuous media \cite{Delplace2017, Perrot2019, Marciani20, Venaille20}. In the absence of a lattice, and thus without a Brillouin zone, there is  no natural close surface to integrate the Berry curvature over, and thus, in general, no well-defined band Chern number. Indeed, the Hamiltonian describing the homogeneous system depends then on the wavenumber ${\mathbf{k}} \in \mathbb{R}^2$, and the fiber bundle defined over this base space cannot always be compactified (into a torus or a sphere for instance). In other words, for a given waveband, the integration of the Berry curvature over $\mathbb{R}^2$ is in general not a Chern number, and when it is, the bulk-edge correspondence is modified in a subtle way \cite{Tauber_anomalous20}. However, in any case, the monopole Chern numbers $\mathcal{C}$ are still well-defined as soon as a band crossing point is found. Put differently, unidirectional interface states are still well defined as topological waves, even though the interface separates two bulk domains that  have ill-defined topology because of their  ill-defined band Chern numbers. In that case, the Berry-Chern monopoles are the natural and sufficient objects to consider.

So, to sum up, the same monopole-spectral flow correspondence naturally captures at the same time the interface states in  Chern insulators,  valley Hall insulators, semi-metallic systems, continuous media, as well as the bulk current of Weyl semi-metals in a magnetic field.

\subsubsection{A polariton-like toy model}

\label{sec:polariton}

The goal of this paragraph is two-fold: First, to illustrate how to construct a toy model of a wave coupled to a resonator that generates a topological spectral flow; and second, to  point out that a spectral flow does not necessarily imply  the existence of uni-directional modes.

Consider a non-dispersive wave of frequency $\omega=c k_x$ that travels through a resonator of frequency $\Omega$. If the wave travels freely in the media, (i.e. if it does not couple to the resonator), the resulting dispersion relation of the system is simply given by the superposition of the two branches $\omega_1=ck_x$ and $\omega_2=\Omega$. This obviously gives a linear band-crossing at $\Omega=ck_x$, but in a 1D parameter space, span by $k_x$ only. To get topologically guided modes out of this simple situation, one needs to embed this two-fold degeneracy point in a larger, i.e. three-dimensional parameter space $(k_x,k_y,m(y))$.  That is, One  needs to couple the wave with the resonator to make new modes (in the spirit of a polariton \cite{Pendry2004,Fano1956, Hopfield1956, Lemoult2016, Yves2017, Baboux2017, Kartashov2019}) but in such a way that it involves the wave number $k_y$ and an external or internal parameter $m$. Moreover, in order to embed the degeneracy point in the three-dimensional space, one sees that such a coupling $\gamma(k_y,m)$  must enter as a complex number, so that the Hamiltonian decomposes over the three Pauli matrices. One ends up with the following ''polariton-like'' toy model Hamiltonian
\begin{align}
\begin{pmatrix}
 c k_x & \gamma^*(k_y,m)  \\
 \gamma(k_y,m)   &  \Omega
\end{pmatrix}
\begin{pmatrix}
\varphi_1\\
\varphi_2
\end{pmatrix}
=\omega
\begin{pmatrix}
\varphi_1\\
\varphi_2
\end{pmatrix} \, .
\label{eq:polariton}
\end{align}
where $\varphi_1$ and $\varphi_2$ are complex numbers.

Its eigenfrequencies $\omega_\pm = \frac{ck_x+\Omega}{2}\pm  \sqrt{ \left( \frac{ck_x- \Omega}{2} \right)^2  +|2\gamma(k_y,m)|^2}$, are degenerated at the point $(k_x^{0} = \frac{\Omega}{c},k_y^{0},m^{0})$ with $\gamma(k_y^{0},m^0)=0$ by construction.
Because this is a two-band crossing point, it generates a non-zero Berry-Chern monopole in $(k_x,k_y,m)$ space. If, moreover, $\gamma(k_y,m)$ varies linearly with $k_y$ and $m$ around the degeneracy point, so that the relation dispersion is linear around the band crossing point, then its topological charge is $|\mathcal{C}_\pm| = 1$. For the sake of simplicity, we choose $\gamma=k_y+{\rm i} m$, and plot the eigenfrequencies in figure  \ref{fig:polariton}. With this choice, this classical toy model looks like the 2D spin $S=1/2$ massive Dirac Hamiltonian, except that the dispersion relation is tilted such that there is no direct gap.  This however does not prevent the existence of a spectral flow which is guarantied  in the dual inhomogeneous operator problem defined after the substitutions $m\rightarrow m(y)$ and $k_y\rightarrow -{\rm i} \partial_y$ provided $m(y)$ changes sign.\footnote{We have assumed for simplicity that the celerity $c$ is a constant. Otherwise, the correspondence with the operator is slightly more involved.} The value of the spectral flow $\Delta \mathcal{N}_\pm$ is given by $-\mathcal{C}_\pm$. The gain or the lost of a  mode (say)  in the band $+$,  is then fixed by the sign of $\mathcal{C}_+$. This sign of the Chern numbers $\mathcal{C}_\pm$ is reversed when reversing the mass term $m\rightarrow -m$ in the model, and accordingly the spectral flow is  reversed if $m(y)$ is a decreasing function of $y$ instead of an increasing one. As shown in figure \ref{fig:polariton} (b) and (c), this has a strong incidence on the physical behavior, since an unidirectional mode is obtained when $m(y)$ is an increasing function of $y$, while a non-propagating flat interface mode is obtained for a decreasing $m(y)$.  Note that in  both cases, there is not a direct gap. This  is the cause of the absence of two counter propagating interface modes when swapping the sign of the mass term.
\begin{figure}[h!]
\centering
\includegraphics[width=14cm]{figure_polariton}
\caption{\footnotesize{Frequency spectra of the polariton-like toy model for (a) the homogeneous (symbol Hamiltonian) problem, and (b) and (c)  the inhomogeneous (operator Hamiltonian) problem, for (b) $m(y)\propto y$ and (c) $m(y)\propto-y$. The dashed line indicates the absence of a direct gap.}}
\label{fig:polariton}
\end{figure}

The general model addressed in the next section also exhibits non-chiral, non uni-directional spectral flows of non-zero energy modes, and spectral flows without a gap in the spectrum of the operator,  in a  natural and systematic way.

\section{ Berry-Chern monopole - spectral flow correspondence beyond linear two-band crossings}
\label{sec:general}

\subsection{ $h.\hat{S}$ models}
\subsubsection{Motivations and generalities}
In the first part of this manuscript, the monopole-spectral flow correspondence was discussed from the example of a spin $1/2$ Hamiltonian, that is, for linear two-band crossings. Actually, the correspondence is neither restricted to two-fold band crossing points, nor to linear dispersion relations. In particular, a spectral flow of Landau levels was computed for a spin $S=1$ generalization of Weyl semi-metals \cite{Bradlyn2016} and even for arbitrarily higher values of $S$ \cite{Ezawa2017}. In the realm of classical physics, the spectral flow of oceanic and atmospheric equatorial waves obtained from the rotating shallow water model, that involves a three-band crossing point, was shown to have a topological interpretation from its symbol \cite{Delplace2017, Faure2019}. Remarkably, the spin $1$ semi-metals under a magnetic field on the one hand, and the rotating shallow water model  on the other hand, display both the exact same spectrum and carry the same Berry-Chern monopoles induce by three-fold degeneracy points. Although these two systems are physically  radically distinct, they are in fact both described with the same spin $S=1$ generalizations of the 3D Weyl semimetals and of the 2D massive Dirac particle examples, respectively. In that case, the Chern numbers take the values $\mathcal{C}=\{0,\pm2\}$, and the spectral flow increases accordingly, so that the correspondence is satisfied. This spectral flow corresponds to a transfer of two Landau levels in the first case, and  to the existence of two eastward propagating waves trapped along the equator, in the second one.

The goal of this section is to discuss the generalization of the monopole-spectral flow correspondence beyond both the spin $1/2$ case, and linear dispersion relations. For that purpose, we will focus on models, that we shall dub $h.\hat{S}$ \textit{models}, as their symbol Hamiltonians read
\begin{align}
{H_S} = \boldsymbol{h}({\boldsymbol{\lambda}})  . \, {\hat{\mathbf{S}}}
\label{eq:h.S}
\end{align}
with $\boldsymbol{h} \in \mathbb{R}^3$.  ${\boldsymbol{\lambda}} \in \mathbb{R}^3$ is a set of three \textit{parameters}, including classical conjugate variables that one needs to quantize to obtain the operator Hamiltonian. ${\hat{\mathbf{S}}}$ is abusively called the spin operator: $\hbar$ has been removed from the usual definition, but its three components form an irreducible representation of the su(2) algebra embedded in su(N)
$[\hat{S}_\alpha,\hat{S}_\beta]= {\rm i} \epsilon_{\alpha\beta\gamma} \hat{S}_\gamma$,
where $\epsilon_{\alpha\beta\gamma}$ is the Levi-Civita symbol and $\{\alpha, \beta, \gamma \} = \{ x,y,z\}$. It is worth stressing that the $h.\hat{S}$ models do not necessarily involve any actual spin. It turns out that such a structure also appears naturally in purely classical waves problems (an example is detailed in the appendix \ref{sec:appendix_SW}). What only matters  is this algebraic structure that is convenient to derive analytical results.

We choose $z$ as the quantization axis of ${\hat{\mathbf{S}}}$. Thus ${\hat{S}_z}$ is diagonal in the spin projection basis $\ket{m_S}$
\begin{align}
{\hat{S}_z} \ket{m_S} = m_S \ket{m_S}
\end{align}
with $m_S$ taking the $2S+1$ integer or half integer values
\begin{align}
m_S = -S, -S+1, \dots , S-1, S \ .
\end{align}
In the following we will also make use of the  spin ladder operators
\begin{align}
{\hat{S}_\pm} = {\hat{S}_x}  \pm {\rm i} {\hat{S}_y}
\label{eq:spin_ladder}
\end{align}
that satisfy
 \begin{align}
{\hat{S}_\pm} \ket{m_S} = \sqrt{S(S+1)-m_S(m_S\pm1)}\ket{m_S\pm1}  \, .
 \label{eq:rel_jpm}
 \end{align}

As we shall see, an elegant general expression of the Chern numbers of $h.\hat{S}$ models can be obtained analytically. After deriving this expression, we shall choose a specific $h.\hat{S}$ model, already introduced by Ezawa \cite{Ezawa2017}, for which we will not only directly obtain the Chern numbers but also investigate the operator Hamiltonian and its spectral flow. To do so, we will introduce a second index,  in the spirit of the index of Dirac operators, and compute it to verify its equality with the Chern number, in agreement by  the Atiyah-Singer theorem, which is the most advanced and abstract formulation of the monopole-spectral flow correspondence. Finally, we will construct the solutions of such a model and plot the different spectral flows.

\subsubsection{Chern numbers of the $h.\hat{S}$ models}
\label{sec:expression_chern}
It is instructive to decompose the computation of the Chern numbers $\mathcal{C}_m$ for the Hamiltonian $\eqref{eq:h.S}$ into two steps, by focusing first on the linear case $\mathbf{h}({\boldsymbol{\lambda}})={\boldsymbol{\lambda}}$ for an arbitrary spin $S$ before generalizing to  arbitrary functions $\mathbf{h}({\boldsymbol{\lambda}})$. This will allow us to distinguish two contributions to the Chern numbers.
\paragraph{Revisiting Berry's original problem.}

Let us start by setting $\mathbf{h}({\boldsymbol{\lambda}})={\boldsymbol{\lambda}}$ in \eqref{eq:h.S}. The Hamiltonian
\begin{align}
{H_S} = {\boldsymbol{\lambda}}  . \, {\hat{\mathbf{S}}}
\end{align}
is then formally equal to that of a quantum spin coupled to a magnetic field ${\boldsymbol{\lambda}}=-\mu_B \mathbf{B}$ with $\mu_B$ the Bohr magneton. This model was discussed by Michael Berry in its seminal paper  \cite{Berry1984} to introduce the geometrical phase accumulated by a quantum state $\ket{\psi_{m_S}}$ after varying   adiabatically the magnetic field along  a loop in space. The adiabaticity invoked here imposes the eigenstates $\ket{\psi_{m_S}}$ of ${H_S}$ to remain eigenstates during the evolution while  $\mathbf{B}$ is varied. In our formal model, ${\boldsymbol{\lambda}}$ is not physically varied in time; it is just a parameter of the Hamiltonian and of its eigenstates.

To compute the geometric phase of the eigenstates of ${H_S}$, Berry computed the 2-form Berry curvature. The Berry phase is then obtained integrating this curvature along a close path, while the Chern number is obtained after  integrating it over a close surface.
Berry's computation of the curvature is based on the expression \eqref{eq:practical_Berry} of the coefficients $F_{ij}^{({m_S})}({\boldsymbol{\lambda}})$, so that the curvature reads
 \begin{align}
 \label{eq:formule_Berry}
 \mathcal{F}_{m_S}({\boldsymbol{\lambda}}) = {\rm i}\sum_{p\neq {m_S}} \frac{\bra{\psi_{m_S}({\boldsymbol{\lambda}})}\partial_{\lambda_i} H_S \ket{\psi_p({\boldsymbol{\lambda}})}\bra{\psi_p({\boldsymbol{\lambda}})}\partial_{\lambda_j} H_S \ket{\psi_{m_S}({\boldsymbol{\lambda}})} }{(E_p({\boldsymbol{\lambda}})-E_{m_S}({\boldsymbol{\lambda}}))^2} {\rm d} \lambda_i \wedge {\rm d} \lambda_j
 \end{align}
This expression of the Berry curvature can be cumbersome to manipulate in cartesian coordinates, especially when one aims at integrating it over the sphere to get the Chern numbers. It is thus much more judicious to work in spherical coordinates, and the symbol Hamiltonian then becomes
 \begin{align}
{H_S} = |{\boldsymbol{\lambda}}|  \, \mathbf{n}(\theta,\phi) \cdot \, {\hat{\mathbf{S}}}
\end{align}
where the normalized vector $ {\boldsymbol{\lambda}} / |{\boldsymbol{\lambda}}| \equiv \mathbf{n}(\theta,\phi)=(\sin\theta \cos\phi,$ $\sin\theta \sin\phi, \cos\theta)$ gives the orientation of ${\boldsymbol{\lambda}}$ in $\mathbb{R}^3$, with $\theta$ and $\phi$ the usual  polar and azimuthal angles in the laboratory frame respectively. Our three abstract parameters $\{ \lambda_i\}$ in \eqref{eq:formule_Berry} must now be understood as the new coordinates $\{ |{\boldsymbol{\lambda}}|, \theta, \phi \}$. One now needs to express explicitly ${H_S}$ with those new variables, so that one can effectively differentiate it according to \eqref{eq:formule_Berry}.

Note that there is an ambiguity here about what $\ket{\psi_{m_S}}$ means. Indeed, $m_S$ refers to the projection of the spin onto a quantization axis. Usually, in problems involving a Zeeman-coupling, the  magnetic field's orientation is fixed, and canonically provides the axis of quantization. Here the situation is more subtle, since the orientation of the magnetic field, given by $\boldsymbol{n}(\theta,\phi)$, is not fixed, so one must precise which axis we refer to when writing $\ket{\psi_{m_S}}$. Let us thus fix the quantization axis, called $z$, in the laboratory frame. Then one writes $\ket{{m_S},\mathbf{n}}$ to refer to an eigenstate of ${H_S}$, and $\ket{m,\mathbf{e}_z}$ to refer to an eigenstate of ${\hat{S}_z}$ i.e.
\begin{align}
\label{eq:jz_Berry}
&{\hat{S}_z} \ket{m_S,\mathbf{e}_z} = m_S \ket{m_S,\mathbf{e}_z} \\
&{H_S} \ket{m_S,\mathbf{n}} = |{\boldsymbol{\lambda}}|\, m_S \ket{m_S,\mathbf{n}}
\end{align}
where ${\hat{S}_z} = \mathbf{e}_z \cdot {\hat{\mathbf{S}}}$. The eigenstates that intervene in the Berry curvature \eqref{eq:formule_Berry} must therefore be understood as $\ket{\psi_{m_S}}\equiv  \ket{m_S,\mathbf{n}}$.

The operator ${\hat{S}_z}$ is related to the Hamiltonian $ \mathbf{n}\cdot {\hat{\mathbf{S}}}$ by a rotation, represented by a unitary operator $U(\theta,\phi) $ that depends on the longitudinal and azimutal angles  as
\begin{align}
 \mathbf{n}(\theta,\phi)\cdot {\hat{\mathbf{S}}} \equiv U(\theta,\phi)\, \mathbf{e}_z \cdot {\hat{\mathbf{S}}} \, U^\dagger(\theta,\phi)
 = U(\theta,\phi)\,{\hat{S}_z} \, U^\dagger(\theta,\phi) \ .
\end{align}
This relates the eigenstates of the ''rotated'' Hamiltonian $H(\theta,\phi)$  to those of ${\hat{S}_z}$ as
\begin{equation}
\begin{aligned}
 \ket{m_S,\mathbf{n}} \equiv U(\theta,\phi)  \ket{m_S,\mathbf{e}_z} \ .
 \label{eq:rotation_ham}
\end{aligned}
\end{equation}
The rotation operator $U(\theta,\phi)$ is not uniquely defined. For instance, one can choose the representation
\begin{align}
U(\theta,\phi) = {\rm e}^{-{\rm i} \phi {\hat{S}_z}}{\rm e}^{-{\rm i} \theta {\hat{S}_y}} {\rm e}^{{\rm i} \phi {\hat{S}_z}} \ .
\label{eq:rotation_spin}
\end{align}
The Berry curvature \eqref{eq:formule_Berry} can now be computed. It implies terms of the form $\bra{m,\mathbf{n}}\partial_\theta {H_S} \ket{p,\mathbf{n}}$, $ \bra{m,\mathbf{n}}\partial_\phi {H_S} \ket{p,\mathbf{n}}$  and  $ \bra{m,\mathbf{n}}\partial_{|{\boldsymbol{\lambda}}|} {H_S} \ket{p,\mathbf{n}}$, but one can notice that all the terms involving a derivative with respect to the radius $|{\boldsymbol{\lambda}}|$ vanish, because they read $ \bra{m,\mathbf{e}_z} p \ket{p,\mathbf{e}_z}$ with $m\neq p$. We are thus left with a pure angle dependence only.
Recalling the expression of a rotation of ${\hat{S}_z}$ by an angle $\theta$ around the $y$ axis
\begin{align}
{\rm e}^{- {\rm i} \theta {\hat{S}_y}} {\hat{S}_z} {\rm e}^{{\rm i} \theta {\hat{S}_y}} = \cos \theta {\hat{S}_z} + \sin \theta {\hat{S}_x}
\end{align}
 one finds
\begin{align}
&\bra{m_S,\mathbf{n}}\partial_\theta {H_S} \ket{p,\mathbf{n}} =|{\boldsymbol{\lambda}}|\,  \bra{m_S,\mathbf{e}_z} {\hat{S}_x} \ket{p,\mathbf{e}_z} \\
&\bra{p,\mathbf{n}}\partial_\phi  {H_S} \ket{m_S,\mathbf{n}} =|{\boldsymbol{\lambda}}|\,  \sin\theta\bra{p,\mathbf{e}_z} {\hat{S}_y} \ket{m_S,\mathbf{e}_z} \ .
\end{align}
 These quantities can be computed easily by using the spin ladder operators
through the substitution
 ${\hat{S}_x}=\frac{1}{2}({\hat{S}_+}+{\hat{S}_-})$ and ${\hat{S}_y} = \frac{1}{2{\rm i}}({\hat{S}_+}-{\hat{S}_-})$, and one ends up with the important result
 \begin{align}
\mathcal{F}_{m_S} =  - m_S\, \sin \theta \, {\rm d} \theta \wedge {\rm d} \phi
\label{eq:curvature_spin}
 \end{align}
 that was first found by Berry \cite{Berry1984}.
In this expression, one recognizes (see appendix \ref{sec:diff}) the surface form $4\pi \Omega_{S^2} = \sin \theta {\rm d} \theta \wedge {\rm d} \phi$ of the sphere $S^2$ in spherical coordinates  that satisfies $\int_{S^2}\Omega_{S^2}=1$.
The Chern numbers are then readily obtained  \eqref{eq:curvature_spin} over $S^2$ so that
 \begin{align}
 \mathcal{C}_{m_S} = -2m_S \ .
 \label{eq:chern_berry_mon}
 \end{align}
 This result was originally derived in a different way by Avron, Sadun, Segert and Simon \cite{Avron1989}.
 In the case of a spin $1/2$, one recovers the result $\mathcal{C}_\pm = \mp 1$   for the two eigenstates $m_S=\pm 1/2$  as derived previously. Higher values of the Chern number can then be reached for higher spins.

\paragraph{The \textit{Hamiltonian map} and its degree.}

Let us now reintroduce the full  $\mathbf{h}({\boldsymbol{\lambda}})$ dependence in the problem.
We shall refer to ${\boldsymbol{\lambda}} \rightarrow \mathbf{h}({\boldsymbol{\lambda}})$ as the  \textit{Hamiltonian map}  $h$ between the parameter space  of the ${\boldsymbol{\lambda}}$'s and a target space $\mathbb{R}^3$ span by $\mathbf{h}$. In others words, the image of ${\boldsymbol{\lambda}}$ by the Hamiltonian map $h$ is represented by a vector $\mathbf{h}\in\mathbb{R}^3$ (see figure \ref{fig:surface}).
As motivated in section \ref{sec:topo_defect}, one can view the degeneracy points ${\boldsymbol{\lambda}}_0$ as defects in parameter space. Those defects can be characterized  by how the Hamiltonian map transforms a surface $\mathcal{S}_{{\boldsymbol{\lambda}}_0}$ enclosing  ${\boldsymbol{\lambda}}_0$, that is by considering
\begin{align}
{\boldsymbol{\lambda}} \in  \textcolor{blue}{\mathcal{S}_{{\boldsymbol{\lambda}}_0}  \in \mathbb{R}^{3}} \rightarrow \mathbf{h}({\boldsymbol{\lambda}})  \in \textcolor{red}{\mathcal{S}_h \in \mathbb{R}^{3}}
\label{eq:S-S}
\end{align}
where $\mathcal{S}_{h}$ is a surface spanned by $\boldsymbol{h}$ in target space \cite{}.
Since those two surfaces are oriented manifold and have the same dimension, the maps $h$ between those spaces are  classified according to their \textit{degree} that tells whether $\mathcal{S}_{h}$ is a closed surface surrounding the origin $\mathbf{h}=0$ and how many times it encloses it \cite{DubrovinBook}. It reads
\begin{align}
\deg h \equiv
 \sum_{{\boldsymbol{\lambda}}^{(i)}} \text{sgn} \ensuremath{\operatorname{det}}
\begin{pmatrix}
\partial_{\lambda_1}h_x & \partial_{\lambda_1}h_y & \partial_{\lambda_1}h_z \\
\partial_{\lambda_2}h_x & \partial_{\lambda_2}h_y & \partial_{\lambda_2}h_z \\
\partial_{\lambda_3}h_x & \partial_{\lambda_3}h_y & \partial_{\lambda_3}h_z \\
\end{pmatrix}\Big|_{\mathbf{h}_0}
\label{eq:degree}
\end{align}
where the matrix whose determinant is computed is the Jacobian matrix between the parameter space and the target space.
The points ${\boldsymbol{\lambda}}^{(i)}$,  which the sum runs over, are called the pre-images of $h$, and are associated to a given arbitrary direction $\mathbf{n}_0 \equiv  \mathbf{h}_0/|\mathbf{h}_0| $, as they satisfy $\mathbf{h}({\boldsymbol{\lambda}}^{(i)})/|\mathbf{h}({\boldsymbol{\lambda}}^{(i)})| =\mathbf{n}_0$. So in other words, the pre-images ${\boldsymbol{\lambda}}^{(i)}$ are all the points of parameter space whose image $\mathbf{h}({\boldsymbol{\lambda}}^{(i)})$ points in a fixed direction $\mathbf{n}_0$ (see figure \ref{fig:surface}).
 The value of the degree does not depend on the choice of that direction, provided the determinant is not singular (i.e. nonzero) for that choice.

\begin{figure}[htb]
 \centering
\includegraphics[width=15cm]{figure_surface_degree}
\caption{\footnotesize{Sketch of the Hamiltonian map $h$ between the parameter space and the target space. The degree of $h$ is defined when considering the image $ \textcolor{red}{\mathcal{S}_h}$ of a surface $ \textcolor{blue}{\mathcal{S}_{{\boldsymbol{\lambda}}_0}}$ that encloses the degeneracy point ${\boldsymbol{\lambda}}_0$. The points ${\boldsymbol{\lambda}}^{(1)}$, ${\boldsymbol{\lambda}}^{(2)}$ and ${\boldsymbol{\lambda}}^{(3)}$ are examples of pre-imagines that enter the formula \eqref{eq:degree} of the degree. }}
\label{fig:surface}
\end{figure}

The degree is obviously an integer number. This is a topological index in that it is invariant under continuous changes of $h$. Two maps $h_1$ and $h_2$ having the same degree can thus be continuously deformed one into each others; they are said to be homotopic. This homotopy invariance makes free the choice of surfaces $\mathcal{S}_{{\boldsymbol{\lambda}}_0}$ that surrounds the degeneracy point ${\boldsymbol{\lambda}}_0$ so that one can choose it as being the unit sphere (see figure \ref{fig:surface}). The degree actually classifies the maps from spheres to spheres
\begin{align}
{\boldsymbol{\lambda}} \in  \textcolor{blue}{S^2} \rightarrow \mathbf{h}/ |\mathbf{h}|  \in \textcolor{red}{S^2}
\label{eq:S2-S2}
\end{align}
according to their homotopy properties. They are the elements of the homotopy group $\pi_2(S^2)=\mathbb{Z}$ (and more generally that of $\pi_n(S^n)=\mathbb{Z}$).
As for the maps $h$ we deal with, the degree is known to have an elegant geometrical interpretation: it counts how many times the unit vector $\mathbf{h}/ |\mathbf{h}|$ in \eqref{eq:S2-S2} wraps the unit sphere (see appendix \ref{sec:diff} for more details).

To see explicitly how the degree alters the value of the Chern number, one must interpret the result of the previous section as if the reasoning was performed in the target space. Namely, the expression of the Berry curvature \eqref{eq:curvature_spin} must  be understood in the space of $\mathbf{h}$ as if the underlying ${\boldsymbol{\lambda}}$ dependence was ignored. However, we want to express the Chern number as an integral of the Berry curvature over a surface (a sphere) that surrounds the degeneracy point in parameter space $\Lambda$ as
\begin{align}
\mathcal{C}_{m_S} =& \frac{1}{2\pi} \int_{\textcolor{blue}{S^2}} \mathcal{F}_{m_S} \\
=& -2m_S \int_{\textcolor{blue}{S^2}} h^{\star}\Omega_{S^2}
\end{align}
where $h^{\star}\Omega_{S^2}$ denotes the \textit{pull-back} of $\Omega_{S^2}$ by $h$; this is a procedure to formally define a differential form in parameter space when it is already defined in target space (see appendix \ref{sec:diff}).
Actually this formal manipulation simplifies easily thanks to the Brouwer theorem, that relates the integrals of a differential form to that of its ''pulled-back'' as
\begin{align}
\mathcal{C}_{m_S}= -2{m_S}  \deg{h}  \int_{\textcolor{red}{S^2}} \Omega_{S^2}
\end{align}
where the integral acts now in target space and where $ \deg{h} $ is precisely the degree of $h$. Again, the integral that remains is $1$ by definition, which leads to the final expression of the topological charges of the Berry-Chern monopoles
 \begin{align}
\mathcal{C}_{m_S} = -2{m_S}\, \deg h \ .
\label{eq:chern_berry_mon_gen}
\end{align}
This is the general and practical expression of the Chern-Berry monopole for $h.\hat{S}$ models \cite{DelplaceHDR, TauberGraf21}.

One can easily check that one recovers $\mathcal{C}_{\pm}=\mp1$ for the $2\times 2$ Hamiltonians (${m_S}=\pm1/2$), since the degree \eqref{eq:degree} in that case is $1$.
It is worth noticing that the formula \eqref{eq:chern_berry_mon_gen} clearly reveals the two concurrent mechanisms that give the Chern numbers their values~: the order of the degeneracy (property of ${\hat{\mathbf{S}}}$) and the wrapping of the Hamiltonian map (property of $\mathbf{h}({\boldsymbol{\lambda}})$).

\subsection{A specific  $h.\hat{S}$ model}

In this section, we propose a  specific $h.\hat{S}$ model to decorticate the correspondence between the spectral flow of an operator Hamiltonian and the Chern numbers of the fiber bundles associated to the eigenstates of the symbol Hamiltonian.  This model generalizes what was discussed in the first part of the manuscript  to (1) arbitrary integer and half-integer spins  $S\geqslant1/2 $ and (2) higher order ladder operators ${\hat{a}} \rightarrow {\hat{a}}^d$ with $d$ a positive integer. The operator Hamiltonian we consider is a composition of spin ladder operators ${\hat{S}_\pm}$ and bosonic ladder operators ${\hat{a}}$ and ${\hat{a}^{\dagger}}$ that reads%
\begin{align}
\label{eq:gen_hop}
{\hat{\mathcal{H}}}_\text{op} =\left(({\hat{a}^{\dagger}})^d \, {\hat{S}_+} +  {\hat{a}}^d  \, {\hat{S}_-} \right)  + \lambda {\hat{S}_z} \ .
\end{align}

Up to a factor $2$, this new ${\hat{\mathcal{H}}}_\text{op}$ coincides with the previously studied operator Hamiltonian (Eqs \eqref{eq:H_Dop} and \eqref{eq:H_Wop}) when $S=1/2$ and $d=1$. It acts on a Hilbert space that decomposes as the product of a Fock space $\mathscr{F}$, where the number states $\ket{n}$ live, and a spin space $\mathscr{S}$ spanned by the spin projections  $\ket{m_S}$. A vector of this space $\mathscr{F} \otimes \mathscr{S}$ thus reads $ \ket{n} \otimes \ket{m_S} \equiv \ket{n,m_S}$. The action of ${\hat{\mathcal{H}}}_\text{op}$ on those states is then to shift the different components $n$ and $m_S$; it thus  somehow describes a dynamics on a lattice of coordinates $(n,m_S)$.

As we shall see, this model exhibits very rich spectral flows, depending on $S$ and $d$. In the limit $|\lambda| \gg1$, which is precisely when the spectral flow modes join a branch, the term $\lambda {\hat{S}_z}$ is dominant, and  one expects therefore the eigenstates $\ket{\Psi}$ of ${\hat{\mathcal{H}}}_\text{op}$ to resemble the spin states $ \ket{m_s}$, because  the typical energy between two consecutive number states $\ket{n}$ and $\ket{n+1}$  in the same branch must tend to zero. There is then a clear separation of energy (or time) scales, such that the spin contribution $\lambda {\hat{S}_z}$ (that one can associate to a fast motion), is dominant over the ${\hat{a}}^d {\hat{S}_-}$ and ${\hat{a}^{\dagger}} {\hat{S}_+}$ terms that imply annihilation and creation operators (and that one can thus associate to a slow motion).
With this in mind, it is tempting, in order to capture the spectral flow occurring  for $|\lambda| \gg1$, to approximate the system by an ''adiabatic'' model, by substituting the operator ${\hat{a}}^d$ by its symbol.
Such a transformation is \textit{a priori} not straightforward, since the symbol of a product of operators is in general not equal to the product of the operators of each symbols  \cite{FollandBook, ZworskiBook}. It turns out that this is however the case in our model, and we have the simple correspondences
\begin{align}
\label{eq:symbol_ad1}
&(\lambda_1+{\rm i} \lambda_2)^d \quad \leftrightarrow \quad \ {\hat{a}}^d\\
&(\lambda_1- {\rm i} \lambda_2)^d \quad \leftrightarrow \quad ({\hat{a}^{\dagger}})^d
\label{eq:symbol_ad2}
\end{align}
where $\lambda_1$ and $\lambda_2$ are canonical conjugate variables in phase space. This mapping is a consequence of a result derived long ago by Mac Coy \cite{McCoy32}, and that is re-derived by induction in the appendix \ref{sec:operator-symbol}.

One thus obtains the symbol Hamiltonian of the expected form ${H_S}=\boldsymbol{h}({\boldsymbol{\lambda}})  . \, {\hat{\mathbf{S}}} $ as announced in \eqref{eq:h.S}, with
\begin{subequations}
\label{eq:gen_symbol2}
\begin{align}
&h_x({\boldsymbol{\lambda}}) = 2\,\text{Re} (\lambda_1-{\rm i} \lambda_2)^d \\
&h_y({\boldsymbol{\lambda}}) = - 2\,\text{Im} (\lambda_1-{\rm i} \lambda_2)^d \\
&h_z({\boldsymbol{\lambda}}) = \lambda \,.
\end{align}
\label{eq:gen_symbol_model}
\end{subequations}

The Hilbert space of the symbol Hamiltonian is reduced compare to that of ${\hat{\mathcal{H}}}_\text{op}$, since only the $2S+1$ spin (fast) degrees of freedom are kept, but is parametrized by the slow degrees of freedom that are accounted by the commuting (classical) variables. In that picture, the space of fast motion is a fiber bundle over the phase space of slow variables \cite{Faure2001, Faure2019}.

In the next section, we introduce the analytical index of the operator Hamiltonian \eqref{eq:gen_hop}, relate it to the spectral flow, compute it and compare its values to that of the Chern numbers  of ${H_S}$ with \eqref{eq:gen_symbol_model}.

\subsection{Analytical and Topological indices, a verification of the Atiyah-Singer theorem}
\subsubsection{Chern numbers of the specific $h.\hat{S}$ model}

The expression \eqref{eq:chern_berry_mon_gen} of the Chern numbers is quite general. In this section we compute it  for the specific symbol Hamiltonian given by Eqs \eqref{eq:gen_symbol2}.
One thus needs to compute the degree of such a Hamiltonian map $h$, which requires the determinant of the Jacobian matrix \eqref{eq:degree} :
\begin{align}
\ensuremath{\operatorname{det}}  \begin{pmatrix}
\partial_{\lambda_1} h_x  & \partial_{\lambda_1} h_y & 0 \\
\partial_{\lambda_2} h_x  & \partial_{\lambda_2} h_y & 0 \\
0  & 0 & 1
\end{pmatrix}
= \frac{\partial h_x}{\partial \lambda_1}  \frac{\partial h_y}{\partial \lambda_2} - \frac{\partial h_x}{\partial \lambda_2}  \frac{\partial h_y}{\partial \lambda_1} \ .
\label{eq:det}
\end{align}
We have
\begin{align}
(\lambda_1-{\rm i} \lambda_2)^d = \sum_{\alpha=0}^d \begin{pmatrix}  d \\ \alpha \end{pmatrix} \lambda_1^{d-\alpha} (-{\rm i} \lambda_2)^\alpha
\end{align}
that we split into even $\alpha=2k$ and odd $\alpha=2k+1$ contributions to get the real and imaginary parts, so that
\begin{align}
h_x = & \, 2\,  \sum_{\substack{k\geqslant 0 \\ 2k  \leqslant  d}} \begin{pmatrix}  d \\ 2k \end{pmatrix} (-1)^k \lambda_1^{d-2k}\lambda_2^{2k} \\
h_y = & \,  2 \sum_{\substack{k\geqslant 0 \\ 2k+1  \leqslant  d}} \begin{pmatrix}  d \\ 2k+1 \end{pmatrix} (-1)^k \lambda_1^{d-2k-1}\lambda_2^{2k+1}  \ .
\end{align}
To compute the determinant \eqref{eq:det}, one then needs to take the derivatives of these terms. One has
\begin{align}
\frac{\partial h_x}{\partial \lambda_2} &= 2 \sum_{\substack{k\geqslant 1 \\ 2k  \leqslant  d}} (-1)^k\begin{pmatrix}  d \\ 2k \end{pmatrix} (2k) \lambda_1^{d-2k}\lambda_2^{2k-1}\\
\frac{\partial h_y}{\partial \lambda_1} &= 2 \sum_{\substack{k\geqslant 0 \\ 2k+1  \leqslant  d-1}} (-1)^k\begin{pmatrix}  d \\ 2k+1 \end{pmatrix} (d-2k-1) \lambda_1^{d-2k-2}\lambda_2^{2k+1}
\end{align}
By setting $k'=k+1$, one can rewrite the second equation as
\begin{align}
\frac{\partial h_y}{\partial \lambda_1} = -2 \sum_{\substack{k'\geqslant 1 \\ 2k'  \leqslant  d}} (-1)^{k'}\begin{pmatrix}  d \\ 2k'-1 \end{pmatrix} (d-(2k'-1))  \lambda_1^{d-2k'}\lambda_2^{2k'-1}
\end{align}
and notice that
\begin{align}
\begin{pmatrix}  d \\ 2k'-1 \end{pmatrix} (d-(2k'-1))
&=\begin{pmatrix}  d \\ 2k \end{pmatrix} 2k
\end{align}
so that
\begin{align}
\frac{\partial h_y}{\partial \lambda_1} = - \frac{\partial h_x}{\partial \lambda_2} \, .
\end{align}
Repeating this algebraic gymnastic for the cross derivatives terms yields  now
\begin{align}
\frac{\partial h_x}{\partial \lambda_1} = \frac{\partial h_y}{\partial \lambda_2}
\end{align}
and thus the determinant \eqref{eq:det} satisfies
\begin{align}
\ensuremath{\operatorname{det}}  \begin{pmatrix}
\partial_{\lambda_1} h_x  & \partial_{\lambda_1} h_y & 0 \\
\partial_{\lambda_2} h_x  & \partial_{\lambda_2} h_y & 0 \\
0  & 0 & 1
\end{pmatrix}
=  \left( \frac{\partial h_x}{\partial \lambda_1} \right)^2 + \left(\frac{\partial h_y}{\partial \lambda_1} \right)^2 >0
\label{eq:det_pos}
\end{align}
Therefore, each pre-image ${\boldsymbol{\lambda}}^{(i)}$ contributes as $+1$ to to the degree \eqref{eq:degree}, that is
\begin{align}
\deg h = \sum_{\lambda^{(i)}} 1 \ .
\end{align}
The computation of the degree now simply consists in counting the number of pre-images.
To complete this calculation, one must choose a direction $\mathbf{n}_0$ in the target space.
The choice $\mathbf{n}_0=\mathbf{e}_z$ is not an option, since it would imply that the expression \eqref{eq:det_pos} of the determinant must be evaluated at $\lambda_1=\lambda_2=0$ (as soon as $d>1$), which yields a vanishing determinant.
Instead, one can choose any direction in the plane $(\mathbf{e}_x,\mathbf{e}_y)$, such that, e.g. $\mathbf{n}_0=\mathbf{e}_y$. Using the polar coordinates $ (\lambda_1-{\rm i} \lambda_2) \equiv z =|z|{\rm e}^{{\rm i} \theta} $, the number of pre-imagines to be determined is the number of values $\theta_p$ that satisfy
\begin{align}
h_x &= 2 |z|^d \cos(\theta_p\, d) = 0\\
h_y &= 2 |z|^d \sin(\theta_p\, d) >0 \ .
\end{align}
This involves
\begin{align}
\theta_p =\frac{\pi}{2d} + 2\pi \frac{p}{d}
\end{align}
with $p$ an integer. It follows that $\theta$ can take $d$ different values on the circle to satisfy those conditions, this means that there are $d$ pre-images.

Finally, the degree of the hamiltonian map is therefore $\deg h = d$ and the values of the Chern numbers of the symbol Hamiltonian for the model \eqref{eq:gen_symbol2} are
\begin{align}
\mathcal{C}_ {m_S} = -2{m_S}\, d \ .
\label{eq:chern_model}
\end{align}

\subsubsection{Spectral flow as an analytical index of ${\hat{\mathcal{H}}}_\text{op}$}

We would like to introduce an index \cite{AS1963, Callias1978, Weinberg1981, Adams2010} that accounts for the spectral flow of the operator Hamiltonian ${\hat{\mathcal{H}}}_\text{op}$ that we have introduced in equation \eqref{eq:gen_hop}. We expect then that this index to be related to the the first Chern numbers \eqref{eq:chern_model}, through the celebrated  Atiyah-Singer index theorem \cite{AS1963, APS1975, GiampieroBook, NakaharaBook}.

In the first part, we introduced the notion of spectral flow from the canonical example of a two-band model (${\hat{\mathcal{H}}}_\text{op}$ with $S=1/2$ and $d=1$) in a handwavy way,  as the number of energy levels that transit from one branch to another when a parameter $\lambda$ is continuously varied.\footnote{See for instance page 99 of \cite{APS1976} for a much more rigorous and abstract definition.}  For this same model, the spectral flow can be defined more formally by a \textit{spectral index} that counts the net number of energy levels that move upward in the vicinity of $\tilde{E}=0$ and $\lambda=0$, when sweeping $\lambda$. The existence of such an index is meaningful provided the spectrum $\tilde{E}$ of ${\hat{\mathcal{H}}}_\text{op}$ is discrete within that window. This discreteness is guaranteed for this model since the symbol Hamiltonian has the particular property of being an \textit{elliptic} operator; its determinant is non-zero except at the degeneracy point \cite{Faure2019, NakaharaBook}.
This spectral index is however not obviously suitable beyond $S>1/2$, since more than two bands are involved. Besides, as we shall see with the specific  ${\hat{\mathcal{H}}}_\text{op}$ we introduced in \eqref{eq:gen_hop}, other spectral flows can exist between the different branches, without crossing the spectral gap around $\tilde{E}=0$ and  $\lambda=0$. We thus need to follow a different strategy.

Alternatively, we could also define an index that counts the number of energy modes that join or leave each branch when $|\lambda| \rightarrow \infty$, without specifying an energy window around the degeneracy point of ${H_S}$.  For that purpose, one can get some insights from another quite standard procedure to define a spectral flow, although more abstract, that was
developed for Dirac Operators \cite{Weinberg1981, BerlineBook, NakaharaBook}.  Dirac operators are first order differential operators whose square is a Laplacian. They decompose over $2N+1$ $\gamma$ matrices which satisfy  the Clifford algebra $\gamma^j  \gamma^k + \gamma^k  \gamma^j=2\delta^{jk} \mathbb{I}_{2N}$. This algebraic structure guarantees a chiral symmetry; namely,  there necessarily exists a unitary operator $\Gamma$ that can be written as $\Gamma=\text{diag}(\mathbb{I}_N,-\mathbb{I}_N)$ and that anticommutes with the Dirac Hamiltonian. It follows that the Dirac Hamiltonian can always take an off-diagonal form as
\begin{align}
\label{eq:dirac}
\mathcal{H}_{Dirac}=
\begin{pmatrix}
0  & D^\dagger \\
D & 0
\end{pmatrix} \ .
\end{align}
For instance, in the particular case $N=1$, the $\gamma$ matrices reduces to the Pauli matrices and $D=\partial_x + A(x)$ with $A(x)$ a real function  \cite{NakaharaBook}. The chiral structure of the Dirac Hamiltonian splits the Hilbert space into two chiral subspaces $E_-\oplus E_+$ defined from the spectral projectors $ \mathcal{P}_+\equiv \frac{1}{2}(\mathbb{I}_{2N}+\Gamma) = \text{diag}(\mathbb{I}_N,0_N)$ and $ \mathcal{P}_-\equiv \frac{1}{2}(\mathbb{I}_{2N}-\Gamma)= \text{diag}(0_N,\mathbb{I}_N)$. A state living in $E_\pm$ is thus an eigenstate of $\Gamma$ with eigenvalue $\pm1$, called the chirality.  The two chiral subspaces are related to each other as
\begin{align}
E_- \overset{ \hat{D}^{\dagger} }{\underset{ \hat{D}}\rightleftarrows} E_+
\label{eq:complexe}
\end{align}
where
\begin{align}
 \hat{D}  =
  \begin{pmatrix}
0  & 0 \\
D & 0
\end{pmatrix}
 , \qquad
 \hat{D}^\dagger=
 \begin{pmatrix}
0  & D^\dagger \\
0 & 0
\end{pmatrix} \ .
\end{align}
Structures such as \eqref{eq:complexe}, with nilpotent operators $\hat{D}$, are referred to as \textit{complexes}  \cite{NakaharaBook}. They are associated to an \textit{index} defined as \cite{Weinberg1981, BerlineBook, NakaharaBook}
\begin{align}
\text{ind}\,  \hat{D} \equiv \text{dim} \, \text{Ker} \,  \hat{D}   -   \text{dim} \, \text{Ker} \,  \hat{D}^{\dagger}
\label{eq:anal_index}
\end{align}
which is well-defined when the right-hand side member is finite. In that case, the index is obviously an integer number, and is known to be robust under homotopic deformations. In other words, multiplying $hat{D}$ by a nonzero number or by an operator which is connected to Identity (such as a unitary) does not change the value of the index.

As we shall see in a moment, $\text{ind}\,  \hat{D}$ yields an important information about the spectrum of the $\mathcal{H}$.
Indeed, this index counts the unbalance of zero energy modes of opposite chiralities.   The occurrence of such so-called zero modes is a fundamental property of fermions coupled to gauge fields.
Their investigation for Dirac operators has stimulated tremendous efforts in high and low energy physics and in mathematics. Importantly for our purpose, the index of Dirac operators is also related to the spectral flow of zero-modes.
It is then tempting to apply this machinery to ${\hat{\mathcal{H}}}_\text{op}$ for $S=1/2$ and $d=1$, that is when the Hamiltonian operator ${\hat{\mathcal{H}}}_\text{op}$ reduces to that of 3D Weyl semimetals  in a magnetic field \eqref{eq:H_Wop} as well as that of  2D Dirac fermions with an anisotropic mass \eqref{eq:H_Dop}. Indeed, in that case, and when we take $\lambda=0$, ${\hat{\mathcal{H}}}_\text{op}$ has a similar form as \eqref{eq:dirac} with $D={\hat{a}}$, $D^\dagger={\hat{a}^{\dagger}}$, and the index $\text{ind}\,  \hat{D}$ can easily be found to be $1$ by a straightforward calculation, in agreement with the spectral flow previously computed.
However, this procedure seems to work only when $\lambda=0$, since the chiral structure is lost in the presence of $\lambda$.
Moreover, and more importantly, ${\hat{\mathcal{H}}}_\text{op}$ does not satisfy anyway such a chiral symmetry beyond $S=1/2$, meaning that ${\hat{\mathcal{H}}}_\text{op}^2$ is not a Laplacian, even when $\lambda=0$. In other words, ${\hat{\mathcal{H}}}_\text{op}$ is not a Dirac operator.

 One can nevertheless circumvent this difficulty and generalize the previous approach by writing the operator Hamiltonian as
\begin{align}
{\hat{\mathcal{H}}}_\text{op} = \mathcal{D} + \mathcal{D}^\dagger + \lambda {\hat{S}_z}
\end{align}
with
\begin{align}
\mathcal{D} \equiv {{\hat{a}}}^d {\hat{S}_-}. \qquad \mathcal{D}^\dagger = ({\hat{a}^{\dagger}})^d {\hat{S}_+}
\end{align}
that allows us to decompose the Hilbert space as
\begin{align}
E_{-S}\overset{ \mathcal{D}^{\dagger} }{\underset{ \mathcal{D}}\rightleftarrows}
\dots
\overset{ \mathcal{D}^{\dagger} }{\underset{ \mathcal{D}}\rightleftarrows}
 E_{m_S-1}\overset{ \mathcal{D}^{\dagger} }{\underset{ \mathcal{D}}\rightleftarrows}  E_{m_S}  \overset{ \mathcal{D}^{\dagger} }{\underset{ \mathcal{D}}\rightleftarrows} E_{m_S+1}.  \overset{ \mathcal{D}^{\dagger} }{\underset{ \mathcal{D}}\rightleftarrows} \dots.  \overset{ \mathcal{D}^{\dagger} }{\underset{ \mathcal{D}}\rightleftarrows}\, E_S
\label{eq:complexe_gen}
\end{align}
where the  subspaces $E_{m_S}$ are span by states $\ket{\Psi}=\ket{p,{m_S}}$ with arbitrary $p\in \mathbb{N}$ and fixed ${m_S}$. Actually, one can also associate an index to such a structure that counts the spectral flow of ${\hat{\mathcal{H}}}_\text{op}(\lambda)$.

Indeed, the spectral flow $\Delta \mathcal{N}_{m_S}$ measures the unbalance of number states $p$ inside a branch ${m_S}$, when $\lambda$ is swept from $-\infty$ to $+\infty$. One can use $\mathcal{D}$ to introduce an index that measures this unbalance, by noticing that ${\hat{\mathcal{H}}}_\text{op}$ is invariant under the joint transformations ${\boldsymbol{\lambda}} \rightarrow -{\boldsymbol{\lambda}}$ and ${\hat{S}_z} \rightarrow -{\hat{S}_z}$. In other words, flipping $\lambda$ has the same effect as flipping ${\hat{S}_z}$. Since this  transformation amounts to perform $m_S\rightarrow -m_S$, then evaluating the spectral flow in a branch $m_S$ thus amounts to counting the pairing of states  $\ket{p,m_S}$ with $\ket{p',-m_S}$. Those states respectively live in $E_{m_S}$ and $E_{-m_S}$ spaces, that are related to each other by $\mathcal{D}^{2m_S}$, as sketched in Fig. \ref{fig:index}, and we have
 \begin{align}
 \label{eq:d2m}
   \mathcal{D}^{2m_S}&\ket{p,m_S} \ \ \ = \prod_{j=1}^d {\sqrt{p+2m_S(j-1)}} \ket{p-2m_Sd,-m_S} \\
   \label{eq:ddag2m}
 (\mathcal{D}^\dagger)^{2m_S}&\ket{p',-m_S}=  \prod_{j=1}^d {\sqrt{p'+2m_Sj}}\ket{p'+2m_S d,m_S} \ .
 \end{align}
 %
Actually, the prefactors in Eqs. \eqref{eq:d2m} and \eqref{eq:ddag2m} are unimportant. What is important is the asymmetry of the action of  $  \mathcal{D}^{2m_S}$ and $(\mathcal{D}^\dagger)^{2m_S}$.
Indeed, Eq. \eqref{eq:ddag2m} indicates that any mode $p'$ in the branch $m_S$ for $\lambda<0$  can be paired with a mode $p=p'+2m_S d$ in the same branch for $\lambda>0$, because $p$ cannot be negative. In other words, $(\mathcal{D}^\dagger)^{2m_S}$ is injective, and we have
 \begin{align}
  \text{dim} \, \text{Ker} (\mathcal{D}^\dagger)^{2m_S} = 0 \ .
 \end{align}
In contrast,
$  \mathcal{D}^{2m_S}$ associates a state $\ket{p',-m_S} \in E_{-m_S}$ for each state $\ket{p,m_S} \in E_{m_S}$ with $p'=p-2m_Sd$, but in a non-injective way, since
 \begin{align}
 \ket{p-2m_S d,-m_S}=0\qquad  \text{for} \quad p=0,\dots,2m_S d-1.
 \end{align}
 It follows that the number of modes at $\lambda>0$ in the branch $m_S$ that do not have a pairing partner in the same branch at $\lambda<0$ through the action of $\mathcal{D}^{2m_S}$, is given by
 \begin{align}
 \text{dim} \, \text{Ker} \mathcal{D}^{2m_S} = 2m_Sd \ .
 \end{align}
The figure \ref{eq:unbalance} summarizes this peculiar structure induced by $\mathcal{D}$ being \textit{almost}  invertible.
The unbalance $\Delta \mathcal{N}_{m_S}$ of modes in a branch $m_S$ between $\lambda>0$ and $\lambda<0$  is finally given by the net number of un-paired modes between $E_{m_S}$ and $E_{-m_S}$ as
\begin{align}
\label{eq:flow_index}
\Delta \mathcal{N}_{m_S} =   \text{dim} \, \text{Ker} \mathcal{D}^{2m_S} -  \text{dim} \, \text{Ker} (\mathcal{D}^\dagger)^{2m_S} \ .
\end{align}

 \begin{figure}[h!]
 \centering
\includegraphics[width=13cm]{figure_Index_spin_complex}
\caption{\footnotesize{Sketch of the mappings $(\mathcal{D}^\dagger)^{2m_S}$ and $\mathcal{D}^{2m_S}$ between the two subspaces $E_{-m_S}$ and $E_{m_S}$. Reversing $\lambda$ amounts to reverse $m_S$. The branch $m_S$  at $\lambda>0$ has $2m_Sd$ modes that do not have a paring partner for $\lambda<0$. There are therefore $2m_Sd$ extra modes in the branch $m_S$ for $\lambda>0$ compared to $\lambda<0$ (in red). All the other number states can be paired with a number state of the space $E_{-m_S}$, or equivalently  $E_{m_S}$ with $\lambda<0$.}
\label{eq:unbalance}}
\label{fig:index}
\end{figure}

The expression \eqref{eq:flow_index} is similar to that of the index of Dirac operators introduced above. It is actually an index more generally associated to a class of operators called \textit{Fredolhm} operators \cite{APS1975, NakaharaBook}, which are invertible operators modulo a compact operator, and whose symbol is invertible. Here $\mathcal{D}^{2m_S}$ is Fredholm  since it is invertible between $E_{m_S}$ and $E_{-m_S}$, up to a finite number ($2m_Sd$) of elements, and we have
\begin{align}
 \Delta \mathcal{N}_{m_S}  = \text{ind} \mathcal{D}^{2m_S}
\label{eq:flow_indice}
\end{align}
with
\begin{align}
\text{ind} \mathcal{D}^{2m_S} = 2m_Sd
\label{eq:index_value}
\end{align}
in the present case. The values of this Fredholm index are found to correspond to those of the Chern numbers $\mathcal{C}_{m_S}$ calculated in \eqref{eq:chern_model} for the corresponding symbol Hamiltonian, for arbitrary spins $S$ and integer $d$, which thus verifies explicitly the Atiyah-Singer theorem in that case, or equivalently themonopole - spectral flow correspondence for this class of $h.\hat{S}$ models.\footnote{The $-$ sign difference between the analytical index and the Chern number comes from the convention $\mathcal{A}=+{\rm i} \bra{\psi}{\rm d} \ket{\psi}$ of the Berry connection. They would have the same sign if one took instead the condensed matter convention  $\mathcal{A}=-{\rm i} \bra{\psi}{\rm d} \ket{\psi}$.}
The  Berry-Chern  monopole-spectral flow correspondence then naturally follows from \eqref{eq:flow_indice} for this class of $h.\hat{S}$ models.

In the next section, we diagonalize ${\hat{\mathcal{H}}}_\text{op}$ to compute explicitly the spectral flows and their associated eigenstates.

\subsection{Eigenstates of ${\hat{\mathcal{H}}}_\text{op}$ and  verification of the monopole - spectral flow correspondence}
\label{sec:spectra_hop}

In this section, we aim at diagonalizing the operator Hamiltonian ${\hat{\mathcal{H}}}_\text{op}$ as ${\hat{\mathcal{H}}}_\text{op} \ket{\Psi} = \tilde{E} \ket{\Psi}$ and plotting its eigenvalues spectrum  $\tilde{E}$ as a function of $\lambda$ beyond the Dirac case  $S>1/2$ and $d=1$, in order to check explicitly the monopole - spectral-flow correspondence as predicted from the values of the Chern numbers \eqref{eq:chern_model} and the analytical index \eqref{eq:index_value}. For that purpose, we propose the following anzatz for the eigenstates
\begin{equation}
\text{Anzatz}\ A0: \qquad \ket{\Psi_0} = \sum_{{m_S}=-S}^{S} {\psi}_{m_S} \ket{(S+{m_S})d+n,{m_S}}
\label{eq:anzatz0}
\end{equation}
where again $\ket{n,{m_S}} \equiv \ket{n} \otimes \ket{{m_S}}$ where $\ket{n}$ is a number state, i.e. that satisfies ${\hat{a}^{\dagger}} {\hat{a}} \ket{n}=n\ket{n}$ and $\ket{{m_S}}$ is a spin state, i.e. it satisfies ${\hat{S}_z}\ket{m_S}=m_S\ket{m_S}$. For each value of $n$, the anzatz $A0$ is a superposition of  such $2S+1$ base states. If we wrote the spin operators in a matrix form, the anzatz $A0$ would take the form of a vector of number states that would read
\begin{align}
\left(
\begin{array}{cl}
\psi_S &\ket{(2S)d+n} \\
\psi_{S-1}& \ket{(2S-1)d+n}\\
\quad \vdots&  \vdots\\
\psi_{-S+1}&\ket{d+n}\\
\psi_{-S}&\ket{n}
\end{array}
\right) \ .
\end{align}
Consistently, one recovers the solution of the 2D Dirac and 3D Weyl operator Hamiltonians \eqref{eq:anzatz_dirac} for the branches $\pm$ when $S=1/2$ and $d=1$.

By applying $\ket{\Psi_0}$  to ${\hat{\mathcal{H}}}_\text{op}$ and then projecting onto the states $\bra{n+S-p,p}$, one obtained the eigenvalue equation
\begin{align}
\underbrace{\begin{pmatrix}
 S\, \lambda    &   \gamma_S  &    &     &   &    &   &   &   \\
\gamma_{S}      &                          &   &     &   &   &    &   &   \\
 &\ddots& \ddots& \ddots & & & & &\\
& &                                         \gamma_{m_S+1} & {m_S}\, \lambda    & \gamma_{m_S}  &\\
&&&\ddots&\ddots&\ddots&   \\
  & & & &  & & \gamma_{-S}   \\
  & & & &  & \gamma_{-S} & -S\, \lambda
\end{pmatrix}}_{H_n^{(2S+1)}}
\begin{pmatrix}
\psi_{S} \\
\psi_{S-1}\\
\vdots \\
\psi_{m_S}  \\
\vdots  \\
\psi_{-S+1}\\
\psi_{-S}
\end{pmatrix}
=\tilde{E}_n^{(m_S)}
\begin{pmatrix}
\psi_{S} \\
\psi_{S-1}\\
\vdots \\
\psi_{m_S}  \\
\vdots  \\
\psi_{-S+1}\\
\psi_{-S}
\end{pmatrix}
\label{eq:matrix}
\end{align}
where we have introduced the coefficients
\begin{align}
 \gamma_{m_S} =  \sqrt{S(S+1)-{m_S}({m_S}-1)} \prod_{j=0}^{d-1}\sqrt{(S+m)d+n-j}
\end{align}
which depends on $S, m_S, n$ and $d$. The upper index $m_S$ in $\tilde{E}_n^{(m_S)}$, which is the spin projection along $z$, can be seen as the branch index, which accordingly varies from $-S$ to $S$ by integer or half integer values. The lower index $n$ labels the discrete levels into each branch $m_S$. The matrix $H_n^{(2S+1)}$ to be diagonalized  has a size $(2S+1)\times (2S+1)$ and is defined for a fixed $n$. Its characteristic polynomial $P_n^{(2S+1)} \equiv \ensuremath{\operatorname{det}}(H_n^{(2S+1)}-\tilde{E}_n^{(m_S)}\, \text{Id})$ is therefore of degree $2S+1$ in $\tilde{E}^{(m_S)}$, and has real coefficients.\footnote{Do not confuse the \textit{degree} of the polynomial, which is simply the higher power of its variable, with the \textit{degree} $\deg h$ of the Hamiltonian map.} For each $n\in\mathbb{N}$, the secular equation $P_n^{(2S+1)}=0$ yields $2S+1$ real solutions, corresponding to a given mode $n$ in each of the energy branches.
The entire spectrum of ${\hat{\mathcal{H}}}_\text{op}$ is finally the list $\{\tilde{E}_n^{(m_S)}\}$ for all $n\in \mathbb{N}$ and for each branch $m_S$.
 Similarly, the coefficients $\psi_{m_S}$ should also depend on $n$. To simplify the notation, and because it is not crucial in the computation of the eigenvalues spectra,  we will abusively  write $\psi_{m_S}$ to refer to amplitudes with possibly different $n$.
The spectra obtained for different values of $S$ and $d$ with the anzatz $A0$ are shown in figure \ref{fig:sf_bulk}.

\begin{figure}[t]
  \centering
\includegraphics[width=14cm]{figure_spectrum_A0}
\caption{\footnotesize{A few examples of eigenvalue spectra $\tilde{E}$ of  ${\hat{\mathcal{H}}}_\text{op}$ obtained with the anzatz $A0$ in \eqref{eq:anzatz0} beyond the case  $S=1/2$ and $d=1$.}}
\label{fig:sf_bulk}
\end{figure}

None of these spectra manifest a spectral flow. The solutions of ${\hat{\mathcal{H}}}_\text{op}$ that do develop a spectral flow follow from different anzatz than $A0$  and can be constructed separately one by one. To do so, we notice that the lowest mode is reached for $n=0$, and has the following form when using the matrix representation for spins degrees of freedom

\begin{align}
\left(
\begin{array}{cl}
\psi_S &\ket{(2S)d} \\
\psi_{S-1}& \ket{(2S-1)d}\\
\quad \vdots&  \vdots\\
\psi_{-S+1}&\ket{d}\\
\psi_{-S}&\ket{0}
\end{array}
\right) \ .
\end{align}
The particularity of $\ket{\Psi_0}$ is that it decomposes as a superposition over the $2S+1$ spin states $\ket{m_S}$.
Actually, other eigenstates can be constructed and that decompose over $2S+1-p$ spin states $\ket{m_S}$ with $p>1$. The lower number state $\ket{0}$ being associated with the lower spin state $\ket{-S}$,  a first iteration for $p=1$ then consists in replacing this number state by $\ket{0} \rightarrow 0$ and reducing the other number states, associated to all the higher spin states $\ket{m_S>-S}$, to $n=-1$  as
\begin{align}
\left(
\begin{array}{cl}
\psi_S^{(1)} &\ket{(2S)d-1} \\
\psi_{S-1}^{(1)}& \ket{(2S-1)d-1}\\
\quad \vdots&  \vdots\\
\psi_{-S+1}^{(1)}&\ket{d-1}\\
0
\end{array}
\right) \ .
\end{align}
Such a form will be referred to as the anzatz $A1$ as it possesses $1$ zero  in its decomposition over spin states, and we can write
\begin{equation}
\text{Anzatz}\  A1: \qquad \ket{\Psi_{1}} = \sum_{m_S=-S+1}^{S} \psi_{m_S}^{(1)} \ket{(S+m_S)d-1,m_S} \ .
\end{equation}
This new Anzatz  yields a new eigenvalues problem where the matrix $H_{-1}^{2S}$ to be diagonalized  is obtained from $H_n^{(2S+1)}$ by removing the last ligne and the last column in \eqref{eq:matrix} and by substituting $n=-1$. It yields $2S$ energy branches  that connect the energy branches obtained previously with anzatz  $A0$. The procedure must be iterated $d$ times by lowering $n$ until the lower number state $\ket{0}$ is reached  as
\begin{align}
\left(
\begin{array}{cl}
\psi_S^{(1)} &\ket{(2S)d-d} \\
\psi_{S-1}^{(1)}& \ket{(2Sd-2d}\\
\quad \vdots&  \vdots\\
\psi_{-S+1}^{(1)}&\ket{0}\\
0
\end{array}
\right) \ .
\end{align}
giving at each iteration a different spectral flow between the same branches $m_S$ and $m_S'$. Once the spin component $\ket{-S+1}$ is finally associated with the fondamental $\ket{n=0}$ in number states, one iterates the anzatz  to look for solutions where both the spin components $m_S=-S$ and $m_S=-S+1$ are zero. More generally, the anzatz with $p$ zeros in the decompositions over the $p$ lower spin states reads
\begin{equation}
\text{Anzatz}\ Ap: \qquad \ket{\Psi_{p}} = \sum_{m=-S+p}^{S} \psi_{m_S}^{(p)} \ket{(S+m_S)d-p,m_S} \ .
\label{eq:anzatz-p}
\end{equation}

One applies this algorithm by keeping decreasing the number states until all the spin components are vanishing but the $m_S=S$ one as
\begin{align}
\left(
\begin{array}{c}
\ket{d-1} \\
0\\
  \vdots\\
0
\end{array}
\right) \ .
\end{align}
that is
\begin{equation}
\text{Anzatz}\ A(2S-1): \qquad \ket{\Psi_{2S-1}} =  \ket{d-1,S} \ .
\end{equation}
The number state $\ket{d-1}$ can still be decreased $d-1$ times to reach the final possibility with the fundamental $\ket{0}$ for the spin $m_S=S$ and $0$ amplitude for every other spin components.
Those $d-1$ 'sub-anzatz', so to speak, all yield the same eigenenergy spectrum with a non-dispersive branch $\tilde{E} = S \lambda$. This branch is common to all spectra for any values of $S$ and $d$, and is the only one existing in the Dirac case for $S=1/2$ and $d=1$.
\begin{figure}[t]
  \centering
\includegraphics[width=14cm]{figure_flow_spectral}
\caption{\footnotesize{A few examples of full eigenvalue spectra $\tilde{E}$ of  ${\hat{\mathcal{H}}}_\text{op}$ beyond the case  $S=1/2$ and $d=1$.  The black branches are obtained with the anzatz $0$ as in figure \ref{fig:sf_bulk} and the colored spectral flows are obtained by applying the successive anzatz $Ap$ \eqref{eq:anzatz-p}~: (red) $p=2S$, (blue) $p=2S-1$, (green) $p=2S-2$ and (yellow) $p=2S-3$.}}
\label{fig:sf_bulk_chiral}
\end{figure}

The complete spectra, that now display various spectral flows,  are shown in figure \ref{fig:sf_bulk_chiral}. Importantly, all those spectral flows satisfy the correspondence with the Berry-Chern monopole of the symbol Hamiltonian as $\mathcal{N}_{m_S} = -\mathcal{C}_{m_S}= 2m_Sd$,  where $\mathcal{N}_{m_S}$ is the number of modes that are gained by the band $m_S$ when sweeping $\lambda$ from $-\infty$ to $+\infty$.

After the canonical Dirac case ($S=1/2$ and $d=1$), the spin $1$ case with $d=1$ is the most encountered in the literature. It naturally appears  in the linear rotating shallow water model  describing rotating fluids such as oceans and atmosphere over large distances, with  a varying Coriolis parameter along the latitude that changes sign at the equator \cite{Matsuno1966,Delplace2017}. In that case, the Coriolis parameter plays the role of the mass term, and the spectral flow parameter $\lambda$ is the wave vector parallel to the equator. The spectral flow then corresponds to unidirectional Eastward equatorial waves  known as the Kelvin and the Yanai waves. The explicit mapping between ${\hat{\mathcal{H}}}_\text{op}$ and the shallow water model is shown in Appendix \ref{sec:appendix_SW}. A similar phenomenology appears with classical waves in various continuous systems such as electromagnetic waves in 2D gyrotropic media and plasmas, both in a varying magnetic field \cite{Marciani20, Parker_Alfven20}. In those physically very different examples, the spectral flows are interpreted as interface sates along a line in real space where a varying mass term (the Coriolis force or the magnetic field)  changes sign. As discussed at the beginning of this manuscript with the Weyl fermions, the same spectrum can emerge from a very different physical mechanism, that is the coupling of a spin 1 quantum particle in 3D with a magnetic field \cite{Bradlyn2016}. Other generalizations of Weyl fermions (i.e. beyond $S=1/2$ and $d=1$) have been discussed theoretically and experimentally in various materials \cite{Huang_PNAS2016, Ezawa2017, Hasan_Weyl2021}.

\section{Take Home message}

Spectral flows are ubiquitous in wave physics and quantum mechanics.
A simple method to look for spectral flows consists in searching for degeneracy points of the symbol Hamiltonian (i.e. of the dispersion relation), and then ''quantize'' the system by either varying in space a parameter controlling the gap amplitude (the mass term) such that it changes sign, or by applying a constant magnetic field if we one starts with massless quantum particles in 3D. The topological description of such spectral flows is based on a correspondence with a ''semiclassical'' description of the system, encoded through the symbol Hamiltonian, where the mass term is assumed to vary sufficiently slowly, that is over a length scale much larger than the wave length of reference. This limit corresponds to a description where the system can be decoupled into two systems with a fast and a slow dynamics, the latest being treated ''classically''. The topological indices of the symbol Hamiltonians (the first Chern numbers), describing the system in that limit, appears as charges associated to the degeneracies of its eigenvalues. They are related to analytical indices of the operator Hamiltonian through the Atiyah-Singer theorem, which are themselves a measure of the spectral flow in each branch. The existence of a spectral flow does not automatically yields the existence of an unidirectional mode, unless a gap separates different branches in the spectrum of the operator Hamiltonian. However, at least on the class of models discussed here, the modes participating to the spectral flow are 1) the most localized (in one direction) in space, and 2) those with the less nodes (the closest to the vacuum mode $\ket{0}$). They are in some sens localized at the ''boundary'' of the Fock space, if one interprets the states $\ket{n}$ as positions on a half infinite lattice labelled by $n$.

To summarize, the general structure of this theory is sketched in figure \ref{fig:summary}.

\begin{figure}[t]
  \centering
\includegraphics[width=15cm]{figure_summary}
\caption{\footnotesize{Quick overview of the general structure of the theory of the monopole-spectral flow correspondence. The blue lines represent an equality relation.}}
\label{fig:summary}
\end{figure}

\section*{Acknowledgement}
I am thankful to Cl\'ement Tauber, Antoine Venaille and Fr\'ed\'eric Faure for inspiring discussions.
This manuscrit was written partially from the lecture I gave at the \'Ecole Normale Sup\'erieure de Lyon from 2018, and in a much less detailed version at Les Houches and Cargèse summer schools in 2018 and 2019, and partially from my Habilitation thesis \cite{DelplaceHDR}.

\bibliographystyle{unsrt}

redpajama_set_name: RedPajamaArXiv