# PT and anti-PT symmetries for astrophysical waves

Armand Leclerc, Guillaume Laibe, Nicolas Perez

► To cite this version:

Armand Leclerc, Guillaume Laibe, Nicolas Perez. PT and anti-PT symmetries for astrophysical waves. Astronomy & Astrophysics - A&A, 2024, 689, pp.A237. 10.1051/0004-6361/202348710. hal-04617008

HAL Id: hal-04617008

https://hal.science/hal-04617008v1

Submitted on 17 Jan 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
PT and anti-PT symmetries for astrophysical waves

Armand Leclerc*, Guillaume Laibe, and Nicolas Perez

Univ Lyon, Univ Lyon1, ENS de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, 69230 Saint-Genis-Laval, France

Received 22 November 2023 / Accepted 28 May 2024

ABSTRACT

*Context.* Discrete symmetries have found numerous applications in photonics and quantum mechanics, but remain little studied in fluid mechanics, particularly in astrophysics.

**Aims.** We aim to show how *PT* and anti-*PT* symmetries determine the behaviour of linear perturbations in a wide class of astrophysical problems. They set the location of 'exceptional points' in the parameter space and the associated transitions to instability, and are associated with the conservation of quadratic quantities that can be determined explicitly.

*Methods.* We study several classical local problems: the gravitational instability of isothermal spheres and thin discs, the Schwarzschild instability, the Rayleigh-Bénard instability and acoustic waves in dust-gas mixtures. We calculate the locations and the order of the exceptional points using the resultant of two univariate polynomials, as well as the conserved quantities in the different regions of the parameter space using Krein theory.

*Results.* All problems studied here exhibit discrete symmetries, even though Hermiticity is broken by different physical processes (self-gravity, buoyancy, diffusion, and drag). This analysis provides genuine explanations for certain instabilities, and for the existence of regions in the parameter space where waves do not propagate. Those two aspects correspond to regions where *PT* and anti-*PT* symmetries are broken respectively. Not all instabilities are associated to symmetry breaking (e.g. the Rayleigh-Benard instability).

**Key words.** diffusion – hydrodynamics – instabilities – waves

# 1. Introduction

Many astrophysical systems, such as stars and discs, are often treated as fluids (e.g. Pringle & King 2007; Armitage 2010). The properties and the stability of the corresponding linear modes are usually analysed with methods that, to our knowledge, do not fully exploit the symmetries of the perturbed system with respect to parity *P* and time *T*.

Some of the most striking properties of *PT*-symmetric systems emerged in the quantum physics community with the seminal paper by Bender & Boettcher (1998), the first non-Hermitian quantum Hamiltonian situation with a real spectrum, which paved the way to open quantum systems. This framework has been extended to situations exhibiting a spontaneous *PT*-symmetry-breaking phase transition (Bender et al. 1999, 2002). Generalisations to quantum field theory were developed by Bender et al. (2004) and considerations of the observables of such systems were put forward by Mostafazadeh & Batal (2004). *PT* symmetry has been a particularly powerful tool in optics and photonics, where it has led to the development of novel designs with special properties. The analogy with quantum mechanics comes from the fact that beam dynamics is governed by a Schrödinger equation, where the optical index behaves as a complex potential, which can be designed to satisfy *PT* symmetry (El-Ganainy et al. 2007; Makris et al. 2008). Klaiman et al. (2008) identified a spontaneous breaking of *PT* symmetry in a waveguide. Musslimani et al. (2008) investigated the non-linear propagation of optical solitons in *PT*-symmetric refraction indexes. (Guo et al. 2009) and (Rüter et al. 2010) made the first observations of *PT* symmetry in optics with a

complex index. A decade of developments of *PT* symmetry in photonics lead to a multitude of applications (e.g. Feng et al. 2017; El-Ganainy et al. 2018; Özdemir et al. 2019). Feng et al. (2014) used spontaneous breaking of *PT* symmetry to design a pure single-mode laser. Peng et al. (2014) propose a design for a low-power optical diode, where the non-linear regime of the *PT* symmetry induces non-reciprocal propagation. Lin et al. (2011), Regensburger et al. (2012) and Feng et al. (2013) used *PT* symmetry in periodic systems (through gratings, periodic crystals, and lattices) to design unidirectional invisibility. *PT*-symmetric degeneracies are now used for their high resonant sensitivity in optical cavities (Hodaei et al. 2017). Anti-*PT* symmetry has also found application in novel optical designs (e.g. Ge & Türeci 2013; Zhang et al. 2020a,b). The topological properties of *PT*-symmetric systems have recently been investigated in the context of fluid dynamics (Fu & Qin 2024). *PT* and anti-*PT* symmetries have only recently been taken into account in fluid systems to directly constrain fluid flows and instabilities by spontaneous symmetry breaking (Qin et al. 2019; Fu & Qin 2020; David et al. 2022). The *PT* symmetry framework has been shown to be particularly adapted to find instabilities in high-dimensional parameter space (e.g. the drift wave instabilities in tokamaks; Qin et al. 2021). It has thus been demonstrated unequivocally that *PT*-symmetric systems exhibit particular properties.

In this study, we aim to show that *PT* and anti-*PT* symmetries also control some of the properties of the linear modes of astrophysical systems, such as instabilities, the absence of propagation, stability exchange, and the conservation of quadratic quantities, all acting in physical processes of different origins. We first define *PT* and anti-*PT* symmetries in Sect. 2. We then demonstrate how they condition the outcome of the linear

<footnote>* Corresponding author; armand.leclerc@ens-lyon.fr</footnote>
analysis of four canonical astrophysical systems, namely the self-gravity instability in discs (Sect. 3), wave propagation in gas-dust mixtures (Sect. 6), and buoyancy instability in stellar (Sect. 4) and planetary interiors (Sect. 5). We then discuss these analyses in a broader context of symmetric systems in Sect. 7.

## 2. Definitions

### 2.1. PT symmetry

The evolution of a monochromatic linear perturbation to the steady state of a fluid can be written under the form of an eigenvalue problem $H(a)X = ωX$, where the operator $H$ that describes the evolution of the perturbation depends continuously on some real parameter $a$. The complex eigenvalues are generically denoted by the frequency $ω$, which can be a complex number $a$ priori. This problem is said to be $PT$ symmetric if there exists a unitary operator $U$ such that

$$U H U^{-1} = H^*, \qquad (1)$$

where $^*$ stands for complex conjugation. In this case, the eigenvalues form pairs: if $ω$ is an eigenvalue of the problem, so is $ω^*$. Indeed, $HX = ωX$ implies $H(U^{-1}X^*) = ω^*(U^{-1}X^*)$. In $PT$-symmetric systems, unstable modes are therefore accompanied by damped modes (imaginary parts of opposite sign). When $ω ≠ ω^*$, $PT$ symmetry is said to be spontaneously broken by $H$. By defining the symmetry operator $PT ≡ UΘ$, where $Θ$ denotes complex conjugation, one has $[PT, H] = 0$. The two operators commute, but do not share $X$ as an eigenvector for non-real eigenvalues.

Mostafazadeh (2002b,a) and Zhang et al. (2020c) showed that $PT$ symmetry is equivalent to another symmetry, pseudo-Hermiticity, for finite-dimensional problems. A problem described by the operator $H$ is pseudo-Hermitian if there exists a Hermitian operator $V$ such that

$$V H V^{-1} = H^†, \qquad (2)$$

where $^†$ denotes the Hermitian conjugate with respect to the appropriate Hermitian scalar product of the problem, which we will write in the following, such that $Y_1 \cdot Y_2 = Y_1^* Y_2$.

Let $X_1$ and $X_2$ be two eigenvectors of $H$ satisfying Eq. (2). The eigenvalues $ω_1$ and $ω_2$ may take the same value: they degenerate. When the eigenvectors also become identical, the degeneracy is called an exceptional point (EP). The linear stability of these systems is described by Krein theory (Krečín 1950; Kirillov 2013): in such systems, the transition between stability and instability necessarily involves the exceptional degeneracy of stable eigenvalues. This phenomenon is known as a Krein collision. Before the collision, the system is stable, $ω_1$ and $ω_2$ are real, $PT$ symmetry is said to be unbroken, and the Krein quantities defined as $X_{1,2} \cdot VX_{1,2}$ are non-zero and of opposite sign. After the collision, the system is unstable, the complex frequencies $ω_1 = ω_2^*$ are complex conjugates of each other, and the two Krein quantities are necessarily zero. The Krein quantity $Y \cdot VY$ is a constant of motion of a general solution $Y(t, a)$ of the evolution equation $i\partial_t Y = H(a)Y$. Indeed,

$$i\partial_t (Y \cdot VY) = Y \cdot (-H^†V)Y + Y \cdot VHY \qquad (3)$$

$$= Y \cdot (VH - H^†V)Y \qquad (4)$$

