\section*{Data availability statement}
The datasets generated during and/or analysed during the current study are available in the Open Science Framework (OSF) repository,

\verb+https://osf.io/jnx8k/?view_only=16426f9a35404264badaaa93162060a7+









\renewcommand{\thefigure}{S\arabic{figure}}
\setcounter{figure}{0}
\renewcommand{\theequation}{S\arabic{equation}}
\setcounter{equation}{0}

\section*{Supplementary Information}
In this Supplementary Information, we provide more details on the behavior of the Hermitian part of the Hamiltonian of our system. We demonstrate how the Hamiltonian parameters are extracted from the experiment, including the non-Hermiticity (linear dichroism). Finally, we demonstrate that the winding of the real effective field constrains the possibility of a non-Hermitian topological transition with annihilation of the exceptional points: The annihilation is only possible if the total winding of the real effective field is zero.



\subsection{Topology of the Hermitian Hamiltonian: Theory}

In this section, we illustrate the topological transitions concerning the Hermitian part of the two-band Hamiltonian (1) of the main text describing the intersection of cavity modes with different numbers. This Hamiltonian is shown below for convenience, together with the definitions:
\begin{equation}
H_{\bf{k}}^{real} = \left( {\begin{array}{*{20}{c}}
{\frac{{E_H^{N+2} + E_V^{N}}}{2} + \frac{{{\hbar ^2}k_x^2}}{{2{m_x}}} + \frac{{{\hbar ^2}k_y^2}}{{2{m_y}}}}&{\Delta  - \beta '{k^2} - \beta {{({k_x} - i{k_y})}^2}}\\
{\Delta  - \beta '{k^2} - \beta {{({k_x} + i{k_y})}^2}}&{\frac{{E_H^{N+2} + E_V^{N}}}{2} + \frac{{{\hbar ^2}k_x^2}}{{2{m_x}}} + \frac{{{\hbar ^2}k_y^2}}{{2{m_y}}}}
\end{array}} \right),
\end{equation}
where $E_H^{N+2}$ and $m_H\sim N+2$ are the energy and mass, of the $N+2$th H-polarised mode and $E_V^{N}$ and $m_V\sim N$, are the energy and mass of the V-polarised mode number $N$. $k_x$, $k_y$ are the 2D wave vector components.
The spin-independent masses $m_x$ and $m_y$ are determined by the birefringence ($n_e$, $n_o$) and the angle of the optical axis $\theta$, according to the following expressions \cite{rechcinska2019engineering}: 
\begin{eqnarray}
    \frac{1}{m_x}&=&\frac{1}{2m'_0}\left(\frac{n_o^2\sin^2\theta+n_e^2\cos^2\theta+n_e^2}{n_e^2 n_o^2}\right),\\
    \frac{1}{m_y}&=&\frac{1}{2m'_0}\left(\frac{n_e^2+n_o^2}{n_e^2 n_o^2}\right),
\end{eqnarray}
where $m'_0$ is the mass at the central frequency (between the two modes in question).
$\beta$ is the magnitude of the TE-TM spin orbit coupling. $\beta'=\frac{\hbar^2(m_V-m_H)}{2m_Hm_V}$ and $\Delta=(E_H^{N+2}-E_V^{N})/2$.
This Hermitian Hamiltonian can be written as a linear combination of identity and Pauli matrices which defines a real effective magnetic field $\Omega_r$ acting on the polarization pseudo-spin. The two non-zero components of the field are $\Omega_r^x=\Delta-\beta'k^2-\beta(k_x^2-k_y^2)$ and  $\Omega_r^y=-2 \beta k_x k_y$.

The pseudospin of a given eigenstate is either aligned or anti-aligned with the effective field (so far as the Hamiltonian remains Hermitian). Experimentally, this pseudospin corresponds to the Stokes vector of light, calculated using the following expressions:
\begin{eqnarray}
    S^1&=&\frac{I_V-I_H}{I_H+I_V},\\
    S^2&=&\frac{I_D-I_A}{I_D+I_A},\nonumber\\
    S^3&=&\frac{I_R-I_L}{I_R+I_L},\nonumber
\end{eqnarray}
where the $I$s are the intensities in each of the polarizations, which can be measured experimentally or calculated theoretically as the absolute value squared of the eigenstate wavefunction in the corresponding basis. Each pseudospin component is an average value of the corresponding spin operator represented by a Pauli matrix. 

The case $\beta'=0$ has already been considered in \cite{Tercas2014,gianfrate2020measurement}. The two parabolas are split at $k=0$, with different possible signs of this splitting (band detuning).

\begin{figure}
    \centering
    \includegraphics[width=0.95\linewidth]{figS1.pdf}
    \caption{\textbf{Pseudospin textures for the Hermitian Hamiltonian}: \textbf{a}  trivial case $\beta<\beta'$, \textbf{b,c}  topologically non-trivial case $\beta>\beta'$ with positive and negative detunings.  The red marks shows the Dirac points of winding number +1 and the blue marks the Dirac points of winding number -1} 
    \label{texture}
\end{figure}

For positive detuning, the TE-TM splitting compensates this splitting along $k_x$ leading to a linear crossing. 
On the contrary, the TE-TM splitting is changing sign along $k_y$ and adds to the $k=0$ splitting. As a result, along $k_y$ the two parabola get away the one with respect to the other and do not cross. The resulting spectrum shows two tilted Dirac cones at $\pm k_{0x}$. The winding number of the effective field around these two points is the same (+1), as shown in Fig.~{\ref{texture}b} showing the distribution of the pseudospin (and effective field) texture in the reciprocal space. This winding is inherited from the winding number 2 of the TE and TM modes. The Berry topological charge is $+1/2$. This topological phase remains while $\beta'$ increases up to $\beta'=\beta$ where two new crossing points between the parabola appear at $k_{0y}$ which sets at infinity when the transition occurs. These crossing points correspond to tilted Dirac cones. The winding number of the in-plane pseudo-spin around these two points is (-1) and the Berry topological charge (-1/2), as shown on Figure~\ref{texture}-b. The sum of the topological charges of the 4 Dirac points is zero and the bands are overall topologically trivial which allows the non-Hermitian topological transition we explore in the main text.

For negative detuning (Fig.~\ref{texture}c), the winding of the two Dirac points does not change, but their positions in reciprocal space are rotated by $\pi/2$ and the particular field distribution changes from monopolar to Rashba.

\subsection{Topology of the Hermitian Hamiltonian: Experiment}


\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.8\linewidth]{figS2.pdf}
    \caption{\textbf{4 diabolical points case}. Experimental dispersion relation along \textbf{a} $k_x$ and \textbf{b} $k_y$ directions observed at 1.39\,V applied to the microcavity, corresponding to the case with total winding 0. Energy of the bands visible in horizontal (vertical) polarization is marked by green (purple) dashed line. \textbf{c} Experimental and \textbf{d} simulated by Berreman method pseudospin texture in $S_1$--$S_2$ plane of lower energy band.}
    \label{exp_N2N}
\end{figure}


\begin{figure}[tbp]
    \centering
    \includegraphics[width=0.8\linewidth]{figS3.pdf}
    \caption{\textbf{2 diabolical points case}. Experimental dispersion relation along \textbf{a} $k_x$ and \textbf{b} $k_y$ directions, observed at 11\,V applied to the microcavity, corresponding to the case with winding 2. Energy of the bands visible in horizontal (vertical) polarization is marked by green (purple) dashed line. \textbf{c} Experimental and \textbf{d} simulated by Berreman method pseudospin texture in $S_1$--$S_2$ plane of lower energy band.}
    \label{exp_NN}
\end{figure}


Cavities filled with birefringent nematic liquid crystal allows for realization of both topologically-distinct cases mentioned above: with non-trivial total winding 2 and  trivial total winding 0 (discussed in the main text), depending on voltage applied to the cell. For large enough energy detuning between the two modes of interest, the relative non-Hermiticity is weak enough so that the measured pseudo-spin is simply aligned with the real effective magnetic field (Hermitian limit). In such a case, the winding of the real field in reciprocal space can be determined experimentally through polarization resolved tomography.


Experiments described in this section were preformed on a sample analogous to the one described in the main text, but with wider LC layer of around 3.9\,$\mu$m.

Fig.\,\ref{exp_N2N} presents realization of the case with winding 0, discussed in detail in the main text. Fig.\,\ref{exp_N2N}a,b presents experimental dispersion relation along $k_x$ and $k_y$ directions observed with 1.39\,V applied to ITO electrodes. At this voltage, the numbers of the modes differ by 2, and both modes cross each other along both wave vector directions.  Experimentally determined pseudospin texture in $S_1$--$S_2$ plane of the lower energy band is plotted in Fig.\,\ref{exp_N2N}c. It  matches the Berreman matrix model presented in Fig.\,\ref{exp_N2N}d, showing 0 total winding.

With amplitude of external voltage of 11\,V modes with the same mode numbers are close to degenerate, as shown on experimental dispersion relations in Fig.\,\ref{exp_NN}. In this case modes crosses each other only along $k_y$ direction. Resulting  pseudospin texture in $S_1$--$S_2$ plane shows two 2D monopoles with total winding number of 2 in both experiment Fig.\,\ref{exp_NN}c and Berreman model simulations Fig.\,\ref{exp_NN}d.

The Dirac points appear along the Y axis because of the negative value of the detuning in the experiment, but their winding is that of the TE-TM field, as expected.



\subsection{Hamiltonian parameters}

\begin{figure}[tbp]
\includegraphics[width=0.95\linewidth]{figS4.pdf}
\caption{\textbf{Fit of the dispersions.} Fit of the polarization-resolved dispersions in the two perpendicular directions. The false color map shows the reflectance as a function of wave vector and energy. Fitting allows to obtain the spin-independent masses $m_x$ and $m_y$, as well as the polarization splittings $\beta$ and $\beta'$. \label{figDisp1}}
\end{figure}
The parameters of the Hermitian Hamiltonian (1) are extracted from the fit of the dispersion in the two perpendicular directions. The fits are shown in Fig.~\ref{figDisp1}.
In a given direction, the average of the two parabolicity coefficients (inverse masses) of the two polarizations, provides the value of $1/m_x$ and $1/m_y$. Then, the difference between the two parabolicites in the two polarizations is $\beta+\beta'$ for $Y$ direction and $\beta-\beta'$ for $X$ direction. Fitting both allows to extract $\beta$ and $\beta'$ separately.



\subsection{Non-Hermitian parameter}

In Fig.~\ref{DGamma}, we show an example of a spectrum of reflectance in two polarizations, clearly demonstrating different linewidths. The experimental points for each of the two polarizations are fitted with the Voigt function, as explained in Methods. We stress that while the non-Hermiticity of the Hamiltonian is a single constant, the linewidth at each particular wave vector is determined by the imaginary part of the corresponding eigenstate. Only at the positions of the Dirac points of the Hermitian Hamiltonian does the linewidth of the eigenstates correspond directly to the non-Hermiticity of the Hamiltonian $\delta \Gamma$. It is also the maximal possible imaginary part of the energy that could be observed for any wave vector. This is what we use to extract the non-Hermitian parameter $\delta\Gamma$.

\begin{figure}
    \centering
    \includegraphics[width=0.8\linewidth]{figS5.pdf}
    \caption{\textbf{Polarization-resolved reflectance spectrum.} The reflectance in polarizations H and V, measured at the position of the diabolical point and plotted as a function of energy. The difference of linewidths is approximately $13\%$.}
    \label{DGamma}
\end{figure}


\subsection{Hermitian winding and annihilation of exceptional points}

We consider a two-band system described by a Hamiltonian depending on a two-dimensional wave vector. The Hermitian part of the Hamiltonian is characterized by a presence of diabolical (Dirac) points. The two-band nature allows writing the Hamiltonian as a superposition of Pauli matrices. For many typical situations, such Hamiltonians can be written using only two Pauli matrices that we can choose to be $\sigma_x$ and $\sigma_y$ (massless Dirac Hamiltonian, Rashba and Dresselhaus spin-orbit couplings, TE-TM spin-orbit coupling for photons), whereas the third Pauli matrix $\sigma_z$, responsible for a symmetry breaking, opens the gap at the diabolical points (the mass term in the Dirac Hamiltonian or the Zeeman splitting for electrons and photonic modes). In this case, it is possible to characterize each of the diabolical points and the whole reciprocal space by winding numbers, which, in turn, determine the topology (the Chern number) of the bands once the gap is opened.

We then consider the evolution of the Hamiltonian with an addition of non-Hermiticity, described by a single constant parameter (independent of the wave vector). Without loss of generality, we consider the non-Hermitian contribution to be described by the Pauli matrix $\sigma_x$. In presence of the non-Hermiticity, each of the diabolical points gives rise to two exceptional points.

Thanks to the use of the Pauli matrices, it is possible to map the system's Hamiltonian to a spin-$1/2$ in a magnetic field, using the so-called pseudospin formalism, where the terms of the Hamiltonian are considered as an effective field. In the presence of a non-Hermiticity, this effective field has both real and imaginary components, $\bm{\Omega}=\bm{\Omega}'+i\bm{\Omega}''$. An exceptional point requires the real and imaginary parts of the effective field to be equal to each other in magnitude and perpendicular in direction: $\Omega'=\Omega''$ and $\bm{\Omega}'\perp \bm{\Omega}''$. 


\begin{figure}[tbp]
\includegraphics[width=1.05\linewidth]{figS6.pdf}
\caption{\textbf{Hermitian winding and annihilation of exceptional points.} (a) Non-zero total winding of the real field, the exceptional points do not annihilate. (b) Zero total winding of the real field, annihilation of exceptional points. The trajectories of the exceptional points are shown with red and blue arrows (corresponding to the winding of the exceptional points themselves). Crosses mark the annihilation points. \label{Sfig1}}
\end{figure}

Since in our case the direction of the imaginary field is fixed $\bm{\Omega}''=\Omega''\bm{e}_x$, and the real field is in the XY plane, the exceptional points can only be found where the x-component of the real field is zero $\bm{\Omega'}=\Omega'\bm{e}_y\Leftrightarrow\Omega'_x=0$. This condition constrains their evolution in the reciprocal space when the imaginary field is increased $\Omega''\to\infty$ (or the real field is decreased $\Omega'\to 0$). The location of the ensemble of the points $\Omega'_x=0$ is determined by the winding of the real field at the diabolical points. With the increase of $\Omega''$, the exceptional points are moving towards infinity. \textbf{If no points, where the condition for the real part $\Omega'_x=0$ is satisfied, can be found at infinity, then the exceptional points must annihilate for some finite critical value $\Omega''_c$.} Two cases are possible:
\begin{enumerate}
    \item The overall winding of the real field is a non-zero integer. In this case, the orientation angle $\theta_r$ of the real field at infinity changes at least once between $\theta_r=0$ and $\theta_r=2\pi$, necessarily going through $\pi/2$ and $3\pi/2$, where $\Omega'_x=0$. Therefore, the exceptional points can be found for any value of $\Omega''$ up to infinity, and thus they cannot annihilate (or if they annihilate accidentally, they must then reappear when $\Omega''$ is increased).
    \item The overall winding of the real field is zero. In this case, the orientation angle of the real field $\theta_r$ at infinity is constrained to an interval which may not include $\pi/2$ or $3\pi/2$. In this case, the exceptional points cannot be present at infinity, and therefore they must annihilate for a finite value of $\Omega''_c$.
\end{enumerate}

The two situations are illustrated in Fig.~\ref{Sfig1}. Panel (a) shows the case with total winding 2 with a dipolar texture typical for the TE-TM field. The trajectories of the exceptional points are shown with red and blue lines, which are defined by the condition $\Omega'_x=0$ (vertical orientation of the real field). At a large scale, the relative position of the two diabolical points does not play any role and the texture is given by $\theta_r=2\varphi$ (where $\varphi$ is the polar angle of the wave vector). The exceptional points at infinity  are  therefore necessarily located at $\varphi=\pm\pi/4, \pm3\pi/4$ (while initially they move from the diabolical points in the vertical direction).

The second situation is illustrated by Fig.~\ref{Sfig1}(b), where the overall winding of the real field is zero. The angle of the real field at infinity is $\theta_r=0$ and there can be no exceptional points. The condition $\Omega'_x=0$ (vertical real field) gives the trajectories of the exceptional points shown in the figure, with the annihilation points marked with  crosses.

Finally, the winding of the exceptional points themselves is shown with a color of the arrows (red and blue). 



The winding of an exceptional point is the winding of the phase of the complex eigenenergies. Indeed, the energy of a second-order exceptional point is given by
\begin{equation}
    E=\pm\sqrt{\alpha q}e^{iw\phi/2}
\end{equation}
where $w=\pm 1$ is precisely the winding number and $q$ is the absolute value of the parameter controlling the deviation from the exceptional point, for example, the wave vector (measured from the exceptional point). This eigenenergy can be obtained as a solution for different Hamiltonians, with different variation of real and imaginary effective fields.
Below, we show that when the imaginary effective field $\bm{\Omega}''$ is constant, the winding of the exceptional point can be easily determined from the texture of the real effective field $\bm{\Omega}'$. The most general form of a Hamiltonian of such type is
\begin{eqnarray}
    \hat{H}&=&\alpha q \left(\sigma_x\cos\varphi+iw\sigma_y\sin\varphi\right)\\
    &+&a(\left(\sigma_x\cos\theta+\sigma_y\sin\theta\right)\nonumber\\
    &+&ia\left(\sigma_x\cos(\theta\pm\pi/2)+\sigma_y\sin(\theta\pm\pi/2)\right)\nonumber
\end{eqnarray}
where the $a$ gives the values of the real and imaginary effective field oriented at an angle $\theta$ and $\theta+\pi/2$ respectively, equal and perpendicular at the exceptional point, $w$ is the winding of the parameter-dependent part of the real effective field and $\alpha$ is its strength. The eigenenergies of this Hamiltonian are given by 
\begin{equation}
    E\approx\pm\sqrt{q}\sqrt{2\alpha a}e^{\pm iw\varphi/2}e^{-i\theta/2}
\end{equation}
We see that the winding of an exceptional point is determined by the winding $w$ of the original diabolical point, modified by the relative angle of the real and imaginary fields at each exceptional point $\pm\pi/2$.




\end{document}