© The Authors 2024

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 𝒫 and time 𝒯.

Some of the most striking properties of 𝒫𝒯-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 𝒫𝒯-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). 𝒫𝒯 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 𝒫𝒯 symmetry (El-Ganainy et al. 2007; Makris et al. 2008). Klaiman et al. (2008) identified a spontaneous breaking of 𝒫𝒯 symmetry in a waveguide. Musslimani et al. (2008) investigated the non-linear propagation of optical solitons in 𝒫𝒯-symmetric refraction indexes. (Guo et al. 2009) and (Rüter et al. 2010) made the first observations of 𝒫𝒯 symmetry in optics with a complex index. A decade of developments of 𝒫𝒯 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 𝒫𝒯 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 𝒫𝒯 symmetry induces non-reciprocal propagation. Lin et al. (2011), Regensburger et al. (2012) and Feng et al. (2013) used 𝒫𝒯 symmetry in periodic systems (through gratings, periodic crystals, and lattices) to design unidirectional invisibility. 𝒫𝒯-symmetric degeneracies are now used for their high resonant sensitivity in optical cavities (Hodaei et al. 2017). Anti-𝒫𝒯 symmetry has also found application in novel optical designs (e.g. Ge & Türeci 2013; Zhang et al. 2020a,b). The topological properties of 𝒫𝒯-symmetric systems have recently been investigated in the context of fluid dynamics (Fu & Qin 2024). 𝒫𝒯 and anti-𝒫𝒯 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 𝒫𝒯 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 𝒫𝒯-symmetric systems exhibit particular properties.

In this study, we aim to show that 𝒫𝒯 and anti-𝒫𝒯 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 𝒫𝒯 and anti-𝒫𝒯 symmetries in Sect. 2. We then demonstrate how they condition the outcome of the linear 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. 𝒫𝒯 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 𝒫𝒯 symmetric if there exists a unitary operator U such that