$$= 0. \qquad (5)$$

Hence,

$$Y \cdot VY = \text{cst.} \qquad (6)$$

Moreover, for a Fourier mode $Y = e^{-iωt}X$, this conserved quantity is exactly zero for unstable and damped modes, and non-zero for propagative modes. This is shown by expressing

$$X \cdot VHX = X \cdot H^†VX = ω^* X \cdot VX \qquad (7)$$

$$= ω X \cdot VX. \qquad (8)$$

According to the identity above, a situation where $ω ≠ ω^*$ implies that $X \cdot VX = 0$.

### 2.2. Anti-PT symmetry

The eigenvalue problem $H(a)X = ωX$ is said to be anti-$PT$-symmetric, or alternatively $CP$-symmetric, if there exists a unitary operator $\tilde{U}$, such that

$$\tilde{U} H \tilde{U}^{-1} = -H^*. \qquad (9)$$

In this case, the eigenvalues come in pairs: $ω$ and $-ω^*$. They have equal imaginary parts, and opposite real parts and describe counter-propagating waves. In particular, $H$ being anti-$PT$ symmetric is equivalent to $iH$ being $PT$ symmetric. The same equivalence between $PT$ symmetry and pseudo-Hermiticity then applies, and guarantees that an anti-$PT$-symmetric operator $H$ is also pseudo-chiral: there exist an Hermitian operator $\tilde{V}$, such that

$$\tilde{V} H \tilde{V}^{-1} = -H^†. \qquad (10)$$

The quantity $X \cdot \tilde{V}X$ is not strictly a Krein quantity, and is not in general a constant of motion of a solution of $i\partial_t X = HX$. However, for a Fourier mode, $Y = e^{-iωt}X$, one has

$$Y \cdot \tilde{V}Y = e^{2\text{Im}(ω)t} X \cdot \tilde{V}X, \qquad (11)$$

and is constant to zero for modes with $\text{Re}(ω) ≠ 0$; that is, for modes that spontaneously break the anti-$PT$ symmetry. Indeed,

$$X \cdot \tilde{V}HX = X \cdot (-H^†\tilde{V})X = -ω^* X \cdot \tilde{V}X \qquad (12)$$

$$= ω X \cdot \tilde{V}X. \qquad (13)$$

In summary, whenever $ω ≠ -ω^*$, $X \cdot \tilde{V}X = 0$, the Krein quantity is zero at all times. An analogue of a Krein collision may occur. When the anti-$PT$ symmetry is unbroken, their exists a pair of eigenvalues of $H$, denoted $ω_1$ and $ω_2$, which are purely imaginary.

Let us consider the situation in which the parameter $a$ passes through an EP, so that after the degeneracy we have $ω_1 = -ω_2^*$ with non-zero real parts. Then, the anti-$PT$ symmetry is spontaneously broken. This phenomenon corresponds to a Krein collision for the operator $iH$. It follows that the quantity $X \cdot \tilde{V}X$ is non-zero and not conserved for non-propagative modes, and strictly zero for propagating modes for all times.

Systems that are both $PT$ and anti-$PT$ symmetric are referred to as bisymmetric. The eigenvalues of bisymmetric systems take the form of $(ω, -ω)$ pairs of real or purely imaginary frequencies. When $ω$ is real, the associated eigenvector spontaneously breaks the anti-$PT$ symmetry. Conversely, when $ω$ is imaginary, the associated eigenvector spontaneously breaks $PT$ symmetry. Thus, one of those two symmetries is necessarily broken for any value of the parameters $a$, and either $X \cdot VX = 0$ (broken $PT$) or $X \cdot \tilde{V}X = 0$ (broken anti-$PT$). In bisymmetric systems, $X \cdot VX$ is always a constant of motion, while $X \cdot \tilde{V}X$ is constant only when it is zero, that is, in the spontaneously broken anti-$PT$ symmetry phase.
### 3. Gravitational instabilities

#### 3.1. Jeans instability

The stability of self-gravitating objects plays a central role in the formation of structures throughout the Universe. The simplest problem to study would be the Jeans instability of a collapsing sphere (Jeans 1902). We consider local radial adiabatic perturbations of a 3D homogeneous self-gravitating, non-rotating, inviscid sphere of constant density $\rho_0$ and sound speed $c_s$. Linear perturbations of mass and momentum conservation combined with a Poisson equation form a 4 × 4 system that can be written after Fourier transforms in time and space:

$$ \omega \begin{pmatrix} \rho' / \rho_0 \\ \mathbf{v}' / c_s \end{pmatrix} = H_J \begin{pmatrix} \rho' / \rho_0 \\ \mathbf{v}' / c_s \end{pmatrix} \equiv \begin{pmatrix} 0 & c_s \mathbf{k}^{\mathrm{T}} \\ \left(1 - \frac{4\pi G \rho_0}{c_s^2 |k|^2}\right) c_s \mathbf{k} & 0 \end{pmatrix} \begin{pmatrix} \rho' / \rho_0 \\ \mathbf{v}' / c_s \end{pmatrix}. \quad (14) $$

This problem is invariant under both parity and time reversal. It is therefore trivially *PT* symmetric, as proven by the statement $H_J = H_J^*$. As such, $H_J$ is pseudo-Hermitian, with the symmetry operator $V = \text{diag}((1 - \frac{4\pi G \rho_0}{c_s^2 |k|^2}), 1, 1, 1)$. The Krein quantity associated to the Krein collision of the Jeans instability is then $X \cdot VX = (1 - \frac{4\pi G \rho_0}{c_s^2 |k|^2}) |\rho'|^2 + |\mathbf{v}'|^2$. This quantity is non-zero for propagating waves, and zero for the unstable modes. Complementarily, $H_J$ is anti-*PT* symmetric and pseudo-chiral; indeed, it satisfies Eq. (10) for $\tilde{V} = \text{diag}(-\left(1 - \frac{4\pi G \rho_0}{c_s^2 |k|^2}\right), 1, 1, 1)$. The quantity $X \cdot \tilde{V}X = -\left(1 - \frac{4\pi G \rho_0}{c_s^2 |k|^2}\right) |\rho'|^2 + |\mathbf{v}'|^2$ is then zero for propagative modes and non-constant and non-zero for unstable modes. The Krein collision occurs at an EP for which $H$ cannot be diagonalised. This EP corresponds to marginal stability, which is reached at the Jeans wavenumber $|k| = k_J \equiv \sqrt{4\pi G \rho_0 / c_s}$. Longer wavelengths are unstable, and shorter wavelengths are propagating.

#### 3.2. Toomre instability

More interesting is the case of a self-gravitating disc (Goodman & Narayan 1988; Goodman 2003; Bertin & Lodato 1999). In this system, *P* and *T* are broken individually, but the overall system of linear perturbations is still *PT* symmetric. Stabilisation of large and small scales by rotation and pressure can be sufficient to render astrophysical discs robust against gravitational collapse (Toomre 1964). Here we discuss the simplest case of axisymmetric perturbations of short radial wavelengths evolving in a razor-thin Keplerian disc. The stability of linear perturbations is given by the celebrated Toomre criterion: $Q \equiv \frac{c_s k}{\pi G \Sigma_0} > 1$, where $k$ denotes the epicyclic frequency of the disc and $\Sigma_0$ its surface density. The razor-thin disc is integrated vertically (Armitage 2010), to give the following set of radial linear perturbations:

$$ H_{\mathrm{sgd}} \begin{pmatrix} h' \\ v'_{r} \\ v'_{\phi} \end{pmatrix} = \begin{pmatrix} 0 & c_s k & 0 \\ -\frac{k}{c_s} \left( \frac{2}{Q|k|} - 1 \right) & 0 & 2i \\ 0 & -\frac{i}{2} & 0 \end{pmatrix} \begin{pmatrix} h' \\ v'_{r} \\ v'_{\phi} \end{pmatrix} = \omega \begin{pmatrix} h' \\ v'_{r} \\ v'_{\phi} \end{pmatrix}, \quad (15) $$

where we use perturbations of the enthalpy $h'$ and the two horizontal velocities $v'_r$ and $v'_\phi$ as variables, and the orbital time $\Omega^{-1}$ and the pressure length $c_s \Omega^{-1}$ as units of time and length.

With $U = \text{diag}(1, 1, -1)$, this problem satisfies Eq. (1) and is therefore *PT* symmetric. Rotation breaks the reflection symmetry in the azimuthal direction and the time reversal symmetry, but not the combination of both. The system is expected to be pseudo-Hermitian, which can be shown explicitly. With

$$ V = \begin{pmatrix} \frac{1}{c_s^2} & 0 & -\frac{4ik}{c_s Q |k|} \\ 0 & 1 & 0 \\ \frac{4ik}{c_s Q |k|} & 0 & 4\left(\frac{2|k|}{Q} + 1\right) \end{pmatrix}, \quad (16) $$