$$\begin{array}{c}\hfill U^H{U}^{-1}=H^{\ast },\end{array}$$(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⁻¹X^*) = ω^*(U⁻¹X^*). In 𝒫𝒯-symmetric systems, unstable modes are therefore accompanied by damped modes (imaginary parts of opposite sign). When ω ≠ ω^*, 𝒫𝒯 symmetry is said to be spontaneously broken by H. By defining the symmetry operator 𝒫𝒯 ≡ UΘ, where Θ denotes complex conjugation, one has [𝒫𝒯, 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 𝒫𝒯 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

$$\begin{array}{c}\hfill VH{V}^{-1}=H^{†},\end{array}$$(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\equiv Y_1^{*\top }{Y}_2$.

Let X₁ and X₂ be two eigenvectors of H satisfying Eq. (2). The eigenvalues ω₁ and ω₂ 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, ω₁ and ω₂ are real, 𝒫𝒯 symmetry is said to be unbroken, and the Krein quantities defined as X_1, 2 ⋅ VX_1, 2 are non-zero and of opposite sign. After the collision, the system is unstable, the complex frequencies ${\omega }_1={\omega }_2^{*}$ are complex conjugates of each other, and the two Krein quantities are necessarily zero. The Krein quantity Y ⋅ VY is a constant of motion of a general solution Y(t, a) of the evolution equation i∂_tY = H(a)Y. Indeed,

$$\begin{array}{cc}\hfill i{\partial }_t(Y·VY)& =Y·(-H^{†}V)Y+Y·VHY\hfill \end{array}$$(3)

$$\begin{array}{cc}\hfill & =Y·(VH-H^{†}V)Y\hfill \end{array}$$(4)

$$\begin{array}{cc}\hfill & =0.\hfill \end{array}$$(5)

Hence,

$$\begin{array}{c}\hfill Y·VY=\mathrm{cst}.\end{array}$$(6)

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

$$\begin{array}{cc}\hfill X·VHX& =X·H^{†}VX={\omega }^{\ast }X·VX\hfill \end{array}$$(7)

$$\begin{array}{cc}& =\omega X·VX.\hfill \end{array}$$(8)

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

2.2. Anti-𝒫𝒯 symmetry

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

$$\begin{array}{c}\hfill {\stackrel{\sim }{U}}^H{\stackrel{\sim }{U}}^{-1}=-H^{\ast }.\end{array}$$(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-𝒫𝒯 symmetric is equivalent to iH being 𝒫𝒯 symmetric. The same equivalence between 𝒫𝒯 symmetry and pseudo-Hermiticity then applies, and guarantees that an anti-𝒫𝒯-symmetric operator H is also pseudo-chiral: there exist an Hermitian operator $\stackrel{\sim }{V,}$ such that

$$\begin{array}{c}\hfill \stackrel{\sim }{V}H{\stackrel{\sim }{V}}^{-1}=-H^{†}.\end{array}$$(10)

The quantity $X·\stackrel{\sim }{V}X$ is not strictly a Krein quantity, and is not in general a constant of motion of a solution of i∂_tX = HX. However, for a Fourier mode, Y = e^−iωtX, one has

$$\begin{array}{c}\hfill Y·\stackrel{\sim }{V}Y={\mathrm{e}}^{2\mathrm{Im}(\omega )t}X·\stackrel{\sim }{V}X,\end{array}$$(11)

and is constant to zero for modes with Re(ω)≠0; that is, for modes that spontaneously break the anti-𝒫𝒯 symmetry. Indeed,

$$\begin{array}{cc}\hfill X·\stackrel{\sim }{V}HX& =X·(-H^{†}\stackrel{\sim }{V})X=-{\omega }^{\ast }X·\stackrel{\sim }{V}X\hfill \end{array}$$(12)

$$\begin{array}{cc}& =\omega X·\stackrel{\sim }{V}X.\hfill \end{array}$$(13)

In summary, whenever ω ≠ −ω^*, $X·\stackrel{\sim }{V}X=0$, the Krein quantity is zero at all times. An analogue of a Krein collision may occur. When the anti-𝒫𝒯 symmetry is unbroken, their exists a pair of eigenvalues of H, denoted ω₁ and ω₂, 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 ${\omega }_1=-{\omega }_2^{*}$ with non-zero real parts. Then, the anti-𝒫𝒯 symmetry is spontaneously broken. This phenomenon corresponds to a Krein collision for the operator iH. It follows that the quantity $X·\stackrel{\sim }{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 𝒫𝒯 and anti-𝒫𝒯 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-𝒫𝒯 symmetry. Conversely, when ω is imaginary, the associated eigenvector spontaneously breaks 𝒫𝒯 symmetry. Thus, one of those two symmetries is necessarily broken for any value of the parameters a, and either X ⋅ VX = 0 (broken 𝒫𝒯) or $X·\stackrel{\sim }{V}X=0$ (broken anti-𝒫𝒯). In bisymmetric systems, X ⋅ VX is always a constant of motion, while $X·\stackrel{\sim }{V}X$ is constant only when it is zero, that is, in the spontaneously broken anti-𝒫𝒯 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 ρ₀ 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:

$$\begin{array}{c}\hfill \omega \left(\begin{array}{c}{\rho }^{\prime }/{\rho }_0\\ {\mathit{v}}^{\prime }/c_{\mathrm{s}}\end{array}\right)=H_{\mathrm{J}}\left(\begin{array}{c}{\rho }^{\prime }/{\rho }_0\\ {\mathit{v}}^{\prime }/c_{\mathrm{s}}\end{array}\right)\equiv \left(\begin{array}{cc}0& c_{\mathrm{s}}{\mathit{k}}^{\top }\\ (1-\frac{4\pi \mathcal{G}{\rho }_0}{c_{\mathrm{s}}^2{|k|}^2})c_{\mathrm{s}}\mathit{k}& 0\end{array}\right)\left(\begin{array}{c}{\rho }^{\prime }/{\rho }_0\\ {\mathit{v}}^{\prime }/c_{\mathrm{s}}\end{array}\right).\end{array}$$(14)

This problem is invariant under both parity and time reversal. It is therefore trivially 𝒫𝒯 symmetric, as proven by the statement $H_{\text{J}}=H_{\text{J}}^{*}$. As such, H_J is pseudo-Hermitian, with the symmetry operator $V=\mathrm{diag}((1-\frac{4\pi \mathcal{G}{\rho }_0}{c_{\mathrm{s}}^2{|k|}^2}),1,1,1).$ The Krein quantity associated to the Krein collision of the Jeans instability is then $X·VX=(1-\frac{4\pi \mathcal{G}{\rho }_0}{c_{\mathrm{s}}^2{|k|}^2}){|{\rho }^{\prime }|}^2+{|{\mathit{v}}^{\prime }|}^2$. This quantity is non-zero for propagating waves, and zero for the unstable modes. Complementarily, H_J is anti-𝒫𝒯 symmetric and pseudo-chiral; indeed, it satisfies Eq. (10) for $\stackrel{\sim }{V}=\mathrm{diag}(-(1-\frac{4\pi \mathcal{G}{\rho }_0}{c_{\mathrm{s}}^2{|k|}^2}),1,1,1)$. The quantity $X·\stackrel{\sim }{V}X=-(1-\frac{4\pi \mathcal{G}{\rho }_0}{c_{\mathrm{s}}^2{|k|}^2}){|{\rho }^{\prime }|}^2+{|{\mathit{v}}^{\prime }|}^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 $|\mathit{k}|=k_{\mathrm{J}}\equiv \sqrt{4\pi \mathcal{G}{\rho }_0}/c_{\mathrm{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, 𝒫 and 𝒯 are broken individually, but the overall system of linear perturbations is still 𝒫𝒯 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_{\mathrm{s}}\kappa }{\pi \mathcal{G}{\mathrm{\Sigma }}_0}>1$, where κ denotes the epicyclic frequency of the disc and Σ₀ its surface density. The razor-thin disc is integrated vertically (Armitage 2010), to give the following set of radial linear perturbations:

$$\begin{array}{c}\hfill H_{\mathrm{sgd}}\left(\begin{array}{c}{h}^{\prime }\\ v_r^{\prime }\\ v_{\varphi }^{\prime }\end{array}\right)=\left(\begin{array}{ccc}0& c_{\mathrm{s}}k& 0\\ -\frac{k}{c_{\mathrm{s}}}(\frac{2}{Q|k|}-1)& 0& 2i\\ 0& -\frac{i}{2}& 0\end{array}\right)\left(\begin{array}{c}h\prime \\ v_r^{\prime }\\ v_{\varphi }^{\prime }\end{array}\right)=\omega \left(\begin{array}{c}{h}^{\prime }\\ v_r^{\prime }\\ v_{\varphi }^{\prime }\end{array}\right),\end{array}$$(15)

where we use perturbations of the enthalpy h′ and the two horizontal velocities $v_r^{\prime }$ and $v_{\varphi }^{\prime }$ as variables, and the orbital time Ω⁻¹ and the pressure length c_sΩ⁻¹ as units of time and length.

With U = diag(1, 1, −1), this problem satisfies Eq. (1) and is therefore 𝒫𝒯 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

$$\begin{array}{c}\hfill V=\left(\begin{array}{ccc}\frac{1}{c_{\mathrm{s}}^2}& 0& -\frac{4ik}{c_{\mathrm{s}}Q|k|}\\ 0& 1& 0\\ \frac{4ik}{c_{\mathrm{s}}Q|k|}& 0& 4(\frac{2|k|}{Q}+1)\end{array}\right),\end{array}$$(16)

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

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill X·VX=0=\frac{{|h^{\prime }|}^2}{c_{\mathrm{s}}^2}+{|v_r^{\prime }|}^2& +4|v_{\varphi }{\prime |}^2(1+\frac{2\pi \mathcal{G}{\mathrm{\Sigma }}_0}{{\mathrm{\Omega }}^2}|k|)\hfill \\ \hfill & -\frac{8\pi \mathcal{G}{\mathrm{\Sigma }}_0}{c_{\mathrm{s}}^2\mathrm{\Omega }}\text{sgn}(k)\mathrm{Im}\left(h{v}_{\varphi }^{\ast }\right),\hfill \end{array}\end{array}$$(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_sgd has degenerated eigenvalues when $Q=\frac{2|k|}{1+k^2}$. These are the positions of EPs for H_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.

Fig. 1. 𝒫𝒯-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 ⋅ VX, orange is $X·\stackrel{\sim }{V}X$, solid blue is Re(ω), and dotted blue is Im(ω), for a destabilising mode.

In addition to being 𝒫𝒯 symmetric, H_sgd is also anti-𝒫𝒯 symmetric, because

$$\begin{array}{c}\hfill \stackrel{\sim }{U}{H}_{\mathrm{sgd}}{\stackrel{\sim }{U}}^{-1}=-H_{\mathrm{sgd}}^{\ast },\end{array}$$(18)

with $\stackrel{\sim }{U}=\mathrm{diag}(1,-1,-1)$. The anti-𝒫𝒯 symmetry is complementary to the 𝒫𝒯 symmetry, and the combination of the two symmetries implies that the eigenvalues always come in pairs: they are always opposite (ω, −ω), 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

$$\begin{array}{cc}\hfill \stackrel{\sim }{V}{H}_{\mathrm{sgd}}{\stackrel{\sim }{V}}^{-1}& =-H_{\mathrm{sgd}}^{†},\hfill \end{array}$$(19)

$$\begin{array}{ccc}\hfill \stackrel{\sim }{V}& \equiv \hfill & \hfill \left(\begin{array}{ccc}\frac{1}{c^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{array}\right).\end{array}$$(20)

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

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill X·\stackrel{\sim }{V}X=& \frac{{|h^{\prime }|}^2}{c_{\mathrm{s}}^2}+|v_r{\prime |}^2+{|v_{\varphi }^{\prime }|}^2(-4+8(\frac{c_{\mathrm{s}}^2{k}^2}{{\mathrm{\Omega }}^2}-\frac{\pi \mathcal{G}{\mathrm{\Sigma }}_0}{{\mathrm{\Omega }}^2}|k|))\hfill \\ \hfill & -\frac{8\pi \mathcal{G}{\mathrm{\Sigma }}_0}{c_{\mathrm{s}}^2\mathrm{\Omega }}\text{sgn}(k)(\frac{c_{\mathrm{s}}^2}{\pi \mathcal{G}{\mathrm{\Sigma }}_0}|k|-1)\mathrm{Im}\left(h{v}_{\varphi }^{\ast }\right).\hfill \end{array}\end{array}$$(21)

Hence, for each mode, either X ⋅ VX or $X·\stackrel{\sim }{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 ω = 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 $\stackrel{\sim }{V}$ identified above are particular cases of more general families of matrices:

$$\begin{array}{c}\hfill V=\left(\begin{array}{ccc}\frac{2c_{\mathrm{s}}(-2k^2+(-1+k^2)Q\left|k\right|)a+k(-2+Q\left|k\right|)c}{4c_{\mathrm{s}}^{2}kQ\left|k\right|}& b\frac{(kQ-2\mathrm{sgn}k)}{2c_{\mathrm{s}}Q}& \mathit{ia}\\ b\frac{(kQ-2\mathrm{sgn}k)}{2c_{\mathrm{s}}Q}& \frac{1}{4}(2c_{\mathrm{s}}ka+c)& \mathit{ib}\\ -ia& -ib& c\end{array}\right),\end{array}$$(22)

where a, b, and c are any real numbers such that detW ≠ 0. The matrix of Eq. (16) is recovered with $a=-\frac{4k}{c_{\mathrm{s}}Q\left|k\right|}$, b = 0, and $c=\frac{8\left|k\right|}{Q}+2$. Similarly,

$$\begin{array}{c}\hfill \stackrel{\sim }{V}=\left(\begin{array}{ccc}\frac{2c_{\mathrm{s}}(-2k^2+(-1+k^2)Q\left|k\right|)a+k(-2+Q\left|k\right|)c}{4c_{\mathrm{s}}^{2}kQ\left|k\right|}& b\frac{i(kQ-2\mathrm{sgn}k)}{2c_{\mathrm{s}}Q}& \mathit{ia}\\ -b\frac{i(kQ-2\mathrm{sgn}k)}{2c_{\mathrm{s}}Q}& -\frac{1}{4}(2c_{\mathrm{s}}ka+c)& b\\ -ia& b& c\end{array}\right).\end{array}$$(23)

The matrix of Eq. (20) is recovered with $a=-\frac{4k(Q|k|-1)}{c_{\mathrm{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² 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,

$$\begin{array}{cc}\hfill \omega X& =H_{\mathrm{b}}X,\hfill \end{array}$$(24)

$$\begin{array}{cc}\hfill H_{\mathrm{b}}& =\left(\begin{array}{cccc}0& 0& 0& L_{\ell }\\ 0& 0& |N|& K_r-iS\\ 0& -|N|& 0& 0\\ L_{\ell }& K_r+iS& 0& 0\end{array}\right),\hfill \end{array}$$(25)

where ω is the complex eigenfrequency of the perturbation, and X ≡ (v w Θ p)^⊤ contains the perturbation’s horizontal velocity, radial velocity, entropy, and pressure after appropriate rescaling. $N^2\equiv -g\frac{\mathrm{d}ln{\rho }_0}{\mathrm{d}r}-\frac{g^2}{c_{\mathrm{s}}^2}$ is the square of the buoyancy frequency, which is negative here, $S\equiv \frac{c_{\mathrm{s}}}{2g}(N^2-\frac{g^2}{c_{\mathrm{s}}^2})-\frac{1}{2}\frac{\mathrm{d}{c}_{\mathrm{s}}}{\mathrm{d}r}+\frac{c_{\mathrm{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_{\mathrm{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:

$$\begin{array}{cc}\hfill V& =\mathrm{diag}(1,1,-1,1),\hfill \end{array}$$(26)

$$\begin{array}{cc}\hfill \stackrel{\sim }{V}& =\mathrm{diag}(1,1,1,-1).\hfill \end{array}$$(27)

Krein theory can therefore be applied to determine how discrete symmetries constrain the partition of energy. Here, X ⋅ 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 ⋅ VX is a constant of motion for any solution of ∂_tX = H_bX, as shown in Sect. 2.1. Moreover, this constant is exactly zero only for unstable modes. In contrast, $X·\stackrel{\sim }{V}X$ also grows exponentially for unstable modes, but is constantly equal to zero for stable modes. As such,

$$\begin{array}{cc}& X·X={|v|}^2+{|w|}^2+{|\mathrm{\Theta }|}^2+{|p|}^2\propto {\mathrm{e}}^{2\mathrm{Im}(\omega )t},\hfill \end{array}$$(28)

$$\begin{array}{cc}& X·VX={|v|}^2+{|w|}^2-{|\mathrm{\Theta }|}^2+{|p|}^2=\mathrm{cst}{\delta }_{\mathrm{Im}(\omega ),0},\hfill \end{array}$$(29)

$$\begin{array}{cc}& X·\stackrel{\sim }{V}X={|v|}^2+{|w|}^2+{|\mathrm{\Theta }|}^2-{|p|}^2\propto {\mathrm{e}}^{2\mathrm{Im}(\omega )t}{\delta }_{\mathrm{Re}(\omega ),0},\hfill \end{array}$$(30)

where δ is the Kronecker delta.

For the radial modes 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 ⋅ VX is zero and $X·\stackrel{\sim }{V}X$ is non-zero. Outside the ring, the modes are stable and propagate as sound waves, and X ⋅ VX is non-zero and $X·\stackrel{\sim }{V}X$ is zero.

Fig. 2. Same as Fig. 1 but for the Schwarzschild instability. Non-hermiticity is broken in this case by buoyancy instead of self-gravity.

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 $\stackrel{\sim }{V}$ have ±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αβL⁴/νκ, which needs to be greater than some critical value of order 10² − 10³ 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 ∂_zT₀ = −β is directed. With dimensionless variables, the problem for Fourier modes exp( − iωt + ik_zz + ik_xx) can be converted into the symmetrical form:

$$\begin{array}{cc}& \omega \left(\begin{array}{c}\stackrel{\sim }{v_z}\\ [0.2cm]\stackrel{\sim }{T}\end{array}\right)=H_{\mathrm{RB}}\left(\begin{array}{c}\stackrel{\sim }{v_z}\\ [0.2cm]\stackrel{\sim }{T}\end{array}\right),\hfill \end{array}$$(31)

$$\begin{array}{cc}& H_{\mathrm{RB}}=-i\left(\begin{array}{cc}{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{array}\right),\hfill \end{array}$$(32)

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

This problem is anti-𝒫𝒯 symmetric and pseudo-chiral, as one reads that

$$\begin{array}{c}\hfill H_{\mathrm{RB}}=-H_{\mathrm{RB}}^{\ast }=-H_{\mathrm{RB}}^{†},\end{array}$$(33)

for Ra ≥ 0. In this case, the anti-𝒫𝒯 symmetry is protected and cannot be spontaneously broken as H_RB is anti-Hermitian (i.e. iH_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

$$\begin{array}{c}\hfill \mathrm{Ra}=(1-\frac{{(1+\mathrm{Pr})}^2}{4\mathrm{Pr}})\frac{{(k_x^2+k_z^2)}^3}{k_x^2},\end{array}$$(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-𝒫𝒯 symmetry is no longer protected and can be broken spontaneously, which would lead to propagating, damped internal waves.

Anti-Hermiticity ensures that the eigenvalues ω = iη are purely imaginary for all k_x, k_z, Ra, and Pr. The transition to instability is given by the change of sign, η < 0 to η > 0. This necessarily occurs for ω = 0, providing the marginal stability criterion for the Rayleigh-Bénard instability. This condition is satisfied when

$$\begin{array}{c}\hfill \mathrm{Ra}=\frac{{(k_x^2+k_z^2)}^3}{k_x^2},\end{array}$$(35)

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

Fig. 3. Different regions of the Rayleigh-Bénard problem in the parameter space. In the upper part Ra > 0, anti-𝒫𝒯 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-𝒫𝒯 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-𝒫𝒯 symmetry is spontaneously broken.

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(ωt − kx)), one obtains the eigenvalue equation

$$\begin{array}{cc}& H_{\mathrm{dw}}X=\omega X,\hfill \end{array}$$(36)

$$\begin{array}{cc}& H_{\mathrm{dw}}\equiv \left(\begin{array}{cccc}0& 0& k& 0\\ 0& 0& 0& k\\ k& 0& iϵ& -i(1-ϵ)\\ 0& 0& -iϵ& i(1-ϵ)\end{array}\right),\hfill \end{array}$$(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 ϵ. Here, H_dw satisfies Eq. (9) with $\stackrel{\sim }{U}=\mathrm{diag}(1,1,-1,-1)$, meaning that H_dw is anti-𝒫𝒯 symmetric. This property stems from the reflection symmetry x ↦ −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, ϵ), which is a curve in a plane. Values of (k, ϵ) 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₁ ≡ ℛ(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-𝒫𝒯 symmetric, R₁ is a real quantity. The values of (k, ϵ) for which R₁ = 0 therefore define a curve. We calculate

$$\begin{array}{c}\hfill R_1=-4k^6{(ϵ-1)}^2(1+k^4-ϵ+k^2(2+9ϵ(-1+\frac{3}{4}ϵ))),\end{array}$$(38)

and show the curve R₁ = 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 ϵ > 8/9, which is indeed the minimum of the curve of EPs.

Fig. 4. Same as Fig. 1 but for the DUSTYWAVE problem. For k = 0.4, the Krein quantity $X·\stackrel{\sim }{V}X$ is represented as a function of ϵ. This Krein quantity provides a quadratic constant of motion for dusty sound waves in the propagating phase (low ϵ).

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, ϵ) space. These two points are located at

$$\begin{array}{c}\hfill k=\pm \frac{1}{\sqrt{3}},ϵ=\frac{8}{9},\end{array}$$(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 $\frac{-i}{3}$. We note that these EP3s carry topological charges, also called winding numbers W₃:

$$\begin{array}{c}\hfill W_3{|}_{{\mathrm{EP}3}_1}=+1,W_3{|}_{{\mathrm{EP}3}_2}=-1,\end{array}$$(40)

where EP3₁ is the point $(k,ϵ)=(+\frac{1}{\sqrt{3}},\frac{8}{9})$ and EP3₂ is the point $(k,ϵ)=(-\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 $\stackrel{\sim }{V}$ can be found explicitly

$$\begin{array}{cc}\hfill \stackrel{\sim }{V}& =\left(\begin{array}{cccc}ϵ& ϵ-1& \mathit{ik}& 0\\ ϵ-1& \frac{{(ϵ-1)}^2}{ϵ}& 0& -\frac{ik(ϵ-1)}{ϵ}\\ -ik& 0& 0& 0\\ 0& \frac{ik(ϵ-1)}{ϵ}& 0& 0\end{array}\right).\hfill \end{array}$$(41)

As has been shown, pseudo-chirality implies that if ω ≠ −ω^*, the quantity $X·\stackrel{\sim }{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·\stackrel{\sim }{V}X:$

$$\begin{array}{cc}\hfill X·\stackrel{\sim }{V}X& =\frac{1-ϵ}{ϵ}|{\rho }_{\mathrm{d}}{|}^2+\frac{ϵ}{1-ϵ}{|{\rho }_{\mathrm{g}}|}^2\hfill \\ \hfill & -2(\mathrm{Re}({\rho }_{\mathrm{g}}{\rho }_{\mathrm{d}})-k\frac{{\rho }_0}{c_{\mathrm{s}}}\mathrm{Im}({\rho }_{\mathrm{g}}{v}_{\mathrm{g}}^{\ast }+{\rho }_{\mathrm{d}}{v}_{\mathrm{d}}^{\ast })),\hfill \end{array}$$(42)

where ρ₀ = ρ_g, 0 + ρ_d, 0 is the total background density. $X·\stackrel{\sim }{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·\stackrel{\sim }{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-𝒫𝒯 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 ρ = ρ_g + ρ_d, v = ρ_gv_g + ρ_dv_d, $\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}$ and Δ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

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill X·{\stackrel{\sim }{V}}_{\mathrm{mono}}X& =\frac{{(1-ϵ)}^3}{ϵ}{\left|{\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\right)}^{\prime }\right|}^2\hfill \\ \hfill & +2k\mathrm{Im}({(1-ϵ)}^2{\left(\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\right)}^{\prime }\mathrm{\Delta }{v}^{\prime \ast }+\frac{{\rho }^{\prime }}{{\rho }_0}{v}^{\prime \ast }).\hfill \end{array}\end{array}$$

This quantity is associated with the symmetry operator

$$\begin{array}{c}\hfill {\stackrel{\sim }{V}}_{\mathrm{mono}}=\left(\begin{array}{cccc}0& 0& -ik& 0\\ 0& \frac{{(ϵ-1)}^3}{ϵ}& 0& -ik{(ϵ-1)}^2\\ ik& 0& 0& 0\\ 0& ik{(ϵ-1)}^2& 0& 0\end{array}\right).\end{array}$$(43)

Squire & Hopkins (2018) show that when grains are additionally streaming though the gas, a generic resonant drag instability develops in the mixture. The system is neither 𝒫𝒯 nor anti-𝒫𝒯 symmetric and Krein theory cannot be applied directly. The instability comes from a more complex mechanism where a resonance occurs between an anti-𝒫𝒯 symmetric correction due to drag and a leading 𝒫𝒯 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

$$\begin{array}{cc}& H_{\mathrm{sgd}}=\hfill \\ \hfill & \left(\begin{array}{cccccc}0& ik{\mathrm{\Sigma }}_g& 0& 0& 0& 0\\ -ik(\frac{2\pi \mathcal{G}}{|k|}-\frac{c_g^2}{{\mathrm{\Sigma }}_g})& \frac{ϵ}{t_s}& -2\mathrm{\Omega }& -\frac{2i\pi \mathcal{G}k}{|k|}& -\frac{ϵ}{t_s}& 0\\ 0& -2B& \frac{ϵ}{t_s}& 0& 0& -\frac{ϵ}{t_s}\\ 0& 0& 0& 0& ik{\mathrm{\Sigma }}_d& 0\\ -\frac{2i\pi \mathcal{G}k}{|k|}& -\frac{1}{t_s}& 0& -ik(\frac{2\pi \mathcal{G}}{|k|}-\frac{c_d^2}{{\mathrm{\Sigma }}_d})& \frac{1}{t_s}& -2\mathrm{\Omega }\\ 0& 0& -\frac{1}{t_s}& 0& -2B& \frac{1}{t_s}\end{array}\right),\hfill \end{array}$$(44)

where Ω 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_gt_s. The characteristic polynomial of the system is of order 5, meaning that the roots ω(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 𝒫𝒯 symmetric, with the operator U = diag(−1, 1, 1, −1, 1, 1). The marginal stability curve is therefore a curve of EPs that corresponds to the spontaneous breaking of 𝒫𝒯 symmetry. From the analysis above, its exact expression is

$$\begin{array}{c}\hfill {\mathcal{R}}_{\mathrm{sgd}}(P,P^{\prime })=0,\end{array}$$(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 𝒫𝒯-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).

Fig. 5. Exceptional Points of the dusty self-gravitating disk are a consequence of the spontaneous breaking of 𝒫𝒯 symmetry. The equation R(P,P′) = 0 is obtained analytically with a resultant method, reproducing the results of Longarini et al. (2022).

8. Conclusion

For several astrophysical objects, the properties of small perturbations around equilibrium are controlled by discrete symmetries, such as the 𝒫𝒯 symmetry and the anti-𝒫𝒯 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-𝒫𝒯 symmetric because of diffusion, and is also anti-Hermitian.

• The propagation of a wave in a dust–gas mixture: the DUSTYWAVE problem is anti-𝒫𝒯 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).

Acknowledgments

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

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 = −gu_z, and has constant density ρ₀. The temperature profile has a constant gradient T(z) = T_h − βz with β = (T_h − T_c)/L. The coefficient of thermal expansion of the fluid is $\alpha =-{\rho }_0^{-1}{\left(\partial \rho /\partial T\right)}_P$. The dynamical viscosity and thermal diffusivity of the fluid are denoted by ν and κ, respectively. The linear system describing the evolution of small perturbations in the Boussinesq approximation around this background state is

$$\begin{array}{cc}\hfill {\partial }_x{v}_x+{\partial }_z{v}_z& =0,\hfill \end{array}$$(A.1)

$$\begin{array}{cc}\hfill ({\partial }_t-\nu \mathrm{\Delta })v_x& =-\frac{1}{{\rho }_0}{\partial }_{x}p,\hfill \end{array}$$(A.2)

$$\begin{array}{cc}\hfill ({\partial }_t-\nu \mathrm{\Delta })v_z& =-\frac{1}{{\rho }_0}{\partial }_{z}p+\alpha gT,\hfill \end{array}$$(A.3)

$$\begin{array}{cc}\hfill ({\partial }_t-\kappa \mathrm{\Delta })T& =\beta v_z.\hfill \end{array}$$(A.4)

Applying double curl, and enforcing the incompressibility equation ∇ ⋅ v = 0, the equation for u_x decouples from the other two equations, leaving

$$\begin{array}{cc}\hfill ({\partial }_t-\nu \mathrm{\Delta })\mathrm{\Delta }{v}_z& =\alpha g{\partial }_{\mathit{xx}}T,\hfill \end{array}$$(A.5)

$$\begin{array}{cc}\hfill ({\partial }_t-\kappa \mathrm{\Delta })T& =\beta v_z.\hfill \end{array}$$(A.6)

For plane wave solutions of the form exp( − iωt + ik_zz + ik_xx), Eqs. A.5 – A.6 become

$$\begin{array}{cc}\hfill (-i\omega -\nu (k_x^2+k_z^2))(k_x^2+k_z^2)v_z& =-\alpha g{k}_x^{2}T,\hfill \end{array}$$(A.7)

$$\begin{array}{cc}\hfill (-i\omega -\kappa (k_x^2+k_z^2))T& =\beta v_z.\hfill \end{array}$$(A.8)

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

$$\begin{array}{cc}\hfill \omega X& =H_{\mathrm{RB}}X,\hfill \end{array}$$(A.9)

$$\begin{array}{cc}\hfill H_{\mathrm{RB}}& =-i\left(\begin{array}{cc}{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{array}\right),\hfill \end{array}$$(A.10)

which corresponds to Eq.32. Consequently, iH_RB is Hermitian if and only if Ra ≥ 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

$$\begin{array}{cc}\hfill ({\partial }_t+v_{\mathrm{g}}{\partial }_x){\rho }_{\mathrm{g}}& =-{\rho }_{\mathrm{g}}{\partial }_x{v}_{\mathrm{g}},\hfill \end{array}$$(B.1)

$$\begin{array}{cc}\hfill ({\partial }_t+v_{\mathrm{d}}{\partial }_x){\rho }_{\mathrm{d}}& =-{\rho }_{\mathrm{d}}{\partial }_x{v}_{\mathrm{d}},\hfill \end{array}$$(B.2)

$$\begin{array}{cc}\hfill ({\partial }_t+v_{\mathrm{g}}{\partial }_x)v_{\mathrm{g}}& =\frac{K}{{\rho }_{\mathrm{g}}}(v_{\mathrm{d}}-v_{\mathrm{g}})-\frac{1}{{\rho }_{\mathrm{g}}}{\partial }_{x}P,\hfill \end{array}$$(B.3)

$$\begin{array}{cc}\hfill ({\partial }_t+v_{\mathrm{d}}{\partial }_x)v_{\mathrm{d}}& =\frac{K}{{\rho }_{\mathrm{d}}}(v_{\mathrm{g}}-v_{\mathrm{d}}),\hfill \end{array}$$(B.4)

where ρ_g, ρ_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₀ + f′. The system of equations can be put in a symmetric form using the rescaled quantities $({\rho }_{\mathrm{g}}^{\prime },{\rho }_{\mathrm{d}}^{\prime },v_{\mathrm{g}}^{\prime },v_{\mathrm{d}}^{\prime })\mapsto X={(\frac{1}{{\rho }_{\mathrm{g},0}+{\rho }_{\mathrm{d},0}}{\rho }_{\mathrm{g}}^{\prime },\frac{1}{{\rho }_{\mathrm{g},0}+{\rho }_{\mathrm{d},0}}{\rho }_{\mathrm{d}}^{\prime },\frac{{\rho }_{\mathrm{g},0}}{c_{\mathrm{s}}({\rho }_{\mathrm{g},0}+{\rho }_{\mathrm{d},0})}{v}_{\mathrm{g}}^{\prime },\frac{{\rho }_{\mathrm{d},0}}{c_{\mathrm{s}}({\rho }_{\mathrm{g},0}+{\rho }_{\mathrm{d},0})}{v}_{\mathrm{d}}^{\prime })}^{\top }$. For isothermal pertubations $P^{\prime }(x,t)=c_{\text{s}}^2{\rho }_{\text{g}}^{\prime }(x,t)$, the evolution of the perturbations are given by the system

$$\begin{array}{c}\hfill {\partial }_{t}X=\left(\begin{array}{cccc}0& 0& -c_{\mathrm{s}}{\partial }_x& 0\\ 0& 0& 0& -c_{\mathrm{s}}{\partial }_x\\ -c_{\mathrm{s}}{\partial }_x& 0& -\frac{K}{{\rho }_{\mathrm{g},0}}& +\frac{(1-ϵ)K}{ϵ{\rho }_{\mathrm{g},0}}\\ 0& 0& \frac{ϵK}{(1-ϵ){\rho }_{\mathrm{d},0}}& -\frac{K}{{\rho }_{\mathrm{d},0}}\end{array}\right)X,\end{array}$$(B.5)

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

$$\begin{array}{cc}\hfill H_{\mathrm{dw}}X& =i{\partial }_{t}X,\hfill \end{array}$$(B.6)

$$\begin{array}{ccc}\hfill H_{\mathrm{dw}}& \equiv \hfill & \hfill \left(\begin{array}{cccc}0& 0& k& 0\\ 0& 0& 0& k\\ k& 0& iϵ& -i(1-ϵ)\\ 0& 0& -iϵ& i(1-ϵ)\end{array}\right).\end{array}$$(B.7)

Appendix C: EPs from resultants

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

$$\begin{array}{cc}\hfill Q_1\left(X\right)& ={\displaystyle \sum _{i=1}^n}{a}_i{X}^i=a_n{\displaystyle \prod _{i=1}^n}(X-{\alpha }_i),\hfill \end{array}$$(C.1)

$$\begin{array}{cc}\hfill Q_2\left(X\right)& ={\displaystyle \sum _{j=1}^m}{b}_j{X}^j=b_m{\displaystyle \prod _{j=1}^m}(X-{\beta }_j),\hfill \end{array}$$(C.2)

where α_i and β_j are their roots. The resultant ℛ(Q₁, Q₂) of these two polynomials is defined as

$$\begin{array}{c}\hfill \mathcal{R}(Q_1,Q_2)\equiv detS(Q_1,Q_2),\end{array}$$(C.3)

where S(Q₁, Q₂) is the Sylvester matrix of Q₁ and Q₂. It satisfies

$$\begin{array}{c}\hfill \mathcal{R}(Q_1,Q_2)=a_n^m{b}_m^n{\displaystyle \prod _{i,j}}({\alpha }_i-{\beta }_j),\end{array}$$(C.4)

which cancels if and only if Q₁ and Q₂ 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 α_i, β_j. We can define P′, the derivative of P(X). The cancellation of the resultant R₁ 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

$$\begin{array}{c}\hfill R_1\equiv \mathcal{R}(P,P^{\prime })=0.\end{array}$$(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, ϵ) space. These points can be located using a procedure that generalises the determination of EP2s. Noting R₁ = ℛ(P, P′) and 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₁ and R₂ cancel simultaneously. For the DUSTYWAVE problem, R₁ is given by Eq. 38, and R₂ by

$$\begin{array}{c}\hfill R_2=400k^8-432k^4(ϵ-1)+96k^6(14-51ϵ+36{ϵ}^2).\end{array}$$(D.1)

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

$$\begin{array}{c}\hfill k=\pm \frac{1}{\sqrt{3}},ϵ=\frac{8}{9}.\end{array}$$(D.2)

Fig. D.1. Exceptional points of the DUSTYWAVE in the parameter space (k, ϵ). Intersections between the red (R1 = 0) and blue (R2 = 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 ±1.

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{array}{c}\hfill \left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}-i\sqrt{3}/2\\ 2i\sqrt{3}\\ -1/2\\ 2\end{array}\right).\end{array}$$(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

$$\begin{array}{c}\hfill W_3=\frac{1}{2\pi }{\displaystyle {\oint }_C}\frac{1}{R_1^2+R_2^2}(R_1{\mathrm{\nabla }}_{\mathit{p}}{R}_2-R_2{\mathrm{\nabla }}_{\mathit{p}}{R}_1)\mathrm{d}\mathit{p},\end{array}$$(D.4)

for ${\mathrm{\nabla }}_{\mathit{p}}=\left(\begin{array}{c}{\partial }_k\\ {\partial }_{ϵ}\end{array}\right)$ and $\mathit{p}=\left(\begin{array}{c}k\\ ϵ\end{array}\right)\in C,$ where C is any closed loop around one of the EP3s. We compute

$$\begin{array}{cc}\hfill W_3{|}_{{\mathrm{EP}3}_1}& =+1,\hfill \end{array}$$(D.5)

$$\begin{array}{cc}\hfill W_3{|}_{{\mathrm{EP}3}_2}& =-1,\hfill \end{array}$$(D.6)

where EP3₁ is the point at $(k,ϵ)=(+\frac{1}{\sqrt{3}},\frac{8}{9})$ and EP3₂ is the point at $(k,ϵ)=(-\frac{1}{\sqrt{3}},\frac{8}{9})$.

Appendix E: DUSTYWAVE: Anti-𝒫𝒯 operator

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

$$\begin{array}{c}\hfill H_{\mathrm{dw}}\mathcal{CP}+\mathcal{CP}{H}_{\mathrm{dw}}=0.\end{array}$$(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₀ an eigenvector of H_dw with pure imaginary eigenvalue. Formally,

$$\begin{array}{c}\hfill H_{\mathrm{dw}}{X}_0=i{\eta }_0{X}_0,\end{array}$$(E.2)

where η₀ is real. Applying the operator 𝒞𝒫 on Eq. E.2 and the commutation rule Eq.(E.1) gives the relation

$$\begin{array}{c}\hfill H_{\mathrm{dw}}(\mathcal{CP}{X}_0)=i{\eta }_0(\mathcal{CP}{X}_0).\end{array}$$(E.3)

If the eigenvalue is not degenerated (i.e. the parameters are chosen for the system not to be at an EP), 𝒞𝒫X₀ is colinear to X₀, i.e.

$$\begin{array}{c}\hfill \mathcal{CP}{X}_0=\alpha X_0,\end{array}$$(E.4)

for some complex number α. Hence, H_dw and 𝒞𝒫 share the same eigenvectors and the anti-𝒫𝒯 symmetry is said to be unbroken. We note that the same reasoning applied to the operator iH_dw shows that iH_dw is 𝒫𝒯 symmetric. 𝒞𝒫 can be straightforwardly diagonalised. Its eigenvectors are

$$\begin{array}{c}\hfill \left(\begin{array}{c}1\\ 0\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 1\\ 0\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 0\\ i\\ 0\end{array}\right),\left(\begin{array}{c}0\\ 0\\ 0\\ i\end{array}\right).\end{array}$$(E.5)