the system satisfies Eq. (2). The Krein quantity associated with this *PT* symmetry provides the following energy partition, which applies to the unstable mode:

$$ X \cdot VX = 0 = \frac{|h'|^2}{c_s^2} + |v'_r|^2 + 4|v'_\phi|^2 \left(1 + \frac{2\pi G \Sigma_0}{\Omega^2} |k|\right) - \frac{8\pi G \Sigma_0}{c_s^2 \Omega} \text{sgn}(k) \text{Im}(hv'_\phi^*), \quad (17) $$

where $k$ denotes the wavenumber in dimensional form. We illustrate this result by numerically computing the value of this quantity through a Krein collision, as $H_{\mathrm{sgd}}$ has degenerated eigenvalues when $Q = \frac{2|k|}{1+k^2}$. These are the positions of EPs for $H_{\mathrm{sgd}}$, which corresponds to positions of a Krein collision. The curve of EPs separates two regions in $(Q, k)$ space: one is the stable region, and the other is the unstable region. Figure 1 shows the EP curve, as well as the Krein collision that occurs when crossing this curve.

In addition to being *PT* symmetric, $H_{\mathrm{sgd}}$ is also anti-*PT* symmetric, because

$$ \tilde{U} H_{\mathrm{sgd}} \tilde{U}^{-1} = -H_{\mathrm{sgd}}^*, \quad (18) $$

with $\tilde{U} = \text{diag}(1, -1, -1)$. The anti-*PT* symmetry is complementary to the *PT* symmetry, and the combination of the two symmetries implies that the eigenvalues always come in pairs: they are always opposite $(\omega, -\omega)$, either real or purely imaginary. They are never general complex numbers, which is a generic property for bisymmetric systems, as one of the two symmetries is necessarily spontaneously broken.

One has

$$ \tilde{V} H_{\mathrm{sgd}} \tilde{V}^{-1} = -H_{\mathrm{sgd}}^{\dagger}, \quad (19) $$

$$ \tilde{V} = \begin{pmatrix} \frac{1}{c_s^2} & 0 & -\frac{4ik(Q|k|-1)}{cQ|k|} \\ 0 & 1 & 0 \\ \frac{4ik(Q|k|-1)}{cQ|k|} & 0 & -\frac{8k^2}{Q|k|} + 8k^2 - 4 \end{pmatrix}. \quad (20) $$

Equation (20) explicitly shows that $H_{\mathrm{sgd}}$ is pseudo-chiral as expected, and Krein theory then applies. When $iH_{\mathrm{sgd}}$ has real eigenvalues, anti-*PT* symmetry is unbroken, and the Krein quantities $X \cdot \tilde{V}X$ are non-zero. This holds when *PT* symmetry is spontaneously broken, that is, in the unstable region of parameters, where the first Krein quantity $X \cdot VX$ is zero. Conversely, in the stable range, anti-*PT* symmetry is broken and *PT* symmetry is unbroken. In this case, $X \cdot \tilde{V}X$ is zero and $X \cdot VX$ is not equal to zero. The Krein quantity $X \cdot \tilde{V}X$ associated with the anti-*PT* symmetry gives the following energy distribution for propagating waves in dimensional form:

$$ X \cdot \tilde{V}X = \frac{|h'|^2}{c_s^2} + |v'_r|^2 + |v'_\phi|^2 \left( -4 + 8 \left( \frac{c_s^2 k^2}{\Omega^2} - \frac{\pi G \Sigma_0}{\Omega^2} |k| \right) \right) - \frac{8\pi G \Sigma_0}{c_s^2 \Omega} \text{sgn}(k) \left( \frac{c_s^2}{\pi G \Sigma_0} |k| - 1 \right) \text{Im}(hv'_\phi^*). \quad (21) $$
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** PT-symmetric properties of the Toomre instability. Top panel: Curve of EPs corresponding to marginal stability. A maximum is reached for $Q = 1$ as expected. Bottom panel: Krein collision of the stability of the Toomre problem when varying the parameters along the black line of the top panel. Depending on the value of $Q$, exactly one of the two Krein quantities identified is zero, depending on which symmetry is broken by the perturbation. The curve in green is $X \cdot VX$, orange is $X \cdot \tilde{V}X$, solid blue is $\mathrm{Re}(\omega)$, and dotted blue is $\mathrm{Im}(\omega)$, for a destabilising mode.</figcaption></figure>

Hence, for each mode, either $X \cdot VX$ or $X \cdot \tilde{V}X$ is necessarily zero: the former is zero for unstable modes and the latter is zero for propagating waves. Both are exactly zero at the degeneracy; that is, for $\omega = 0$, which is true for $Q = \frac{2|k|}{1+k^2}$. Figure 1 shows this exchange in non-zero quantities when crossing the Krein collision. Krein quantities simply provide energy partitions of the system in both regimes.

The matrices $V$ and $\tilde{V}$ identified above are particular cases of more general families of matrices:

$$V = \begin{pmatrix} \frac{2c_s(-2k^2 + (-1+k^2)Q|k|)a + k(-2+Q|k|)c}{4c_s^2 k Q |k|} & b \frac{(kQ-2\operatorname{sgnk})}{2c_s Q} & ia \\ b \frac{(kQ-2\operatorname{sgnk})}{2c_s Q} & \frac{1}{4}(2c_s ka + c) & ib \\ -ia & -ib & c \end{pmatrix}, \quad (22)$$

where $a, b$, and $c$ are any real numbers such that $\det W \neq 0$. The matrix of Eq. (16) is recovered with $a = -\frac{4k}{c_s Q|k|}$, $b = 0$, and

$$c = \frac{8|k|}{Q} + 2. \text{ Similarly,}$$

$$\tilde{V} = \begin{pmatrix} \frac{2c_s(-2k^2 + (-1+k^2)Q|k|)a + k(-2+Q|k|)c}{4c_s^2 k Q |k|} & b \frac{i(kQ-2\operatorname{sgnk})}{2c_s Q} & ia \\ -b \frac{i(kQ-2\operatorname{sgnk})}{2c_s Q} & -\frac{1}{4}(2c_s ka + c) & b \\ -ia & b & c \end{pmatrix}. \quad (23)$$

The matrix of Eq. (20) is recovered with $a = -\frac{4k(Q|k|-1)}{c_s Qk}$, $b = 0$, and $c = -\frac{8}{Q}|k| + 8k^2 - 4$. The same holds for the different problems addressed in this study.

## 4. Buoyancy instability in stars

We now consider a non-rotating star in static equilibrium, balanced by thermal pressure and self-gravity. The steady state is spherically symmetric, and the density is stratified and decreases towards the surface. This equilibrium can be unstable, as a perturbation can be amplified by buoyancy if the square of the buoyancy frequency $N^2$ is negative, as given by the Schwarzschild criterion (Schwarzschild 1906).

This problem was revisited by Leclerc et al. (2023) in inhomogeneous media, where the system proved to be pseudo-Hermitian and pseudo-chiral. Pseudo-Hermiticity was used to determine the Krein quantity. We now complete the picture by determining the Krein quantity associated to pseudo-chirality, and show the Krein collision in this bisymmetric problem.

Neglecting the self-gravity of the perturbations (Cowling’s approximation Cowling 1941), we start from Eq. (5) of Leclerc et al. (2023), which is the eigenvalue equation,

$$\omega X = H_b X, \quad (24)$$

$$H_b = \begin{pmatrix} 0 & 0 & 0 & L_\ell \\ 0 & 0 & |N| & K_r - iS \\ 0 & -|N| & 0 & 0 \\ L_\ell & K_r + iS & 0 & 0 \end{pmatrix}, \quad (25)$$

where $\omega$ is the complex eigenfrequency of the perturbation, and $X = (v \ w \ \Theta \ p)^\mathrm{T}$ contains the perturbation’s horizontal velocity, radial velocity, entropy, and pressure after appropriate rescaling. $N^2 = -g \frac{d\ln \rho_0}{dr} - \frac{g^2}{c_s^2}$ is the square of the buoyancy frequency, which is negative here, $S = \frac{c_s}{2g} (N^2 - \frac{g^2}{c_s^2}) - \frac{1}{2} \frac{dc_s}{dr} + \frac{c_s}{r}$ is another characteristic frequency called the buoyant-acoustic frequency, which quantifies the coupling between g-modes and p-modes in asteroseismic problems, $L_\ell = \frac{c_s}{r} \sqrt{\ell(\ell+1)}$ is the Lamb frequency, and $K_r$ represents the local radial wavenumber of the wave (see Appendix E of Leclerc et al. 2023 for details).

This problem is both pseudo-Hermitian and pseudo-chiral, and is associated with the two matrices:

$$V = \mathrm{diag}(1, \ 1, -1, \ 1), \quad (26)$$

$$\tilde{V} = \mathrm{diag}(1, \ 1, \ 1, -1). \quad (27)$$

Krein theory can therefore be applied to determine how discrete symmetries constrain the partition of energy. Here, $X \cdot X$ is generally not a conserved quantity, as $H_b$ is not Hermitian; it grows exponentially for unstable modes and is conserved only for stable modes. However, $X \cdot VX$ is a constant of motion for any solution of $\partial_t X = H_b X$, as shown in Sect. 2.1. Moreover, this constant is exactly zero only for unstable modes. In contrast, $X \cdot \tilde{V}X$ also grows exponentially for unstable modes, but is constantly equal
to zero for stable modes. As such,

$$X \cdot X = |v|^2 + |w|^2 + |\Theta|^2 + |p|^2 \propto e^{2\mathrm{Im}(\omega)t}, \quad (28)$$

$$X \cdot VX = |v|^2 + |w|^2 - |\Theta|^2 + |p|^2 = \mathrm{cst} \, \delta_{\mathrm{Im}(\omega),0}, \quad (29)$$

$$X \cdot \tilde{V}X = |v|^2 + |w|^2 + |\Theta|^2 - |p|^2 \propto e^{2\mathrm{Im}(\omega)t} \delta_{\mathrm{Re}(\omega),0}, \quad (30)$$

where $\delta$ is the Kronecker delta.

For the radial modes $L_l = 0$, the modes degenerate over the ring parameterised by $K_r^2 + S^2 = -N^2 > 0$. These parameter values are EPs, as $H_b$ can no longer be diagonalised there. The ring separates two ranges of parameters, the stable and the unstable part. Figure 2 shows the Krein quantities, over the ring of EPs in which the collision takes place. The Krein collision behaves as expected: inside the ring, the system is unstable, and $X \cdot VX$ is zero and $X \cdot \tilde{V}X$ is non-zero. Outside the ring, the modes are stable and propagate as sound waves, and $X \cdot VX$ is non-zero and $X \cdot \tilde{V}X$ is zero.

This structure is equivalent to the bisymmetric problem of self-gravitating discs discussed in Sect. 3, and shown on Fig. 1, even if the physical situation is very different. The approach followed in this problem is the reverse of that followed in Section 3 for the Toomre problem: it is easier to first establish that the system is pseudo-Hermitian and pseudo-chiral, as the matrices $V$ and $\tilde{V}$ have $\pm 1$ on the diagonal.

## 5. Rayleigh-Bénard instability

Not all instabilities correspond to spontaneous breaking of symmetries. The Rayleigh-Bénard instability (Bénard 1900; Rayleigh 1916) is such a case. Interiors of rocky planets are often treated as incompressible with both significant viscous and thermal diffusion (e.g. Bergé & Dubois 1984; Bodenschatz et al. 2000; Brandenburg 2021). This results in a regime of buoyancy instability that is controlled by diffusion effects and differs fundamentally from that discussed in Sect. 4 for adiabatic perturbations in compressible stars or planetary atmospheres. In the Rayleigh-Bénard regime, the instability criterion is given by the value of the Rayleigh number $Ra = g\alpha\beta L^4 / \nu\kappa$, which needs to be greater than some critical value of order $10^2-10^3$ depending on the boundary conditions. This famous problem admits analytical solutions for a homogeneous background configuration with rigid-lid isothermal boundary conditions, which provides us with an illustrative example of symmetries involving diffusive effects.

The basic set of equations consists of the incompressibility condition, a Navier-Stokes equation, and thermal heat diffusion for a 2D $(x, z)$ fluid (see Appendix A). The $x$ direction is invariant by translation and infinite, and the $z$ direction is of length $L$, along which a constant temperature gradient $\partial_z T_0 = -\beta$ is directed. With dimensionless variables, the problem for Fourier modes $\exp(-i\omega t + ik_z z + ik_x x)$ can be converted into the symmetrical form:

$$\omega \begin{pmatrix} \tilde{v}_z \\ \tilde{T} \end{pmatrix} = H_{\mathrm{RB}} \begin{pmatrix} \tilde{v}_z \\ \tilde{T} \end{pmatrix}, \quad (31)$$

$$H_{\mathrm{RB}} = -i \begin{pmatrix} k^2 & -\frac{k_x}{k} \left(\frac{\mathrm{Ra}}{\mathrm{Pr}}\right)^{1/2} \\ -\frac{k_x}{k} \left(\frac{\mathrm{Ra}}{\mathrm{Pr}}\right)^{1/2} & k^2/\mathrm{Pr} \end{pmatrix}, \quad (32)$$

where $\mathrm{Pr} = \nu/\kappa$ is the Prandtl number and $k = \sqrt{k_x^2 + k_z^2}$.

This problem is anti-$\mathcal{PT}$ symmetric and pseudo-chiral, as one reads that

$$H_{\mathrm{RB}} = -H_{\mathrm{RB}}^* = -H_{\mathrm{RB}}^{\dagger}, \quad (33)$$

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>Fig. 2. Same as Fig. 1 but for the Schwarzschild instability. Non-hermiticity is broken in this case by buoyancy instead of self-gravity.</figcaption></figure>

for $Ra \ge 0$. In this case, the anti-$\mathcal{PT}$ symmetry is protected and cannot be spontaneously broken as $H_{\mathrm{RB}}$ is anti-Hermitian (i.e. $iH_{\mathrm{RB}}$ is Hermitian). Indeed, the spectral theorem guarantees that the eigenvalues are purely imaginary and the eigenvectors are orthogonal to each other.

The location of EPs is given by

$$Ra = \left(1 - \frac{(1 + \mathrm{Pr})^2}{4\mathrm{Pr}}\right) \frac{(k_x^2 + k_z^2)^3}{k_x^2}, \quad (34)$$

which requires that $Ra$ be negative. This condition can only be fulfilled for an inverted temperature gradient, which corresponds to a stably stratified liquid. In this case, the anti-$\mathcal{PT}$ symmetry is no longer protected and can be broken spontaneously, which would lead to propagating, damped internal waves.

Anti-Hermiticity ensures that the eigenvalues $\omega = i\eta$ are purely imaginary for all $k_x, k_z, \mathrm{Ra}$, and $\mathrm{Pr}$. The transition to instability is given by the change of sign, $\eta < 0$ to $\eta > 0$. This necessarily occurs for $\omega = 0$, providing the marginal stability criterion for the Rayleigh-Bénard instability. This condition is satisfied
<figure><img src="image_5.png" /><figcaption>Fig. 3. Different regions of the Rayleigh-Bénard problem in the parameter space. In the upper part Ra > 0, anti-PT symmetry is protected, frequencies are purely imaginary, and the system exhibits an instability (black dashed) and no EPs. In the lower part Ra < 0, anti-PT is no longer protected by anti-Hermiticity, and the system exhibits a curve of EP2s (red) delimiting a region where perturbations propagate. It is a region where anti-PT symmetry is spontaneously broken.</figcaption></figure>

when

$$Ra = \frac{(k_x^2 + k_z^2)^3}{k_x^2}, \qquad (35)$$

whose minimum for each $k_x$, $k_z$ compatible with the boundary conditions gives the critical value of the Rayleigh number (e.g. $Ra_c = \frac{27}{4}\pi^4$ for stress-free, impenetrable, and isothermal boundaries). Figure 3 shows the regions of this problem.

## 6. Waves in dusty mixtures

Mixtures of pressureless dust and inviscid non-magnetised gas are generically used to model basic hydrodynamical properties of the dusty interstellar medium (e.g. Saffman 1962; Baines et al. 1965; Ahuja 1973; Gumerov et al. 1988). Exchange of momentum between gas and dust is modelled as a drag force proportional to the difference between their velocities. The DUSTYWAVE problem consist of studying the propagation of a 1D acoustic wave in a homogeneous dusty medium. The question is whether small disturbances propagate or are simply attenuated, as both regimes are possible (e.g. Laibe & Price 2011; David-Cléris & Laibe 2021). The mixture is initially homogeneous and at rest. Looking for small perturbative solutions as Fourier modes exp ($i(\omega t - kx)$), one obtains the eigenvalue equation

$$H_{dw}X = \omega X, \qquad (36)$$

$$H_{dw} = \begin{pmatrix} 0 & 0 & k & 0 \\ 0 & 0 & 0 & k \\ k & 0 & i\epsilon & -i(1-\epsilon) \\ 0 & 0 & -i\epsilon & i(1-\epsilon) \end{pmatrix}, \qquad (37)$$

after appropriate rescaling and choice of dimensionless parameters, as described in Appendix B. The operator $H_{dw}$ formally depends on two parameters: the wavenumber $k$ and the dust density fraction of the mixture $\epsilon$. Here, $H_{dw}$ satisfies Eq. (9) with $\tilde{U} = \text{diag}(1, 1, -1, -1)$, meaning that $H_{dw}$ is anti-PT symmetric. This property stems from the reflection symmetry $x \mapsto -x$ of

the fluid in physical space, as shown generically by David et al. (2022). Hence, the eigenvalues of $H_{dw}$ are either imaginary or consist of pairs of complex numbers with opposite real parts. In parameter space, the transition between propagating and non-propagative regions must then be an EP at which the solution of two counter-propagating waves degenerates.

The characterisation of these EPs is carried out according to the procedure described in Delplace et al. (2021) (see Appendix C for details). Exceptional points of order 2 (or EP2s) are located on a 1D manifold in the parameter space ($k, \epsilon$), which is a curve in a plane. Values of ($k, \epsilon$) for which an eigenvalue has a multiplicity of two are characterised by the fact that the characteristic polynomial $P(X)$ of $H_{dw}$ has a root of multiplicity of two. Such a root is therefore also a root of its derivative $P'(X)$, which means that the resultant $R_1 \equiv \mathcal{R}(P, P')$ between the two polynomials is 0. The resultant between two polynomials only cancels out if they have a common root. As $H_{dw}$ is anti-PT symmetric, $R_1$ is a real quantity. The values of ($k, \epsilon$) for which $R_1 = 0$ therefore define a curve. We calculate

$$R_1 = -4k^6(\epsilon - 1)^2(1 + k^4 - \epsilon + k^2(2 + 9\epsilon(-1 + \frac{3}{4}\epsilon))), \qquad (38)$$

and show the curve $R_1 = 0$ (top panel, red) in Fig. 4. This curve provides a simple alternative derivation of the result of David-Cléris & Laibe (2021). Non-propagating waves exist for $\epsilon > 8/9$, which is indeed the minimum of the curve of EPs.

Interestingly, in addition to the EP2s, the above method reveals the existence of third-order EPs for which three eigenvectors merge. These EP3s form a fold of dimension 0, which are two points in ($k, \epsilon$) space. These two points are located at

$$k = \pm \frac{1}{\sqrt{3}}, \quad \epsilon = \frac{8}{9}, \qquad (39)$$

which are the points at which the curve of the EP2s loses its regularity. At these EP3s, $H_{dw}$ has only two eigenvectors corresponding to the two eigenvalues 0 and $-i/3$. We note that these EP3s carry topological charges, also called winding numbers $W_3$:

$$W_3|_{EP3_1} = +1, \quad W_3|_{EP3_2} = -1, \qquad (40)$$

where $EP3_1$ is the point ($k, \epsilon$) = $(+\frac{1}{\sqrt{3}}, \frac{8}{9})$ and $EP3_2$ is the point ($k, \epsilon$) = $(-\frac{1}{\sqrt{3}}, \frac{8}{9})$ (see Appendix D). The question of connecting these topological charges to edge modes or particular global modes remains open (e.g. Leclerc et al. 2023).

According to the discussion in Sect. 2.2, $H_{dw}$ must be pseudo-chiral. The operator $\tilde{V}$ can be found explicitly

$$\tilde{V} = \begin{pmatrix} \epsilon & \epsilon - 1 & ik & 0 \\ \epsilon - 1 & \frac{(\epsilon - 1)^2}{\epsilon} & 0 & -\frac{ik(\epsilon - 1)}{\epsilon} \\ -ik & 0 & 0 & 0 \\ 0 & \frac{ik(\epsilon - 1)}{\epsilon} & 0 & 0 \end{pmatrix}. \qquad (41)$$

As has been shown, pseudo-chirality implies that if $\omega \neq -\omega^*$, the quantity $X \cdot \tilde{V}X$ must be zero. This condition is fulfilled when sound propagates. Although the system is not conservative, it still possesses an associated Krein quantity $X \cdot \tilde{V}X$:

$$X \cdot \tilde{V}X = \frac{1 - \epsilon}{\epsilon} |\rho_d|^2 + \frac{\epsilon}{1 - \epsilon} |\rho_g|^2 - 2 \left( \text{Re}(\rho_g \rho_d) - k \frac{\rho_0}{c_s} \text{Im}(\rho_g v_g^* + \rho_d v_d^*) \right), \qquad (42)$$
<figure><img src="image_6.png" /><figcaption>Fig. 4. Same as Fig. 1 but for the DUSTYWAVE problem. For $k = 0.4$, the Krein quantity $X \cdot \tilde{V}X$ is represented as a function of $\epsilon$. This Krein quantity provides a quadratic constant of motion for dusty sound waves in the propagating phase (low $\epsilon$).</figcaption></figure>

where $\rho_0 = \rho_{g,0} + \rho_{d,0}$ is the total background density. $X \cdot \tilde{V}X$ must be zero for propagating waves and non-zero for non-propagating waves. Figure 4 shows the numerical confirmation of this result. $X \cdot \tilde{V}X$ is only a constant of motion if it is zero, which is only the case for sound waves. Another technique, where the anti-$\mathcal{PT}$ operator is diagonalised, can be used to extract constants of motion in the unbroken phase, as shown for the Kelvin-Helmholtz instability by Qin et al. (2019) (see Appendix E for details).

Finally, it should be noted that in a monofluid description of the mixture, the variables are $\rho = \rho_g + \rho_d$, $v = \rho_g v_g + \rho_d v_d$, $\frac{\rho_d}{\rho_g}$ and $\Delta v = v_d - v_g$ (Laibe & Price 2014). If the perturbed quantities are denoted by $'$, the Krein quantity that cancels out for propagating perturbations is

$$X \cdot \tilde{V}_{\text{mono}} X = \frac{(1 - \epsilon)^3}{\epsilon} \left| \left( \frac{\rho_d}{\rho_g} \right)' \right|^2 + 2k \operatorname{Im} \left( (1 - \epsilon)^2 \left( \frac{\rho_d}{\rho_g} \right)' \Delta v'^* + \frac{\rho'}{\rho_0} v'^* \right).$$

This quantity is associated with the symmetry operator

$$\tilde{V}_{\text{mono}} = \begin{pmatrix} 0 & 0 & -ik & 0 \\ 0 & \frac{(\epsilon-1)^3}{\epsilon} & 0 & -ik(\epsilon-1)^2 \\ ik & 0 & 0 & 0 \\ 0 & ik(\epsilon-1)^2 & 0 & 0 \end{pmatrix}. \quad (43)$$

Squire & Hopkins (2018) show that when grains are additionally streaming through the gas, a generic resonant drag instability develops in the mixture. The system is neither $\mathcal{PT}$ nor anti-$\mathcal{PT}$ symmetric and Krein theory cannot be applied directly. The instability comes from a more complex mechanism where a resonance occurs between an anti-$\mathcal{PT}$ symmetric correction due to drag and a leading $\mathcal{PT}$ symmetric perturbation that arises from the streaming (e.g. Zhuravlev 2019; Magnan et al. 2024). Analysis of these problems from the point of view of discrete symmetries requires further specific study.

## 7. Gravitational instability of a dusty disc

To our knowledge, there are still no published theoretical predictions about what happens when systems with different discrete symmetries are combined. To illustrate this, we consider a final example, drawing on the analysis of Longarini et al. (2022), which deals with the gravitational instability of a dusty razor-thin disc. In the context of planet formation, the aim is to quantify the ways in which dust can favour the local collapse of gas or even clump itself into planetary embryos. This problem is a combination of the two problems presented in Sects. 3 and 6 and as such, is of order 6. With the notations used above, the matrix of linear perturbation in the local shearing box is

$$H_{\text{sgd}} = \begin{pmatrix} 0 & ik\Sigma_g & 0 & 0 & 0 & 0 \\ -ik \left( \frac{2\pi G}{|k|} - \frac{c_g^2}{\Sigma_g} \right) & \frac{\epsilon}{t_s} & -2\Omega & \frac{-2\pi Gk}{|k|} & -\frac{\epsilon}{t_s} & 0 \\ 0 & -2B & \frac{\epsilon}{t_s} & 0 & 0 & -\frac{\epsilon}{t_s} \\ 0 & 0 & 0 & 0 & ik\Sigma_d & 0 \\ -\frac{2\pi Gk}{|k|} & -\frac{1}{t_s} & 0 & -ik \left( \frac{2\pi G}{|k|} - \frac{c_g^2}{\Sigma_d} \right) & \frac{1}{t_s} & -2\Omega \\ 0 & 0 & -\frac{1}{t_s} & 0 & -2B & \frac{1}{t_s} \end{pmatrix}, \quad (44)$$

where $\Omega$ is the orbital frequency, $B$ is the local Oort parameter, and $c_g$ and $c_d$ denote the sound speeds of the gas and dust, respectively. Time and lengths are rescaled to the stopping time $t_s$ and the stopping length $c_g t_s$. The characteristic polynomial of the system is of order 5, meaning that the roots $\omega(k)$ cannot be determined analytically. In the regime of weak drag, Longarini et al. (2022) determines the marginal stability numerically.

On the other hand, this problem is $\mathcal{PT}$ symmetric, with the operator $U = \operatorname{diag}(-1, 1, 1, -1, 1, 1)$. The marginal stability curve is therefore a curve of EPs that corresponds to the spontaneous breaking of $\mathcal{PT}$ symmetry. From the analysis above, its exact expression is

$$\mathcal{R}_{\text{sgd}}(P, P') = 0, \quad (45)$$

a polynomial of high order with respect to the parameters of the problem. The explicit expression of Eq. (45) extends over two pages and is given in the worksheet in the Acknowledgements of the present paper. Hence, from a general perspective, the $\mathcal{PT}$-symmetric characterisation of a marginal stability criterion allows us to derive directly its analytical expression, albeit in a cumbersome manner, without approximations relying on asymptotically strong or weak drag regimes. Following Qin et al. (2021), this astrophysical example shows that the analysis of discrete symmetries therefore provides a powerful way to make analytical predictions about the marginal stability of symmetric systems with high dimensionalities (see Fig. 5).
$ε = 0.15, ξ = 0.005$

<figure><img src="image_7.png" /><figcaption>**Fig. 5.** Exceptional Points of the dusty self-gravitating disk are a consequence of the spontaneous breaking of $PT$ symmetry. The equation $R(P, P') = 0$ is obtained analytically with a resultant method, reproducing the results of Longarini et al. (2022).</figcaption></figure>

## 8. Conclusion

For several astrophysical objects, the properties of small perturbations around equilibrium are controlled by discrete symmetries, such as the $PT$ symmetry and the anti-$PT$ symmetry. For example, we show that such symmetries are relevant in the following contexts:

- The stability of isothermal spheres and thin self-gravitating discs: the Jeans and the Toomre problem are bisymmetric because of self-gravity.

- The propagation of waves and the onset of convection in stratified fluids for compressible and adiabatic perturbations: the Schwarzschild problem is bisymmetric because of buoyancy.

- The propagation of waves and the onset of convection in stratified fluids for incompressible and diffusive perturbations: the Rayleigh-Bénard problem is anti-$PT$ symmetric because of diffusion, and is also anti-Hermitian.

- The propagation of a wave in a dust-gas mixture: the DUSTYWAVE problem is anti-$PT$ symmetric because of drag.

In these systems, the parameter space can be divided into different regions, in which the discrete symmetries are either broken or not. These regions are separated by exceptional points (EPs) for which the system of linear perturbation can no longer be diagonalised.

From a methodology perspective, the analysis of discrete symmetries allows – when relevant – the direct determination of physical quantities that may be harder to obtain alternatively:

- The location of the EPs in the parameter space can be determined analytically from the symmetries without knowing the dispersion relation. This method is particularly effective for systems with high dimensions, where the dispersion relations may be untractable for analytical techniques, as in the case of the dusty self-gravitating disc.

- Krein theory then makes it possible to directly determine the exchange of stability when EPs are crossed, as well as the partition of energy enforced by Krein invariants associated with the discrete symmetries, even for dissipative systems.

Some of these EPs can be of high order and carry topological charges; this is the case for the DUSTYWAVE problem. We note

that not all physical systems allow such discrete symmetries. In particular, their spontaneous breaking is not a necessary condition for the system to be linearly unstable, as demonstrated by the Rayleigh-Bénard problem. We have focused here on the analysis of local stability, which applies to homogeneous systems. The study of global modes for inhomogeneous systems requires an extension of the Wigner-Weyl framework to non-Hermitian systems (e.g. Onuki 2020; Perez et al. 2021, 2022; Perez 2022; Leclerc et al. 2023).

**Acknowledgements.** We thank A. Venaille, P. Delplace, T. David-Cléris, E. Lynch, C. Longarini and S. Labrosse for stimulating discussions and comments. We thank the conference STAR@LYON for stimulation of this study. GL, AL and NP acknowledge the support of the European Research Council (ERC) CoG project PODCAST No. 864965. AL is funded by Contrat Doctoral Spécifique Normalien. We have used MATHEMATICA (Wolfram Research Inc. 2021), worksheets can be found at https://github.com/ArmandLeclerc/PTsymInAstroWaves.

## References

- Ahuja, A. S. 1973, J. Appl. Phys., 44, 4863

- Armitage, P. J. 2010, Astrophysics of Planet Formation (Cambridge University Press)

- Baines, M., Williams, I., Asebiomo, A., & Agacy, R. 1965, MNRAS, 130, 63

- Bénard, H. 1900, Revue Gen. Sci. Pure Appl., 11, 1261

- Bender, C. M., & Boettcher, S. 1998, Phys. Rev. Lett., 80, 5243

- Bender, C. M., Boettcher, S., & Meisinger, P. N. 1999, J. Math. Phys., 40, 2201

- Bender, C. M., Brody, D. C., & Jones, H. F. 2002, Phys. Rev. Lett., 89, 270401

- Bender, C. M., Brody, D. C., & Jones, H. F. 2004, Phys. Rev. D, 70, 025001

- Bergé, P., & Dubois, M. 1984, Contemp. Phys., 25, 535

- Bertin, G., & Lodato, G. 1999, A&A, 350, 694

- Bodenschatz, E., Pesch, W., & Ahlers, G. 2000, Ann. Rev. Fluid Mech., 32, 709

- Brandenburg, A. 2021, Lecture 4: Rayleigh-Bénard Problem http://norlx51.nordita.org/~brandenb/teach/AdvAstroFluids/4_Convection/notes.pdf

- Cowling, T. G. 1941, MNRAS, 101, 367

- David, T. W., Delplace, P., & Venaille, A. 2022, Phys. Fluids, 34, 056605

- David-Cléris, T., & Laibe, G. 2021, MNRAS, 504, 2889

- Delplace, P., Yoshida, T., & Hatsugai, Y. 2021, Phys. Rev. Lett., 127, 186602

- El-Ganainy, R., Makris, K. G., Christodoulides, D. N., & Musslimani, Z. H. 2007, Opt. Lett., 32, 2632

- El-Ganainy, R., Makris, K. G., Khajavikhan, M., et al. 2018, Nat. Phys., 14, 11

- Feng, L., Xu, Y.-L., Fegadolli, W. S., et al. 2013, Nat. Mater., 12, 108

- Feng, L., Wong, Z. J., Ma, R.-M., Wang, Y., & Zhang, X. 2014, Science, 346, 972

- Feng, L., El-Ganainy, R., & Ge, L. 2017, Nat. Photonics, 11, 752

- Fu, Y., & Qin, H. 2020, New J. Phys., 22, 083040

- Fu, Y., & Qin, H. 2024, Phys. Rev. Res., 6, 023273

- Ge, L., & Türeci, H. E. 2013, Phys. Rev. A, 88, 053810

- Goodman, J. 2003, MNRAS, 339, 937

- Goodman, J., & Narayan, R. 1988, MNRAS, 231, 97

- Gumerov, N., Ivandaev, A., & Nigmatulin, R. 1988, J. Fluid Mech., 193, 53

- Guo, A., Salamo, G. J., Duchesne, D., et al. 2009, Phys. Rev. Lett., 103, 093902

- Hodaei, H., Hassan, A. U., Wittek, S., et al. 2017, Nature, 548, 187

- Jeans, J. H. 1902, Philos. Trans. R. Soc. London Ser. A, 199, 1

- Kirillov, O. N. 2013, Nonconservative Stability Problems of Modern Physics (de Gruyter)

- Klaiman, S., Günther, U., & Moiseyev, N. 2008, Phys. Rev. Lett., 101, 080402

- Krein, M. 1950, Dokl. Akad. Nauk SSSR A, 73, 445

- Laibe, G., & Price, D. J. 2011, MNRAS, 418, 1491

- Laibe, G., & Price, D. J. 2014, MNRAS, 440, 2136

- Leclerc, A., Jezequel, L., Perez, N., et al. 2023, The Exceptional Ring of buoyancy instability in stars

- Lin, Z., Ramezani, H., Eichelkraut, T., et al. 2011, Phys. Rev. Lett., 106, 213901

- Longarini, C., Lodato, G., Bertin, G., & Armitage, P. J. 2022, MNRAS, 519, 2017

- Magnan, N., Heinemann, T., & Latter, H. N. 2024, MNRAS, 529, 688

- Makris, K. G., El-Ganainy, R., Christodoulides, D. N., & Musslimani, Z. H. 2008, Phys. Rev. Lett., 100, 103904

- Mostafazadeh, A. 2002a, J. Math. Phys., 43, 3944

- Mostafazadeh, A. 2002b, J. Math. Phys., 43, 205

- Mostafazadeh, A., & Batal, A. 2004, J. Phys. A: Math. General, 37, 11645
Musslimani, Z. H., Makris, K. G., El-Ganainy, R., & Christodoulides, D. N. 2008, Phys. Rev. Lett., 100, 030402

Onuki, Y. 2020, J. Fluid Mech., 883, A56

Özdemir, S. K., Rotter, S., Nori, F., & Yang, L. 2019, Nat. Materials, 18, 783

Peng, B., Özdemir, S. K., Lei, F., et al. 2014, Nat. Phys., 10, 394

Perez, N. 2022, Ph.D. Thesis, Ecole normale supérieure de lyon-ENS LYON

Perez, N., Delplace, P., & Venaille, A. 2021, Proc. R. Soc. A, 477, 20200844

Perez, N., Delplace, P., & Venaille, A. 2022, Phys. Rev. Lett., 128, 184501

Pringle, J. E., & King, A. 2007, Astrophysical Flows (Cambridge University Press)

Qin, H., Fu, Y., Glasser, A. S., & Yahalom, A. 2021, Phys. Rev. E, 104, 015215

Qin, H., Zhang, R., Glasser, A. S., & Xiao, J. 2019, Phys. Plasmas, 26, 032102

Rayleigh, L. 1916, Lond. Edinburgh Dublin Philos Mag. J. Sci., 32, 529

Regensburger, A., Bersch, C., Miri, M.-A., et al. 2012, Nature, 488, 167

Rüter, C. E., Makris, K. G., El-Ganainy, R., et al. 2010, Nat. Phys., 6, 192

Saffman, P. 1962, J. Fluid Mech., 13, 120

Schwarzschild, K. 1906, Nachrichten von der Königlichen Gesellschaft der Wissenschaften zu Göttingen Math.-phys. Klasse, 195, 41

Squire, J., & Hopkins, P. F. 2018, ApJ, 856, L15

Toomre, A. 1964, ApJ, 139, 1217

Wolfram Research Inc. 2021, Mathematica, Version 13.0.0, champaign, IL.

Zhang, F., Feng, Y., Chen, X., Ge, L., & Wan, W. 2020a, Phys. Rev. Lett., 124, 053901

Zhang, H., Huang, R., Zhang, S.-D., et al. 2020b, Nano Lett., 20, 7594

Zhang, R., Qin, H., & Xiao, J. 2020c, J. Math. Phys., 61, 012101

Zhuravlev, V. 2019, MNRAS, 489, 3850
# Appendix A: Rayleigh-Bénard: Linear perturbations

Here, we consider a two-dimensional fluid parametrised by the coordinates $(x, z)$. The fluid is at rest in a constant gravity field $g = -g\mathbf{u}_z$, and has constant density $\rho_0$. The temperature profile has a constant gradient $T(z) = T_h - \beta z$ with $\beta = (T_h - T_c)/L$. The coefficient of thermal expansion of the fluid is $\alpha = -\rho_0^{-1} (\partial\rho/\partial T)_p$. The dynamical viscosity and thermal diffusivity of the fluid are denoted by $\nu$ and $\kappa$, respectively. The linear system describing the evolution of small perturbations in the Boussinesq approximation around this background state is

$$ \partial_x v_x + \partial_z v_z = 0, \qquad (A.1) $$

$$ (\partial_t - \nu\Delta)v_x = -\frac{1}{\rho_0}\partial_x p, \qquad (A.2) $$

$$ (\partial_t - \nu\Delta)v_z = -\frac{1}{\rho_0}\partial_z p + \alpha g T, \qquad (A.3) $$

$$ (\partial_t - \kappa\Delta)T = \beta v_z. \qquad (A.4) $$

Applying double curl, and enforcing the incompressibility equation $\nabla \cdot \mathbf{v} = 0$, the equation for $u_x$ decouples from the other two equations, leaving

$$ (\partial_t - \nu\Delta)\Delta v_z = \alpha g \partial_{xx} T, \qquad (A.5) $$

$$ (\partial_t - \kappa\Delta)T = \beta v_z. \qquad (A.6) $$

For plane wave solutions of the form $\exp(-i\omega t + ik_z z + ik_x x)$, Eqs. A.5 – A.6 become

$$ (-i\omega - \nu(k_x^2 + k_z^2))(k_x^2 + k_z^2)v_z = -\alpha g k_x^2 T, \qquad (A.7) $$

$$ (-i\omega - \kappa(k_x^2 + k_z^2))T = \beta v_z. \qquad (A.8) $$

Rescaling the time by $L^2/\nu$, the lengths by $L$ and the temperature with $\beta L$ leads to the introduction of the Rayleigh and Prandtl numbers, $Ra = g\alpha\beta L^4/\nu\kappa$ and $Pr = \nu/\kappa$. The following rescaling of the equations $T \mapsto \frac{k_x}{k} \left(\frac{Ra}{Pr}\right)^{1/2} T$ with $k = \sqrt{k_x^2 + k_z^2}$ provides a system of eigenvalues that has the symmetric form

$$ \omega X = H_{RB}X, \qquad (A.9) $$

$$ H_{RB} = -i \begin{pmatrix} k^2 & -\frac{k_x}{k} \left(\frac{Ra}{Pr}\right)^{1/2} \\ -\frac{k_x}{k} \left(\frac{Ra}{Pr}\right)^{1/2} & k^2/Pr \end{pmatrix}, \qquad (A.10) $$

which corresponds to Eq. 32. Consequently, $iH_{RB}$ is Hermitian if and only if $Ra \ge 0$.

# Appendix B: DUSTYWAVE: linear perturbations

Here, we consider a dusty mixture where gas and dust are treated as two fluids interacting with each other via a drag force. For pressureless dust in 1D, the equations for mass and momentum conservation are

$$ (\partial_t + v_g \partial_x) \rho_g = -\rho_g \partial_x v_g, \qquad (B.1) $$

$$ (\partial_t + v_d \partial_x) \rho_d = -\rho_d \partial_x v_d, \qquad (B.2) $$

$$ (\partial_t + v_g \partial_x) v_g = \frac{K}{\rho_g} (v_d - v_g) - \frac{1}{\rho_g} \partial_x P, \qquad (B.3) $$

$$ (\partial_t + v_d \partial_x) v_d = \frac{K}{\rho_d} (v_g - v_d), \qquad (B.4) $$

where $\rho_g$, $\rho_d$ denote the gas and the dust densities, $v_g$, $v_d$ denote the gas and the dust velocities, and $K$ is the drag coefficient. The mixture is initially homogeneous and at rest. We

can write each quantity as a background and perturbation component $f = f_0 + f'$. The system of equations can be put in a symmetric form using the rescaled quantities $(\rho'_g, \rho'_d, v'_g, v'_d) \mapsto X = \left(\frac{1}{\rho_{g,0}+\rho_{d,0}}\rho'_g, \frac{1}{\rho_{g,0}+\rho_{d,0}}\rho'_d, \frac{\rho_{g,0}}{c_s(\rho_{g,0}+\rho_{d,0})}v'_g, \frac{\rho_{d,0}}{c_s(\rho_{g,0}+\rho_{d,0})}v'_d\right)^T$. For isothermal perturbations $P'(x, t) = c_s^2 \rho'_g(x, t)$, the evolution of the perturbations are given by the system

$$ \partial_t X = \begin{pmatrix} 0 & 0 & -c_s \partial_x & 0 \\ 0 & 0 & 0 & -c_s \partial_x \\ -c_s \partial_x & 0 & \frac{K}{\epsilon K} + \frac{(1-\epsilon)K}{\epsilon \rho_{g,0}} & 0 \\ 0 & 0 & \frac{1-\epsilon}{(1-\epsilon)\rho_{d,0}} & -\frac{\rho_{g,0}}{\rho_{d,0}} \end{pmatrix} X, \qquad (B.5) $$

where $\epsilon \equiv \frac{\rho_{d,0}}{\rho_{g,0}+\rho_{d,0}}$ is the dust fraction, and $1 - \epsilon = \frac{\rho_{g,0}}{\rho_{g,0}+\rho_{d,0}}$ is the gas fraction. The barycentric stopping time $t_s \equiv \frac{\rho_{g,0}}{K(\rho_{g,0}+\rho_{d,0})}$ corresponds to the time needed for drag to damp the differential velocity between gas and dust. The stopping length $c_s t_s$ is the typical spatial separation between the fluids. Using $t_s$ and $c_s t_s$ as units of time and length, we obtain the dimensionless system

$$ H_{dw} X = i \partial_t X, \qquad (B.6) $$

$$ H_{dw} \equiv \begin{pmatrix} 0 & 0 & k & 0 \\ 0 & 0 & 0 & k \\ k & 0 & i\epsilon & -i(1-\epsilon) \\ 0 & 0 & -i\epsilon & i(1-\epsilon) \end{pmatrix}. \qquad (B.7) $$

# Appendix C: EPs from resultants

Here, we consider two polynomials $Q1(X)$ and $Q2(X)$ defined by

$$ Q_1(X) = \sum_{i=1}^{n} a_i X^i = a_n \prod_{i=1}^{n} (X - \alpha_i), \qquad (C.1) $$

$$ Q_2(X) = \sum_{j=1}^{m} b_j X^j = b_m \prod_{j=1}^{m} (X - \beta_j), \qquad (C.2) $$

where $\alpha_i$ and $\beta_j$ are their roots. The resultant $\mathcal{R}(Q_1, Q_2)$ of these two polynomials is defined as

$$ \mathcal{R}(Q_1, Q_2) \equiv \det S(Q_1, Q_2), \qquad (C.3) $$

where $S(Q_1, Q_2)$ is the Sylvester matrix of $Q_1$ and $Q_2$. It satisfies

$$ \mathcal{R}(Q_1, Q_2) = a_n^m b_m^n \prod_{i,j} (\alpha_i - \beta_j), \qquad (C.4) $$

which cancels if and only if $Q_1$ and $Q_2$ share a common root. The quantity (C.3) is determined by the coefficients $a_i, b_j$ and does not require that we find the roots $\alpha_i, \beta_j$. We can define $P'$, the derivative of $P(X)$. The cancellation of the resultant $R_1$ of $P$ and $P'$ characterise the fact that $P$ has a root with double multiplicity. In the context of this study, an EP of order at least 2 is the space formed by parameters for which

$$ R_1 \equiv \mathcal{R}(P, P') = 0. \qquad (C.5) $$

# Appendix D: DUSTYWAVE: topological charges

In addition to EPs of order 2, the DUSTYWAVE is also found to have EPs of order 3. The structure of an EP3 is a fold of dimension 0, which is a point in the $(k, \epsilon)$ space. These points can be located using a procedure that generalises the determination of
For the other two points, $P$ and $P''$ share a root different from the one shared by $P$ and $P'$. As such, they are EP2s only. At the EP3s, $H$ has two eigenvalues, 0 and $-\frac{i}{3}$. It cannot be diagonalised, as it only has two eigenvectors,

$$ \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} -i\sqrt{3}/2 \\ 2i\sqrt{3} \\ -1/2 \\ 2 \end{pmatrix}. \qquad (D.3) $$

The EP3s have a non-zero topological charge, as defined in (Delplace et al. 2021). These topological invariants are defined as the winding numbers

$$ W_3 = \frac{1}{2\pi} \oint_C \frac{1}{R_1^2 + R_2^2} (R_1 \nabla_p R_2 - R_2 \nabla_p R_1) d\mathbf{p}, \qquad (D.4) $$

for $\nabla_p = (\partial_k)$ and $\mathbf{p} = (k/\epsilon) \in C$, where $C$ is any closed loop around one of the EP3s. We compute

$$ W_3|_{EP3_1} = +1, \qquad (D.5) $$

$$ W_3|_{EP3_2} = -1, \qquad (D.6) $$

where $EP3_1$ is the point at $(k, \epsilon) = (+\frac{1}{\sqrt{3}}, \frac{8}{9})$ and $EP3_2$ is the point at $(k, \epsilon) = (-\frac{1}{\sqrt{3}}, \frac{8}{9})$.

## Appendix E: DUSTYWAVE: Anti-$\mathcal{PT}$ operator

We define $\Theta$ as the operator of complex conjugation and $CP \equiv \tilde{U}\Theta$ as the operator for anti-$\mathcal{PT}$ symmetry. Equation (9) can be rewritten in the form

$$ H_{dw}CP + CPH_{dw} = 0. \qquad (E.1) $$

Equation E.1 constrains the phase between the components of the eigenvectors of $H_{dw}$, which are also called relations of polarisation. Here, we consider $X_0$ an eigenvector of $H_{dw}$ with pure imaginary eigenvalue. Formally,

$$ H_{dw}X_0 = i\eta_0 X_0, \qquad (E.2) $$

where $\eta_0$ is real. Applying the operator $CP$ on Eq. E.2 and the commutation rule Eq.(E.1) gives the relation

$$ H_{dw}(CPX_0) = i\eta_0(CPX_0). \qquad (E.3) $$

If the eigenvalue is not degenerated (i.e. the parameters are chosen for the system not to be at an EP), $CPX_0$ is colinear to $X_0$, i.e.

$$ CPX_0 = \alpha X_0, \qquad (E.4) $$

for some complex number $\alpha$. Hence, $H_{dw}$ and $CP$ share the same eigenvectors and the anti-$\mathcal{PT}$ symmetry is said to be unbroken. We note that the same reasoning applied to the operator $iH_{dw}$ shows that $iH_{dw}$ is $\mathcal{PT}$ symmetric. $CP$ can be straightforwardly diagonalised. Its eigenvectors are

$$ \begin{pmatrix} 1 \\ 0 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 1 \\ 0 \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ i \\ 0 \end{pmatrix}, \begin{pmatrix} 0 \\ 0 \\ 0 \\ i \end{pmatrix}. \qquad (E.5) $$

<figure><img src="image_8.png" /><figcaption>Fig. D.1. Exceptional points of the DUSTYWAVE in the parameter space $(k, \epsilon)$. Intersections between the red ($R_1 = 0$) and blue ($R_2 = 0$) curves indicate the positions of EP3s. The two bottom panels show zoomed-in images of two regions. The EP3s are marked by black points, with associated topological charge of $\pm 1$.</figcaption></figure>

EP2s. Noting $R_1 = \mathcal{R}(P, P')$ and $R_2 = \mathcal{R}(P, P'')$, the resultants of the characteristic polynomial and its derivatives, and EP3s correspond to points in the parameter space for which $P, P'$ and $P''$ have a common root. These points are those where both $R_1$ and $R_2$ cancel simultaneously. For the DUSTYWAVE problem, $R_1$ is given by Eq. 38, and $R_2$ by

$$ R_2 = 400k^8 - 432k^4(\epsilon - 1) + 96k^6(14 - 51\epsilon + 36\epsilon^2). \qquad (D.1) $$

Fig.D.1 shows the curves $R_1 = 0$ (red) and $R_2 = 0$ (blue). These curves cross four times, and two of these points are EP3s. They are located at

$$ k = \pm \frac{1}{\sqrt{3}}, \quad \epsilon = \frac{8}{9}. \qquad (D.2) $$
The condition given by Eq. E.4 therefore imposes the particular
polarisation relations on non-propagating waves:

$$
\left( \frac{\rho_g'}{\rho_d'} \right)^* = + \frac{\rho_g'}{\rho_d'}, \qquad (E.6)
$$

$$
\left( \frac{\rho_g'}{v_g'} \right)^* = - \frac{\rho_g'}{v_g'}, \qquad (E.7)
$$

$$
\left( \frac{\rho_d'}{v_d'} \right)^* = - \frac{\rho_d'}{v_d'}, \qquad (E.8)
$$

and other linear combinations of Eqs. E.6 – E.8. Therefore, the
analysis of the discrete symmetries of the problem shows with a
very small number of calculations that gas and dust densities are
in phase, as are their velocities, while densities and velocities are
in phase quadrature. Similar relations exist between fields in the
unbroken *PT* symmetry phase of the Kelvin-Helmholtz problem
discussed by Fu & Qin (2020).

This result can be reformulated with real fields. In the unbro-
ken anti-PT symmetry phase, the perturbations are damped
without oscillating, and are of the form

$$
\rho_g' \propto e^{-\eta_0 t} \cos(kx), \qquad (E.9)
$$

$$
\rho'_{d} \propto e^{-\eta_{0} t} \cos(kx), \qquad (E.10)
$$

$$
v'_{g} \propto e^{-\eta_{0} t} \sin(kx), \qquad (E.11)
$$

$$
v'_{d} \propto e^{-\eta_{0} t} \sin(kx), \qquad (E.12)
$$

which imposes that quantities be zero when they are spatially
averaged. Formally, with

$$
\langle f \rangle(t) \equiv \lim_{L \to \infty} \frac{1}{L} \int_{-L/2}^{L/2} f(x, t) dx, \qquad (E.13)
$$

one obtains for example

$$
\langle v'_{g} \rho'_{g} \rangle = 0, \qquad (E.14)
$$

$$
\langle v'_{d} \rho'_{d} \rangle = 0, \qquad (E.15)
$$

for all times. On the contrary, in the phase where anti-PT sym-
metry is broken, the frequency of the perturbations has a real
part. Waves are propagating, and the polarisation relations are
different from those in Eqs.(E.6)-(E.8), and the time-averages
of quantities defined in Eqs. E.14 – E.15 are non-zero for
this phase. Only propagating waves have non-zero macroscopic
mass flux.