Exceptional points in 1D arrays of quantum harmonic oscillators^(a)

Let us introduce the vector of moments via the orthogonal transformation $\widetilde{\vec{t}{u}}=\mathsf{P}^T\vec{t}{u}$, where $\widetilde{\vec{t}{u}}$ contains the k-modes $\hat{a}_{k_l}, \hat{b}_{k_l}$ defined by

$$\begin{eqnarray} \hat{a}_{k_l}&= \sum_{n=1}^N\hat{a}_n \mathcal{S}^{(a)}_{n,k_l},\quad \hat{a}_n=\sum_{l=1}^N\hat{a}_{k_l}\mathcal{S}^{(a)}_{n,k_l}, \end{eqnarray} \tag{ 5 }$$

$$\begin{eqnarray} \mathcal{S}^{(a)}_{n,k_l}&= \sqrt{\frac{2}{N+\frac{1}{2}}}\sin\bigg[k_l\bigg(n-\frac{1}{2}\bigg)\bigg],\nonumber\\ \mathcal{S}^{(b)}_{n,k_l}&= \sqrt{\frac{2}{N+\frac{1}{2}}}\sin(k_ln),\nonumber \end{eqnarray}$$

with $k_l={\pi l}/{(N+1/2)}$, $l=1,\ldots,N$ and similarly for $\hat{b}_n, \hat{b}_{k_l}$. This transformation preserves the commutation relations because $\sum_{n=1}^N \mathcal{S}^{(x)}_{n,k_l}\mathcal{S}^{(x)}_{n,k_{l'}}=\delta_{l,l'}$ and $\sum_{l=1}^N \mathcal{S}^{(x)}_{n,k_l}\mathcal{S}^{(x)}_{n',k_{l}}=\delta_{n,n'}$ ($x=a,b$), as one can readily show. Then $\widetilde{\mathsf{M}}=\mathsf{P}^T\mathsf{M}\mathsf{P}$ with $\widetilde{\mathsf{M}}=\bigoplus_{l=1}^N \widetilde{\mathsf{M}}_{k_l}$ and

$$\begin{equation} \widetilde{\mathsf{M}}_{k_l}=\left(\begin{array}{c@{\quad}c} -\gamma_a & -i\lambda_{k_l}\\ -i\lambda_{k_l} & -\gamma_b\\ \end{array}\right)\!, \end{equation} \tag{ 6 }$$

with $\lambda_{k_l}=2\lambda \cos(k_l/2)$. The eigenvalues of the whole system are obtained by diagonalizing the $2\times 2$ blocks and read

$$\begin{equation} \Omega^\pm_{k_l}=\frac{\bar{\gamma}}{2}\pm\frac{1}{2}\sqrt{(\gamma_a-\gamma_b)^2-4\lambda_{k_l}^2}, \end{equation} \tag{ 7 }$$

with $\bar{\gamma}=\gamma_a+\gamma_b$. We denote the eigenvalues real and imaginary parts as $\gamma^\pm_{k_l}$ and $\nu^\pm_{k_l}$, respectively (fig. 1(a)). We obtain that, as the coupling λ is varied, below a certain threshold the eigenvalues are real, while above they become complex. When the imaginary parts vanish (moving from the rotating frame to the lab one) all units oscillate at their bare frequencies, $\omega_0$, in spite of their mutual interactions. This is in stark contrast to what one observes in the corresponding closed system where, for any coupling strength, the eigenfrequencies differ from the intrinsic frequency $\omega_0$. At the symmetry breaking threshold point the square root in eq. (7) vanishes and the two eigenvalues with the corresponding eigenvectors coalesce [1,2]. This point is thus an EP, at which $\widetilde{\mathsf{M}}_{k_l}$ (and thus $\mathsf{M}$) becomes non-diagonalizable. Therefore, our system displays N EPs (fig. 1(a)), one for each k_l. The EPs separate different regimes in which the number of frequencies varies from 1 to 2N, a feature that cannot be observed for any closed or homogeneously dissipating system. Indeed, in the case of balanced loss rates $\gamma_a=\gamma_b$, the bath would only introduce a uniform dissipation decay, while the frequencies would depart from $\omega_0$ and be modified as in the absence of the environment.

To study the decay dynamics near an EP, we proceed to solve directly the block-diagonal system of equations using the Laplace transform method. The equations for the k-modes are

$$\begin{eqnarray} \frac{\textrm{d}}{\textrm{d}t}\langle{\hat{a}}_{k_l}\rangle&= -\gamma_a \langle\hat{a}_{k_l}\rangle -i\lambda_{k_l} \langle\hat{b}_{k_l}\rangle, \end{eqnarray} \tag{ 8 }$$

$$\begin{eqnarray} \frac{\textrm{d}}{\textrm{d}t}\langle{\hat{b}}_{k_l}\rangle&= -\gamma_b \langle\hat{b}_{k_l}\rangle -i\lambda_{k_l} \langle\hat{a}_{k_l}\rangle. \end{eqnarray} \tag{ 9 }$$

Solving their Laplace transformed version we obtain in the s-domain:

$$\begin{eqnarray} \langle \hat{a}_{k_l}(s)\rangle&= \frac{(s+\gamma_b)\langle\hat{a}_{k_l}\rangle_0-i\lambda_{k_l} \langle\hat{b}_{k_l}\rangle_0}{(s+\Omega^+_{k_l})(s+\Omega^-_{k_l})}, \end{eqnarray} \tag{ 10 }$$

$$\begin{eqnarray} \langle\hat{b}_{k_l}(s)\rangle&= \frac{(s+\gamma_a)\langle\hat{b}_{k_l}\rangle_0-i\lambda_{k_l} \langle\hat{a}_{k_l}\rangle_0}{(s+\Omega^+_{k_l})(s+\Omega^-_{k_l})}, \end{eqnarray} \tag{ 11 }$$

and the initial condition given by $\langle\hat{a}_{k_l}(\hat{b}_{k_l})\rangle_0$. The poles $\Omega^\pm_{k_l}$ correspond to the eigenvalues of $\mathsf{M}$, which display an EP in each k-block at $|2\lambda_k|=|\gamma_a-\gamma_b|$. Here, as well as in the scattering matrix, an EP appear as a double pole, which makes clear the emergence of polynomial deviations from the typical exponential-decay behavior. From eqs. (7), (10) and (11) we identify two dynamical regimes for each k-mode, separated by an EP. i) Single-frequency regime $|\gamma_a-\gamma_b|>|2\lambda_{k_l}|$: we have two different real poles and no oscillations (in the rotating frame). ii) Exceptional point $|\gamma_a-\gamma_b|=|2\lambda_{k_l}|$: we have a real second-order pole, yielding polynomial corrections to the exponential decay without oscillations. iii) Two-frequency regime $|\gamma_a-\gamma_b|<|2\lambda_{k_l}|$: we have two different complex poles, yielding oscillations with equal damping. In the time domain the dynamics in each regime read as follows.

i) Single-frequency regime:

$$\begin{eqnarray} \langle \hat{a}_{k_l}\rangle=\,& \bigg\{\langle\hat{a}_{k_l}\rangle_0\frac{(\gamma_b-\gamma_{k_l}^+)e^{-\gamma_{k_l}^+ t}-(\gamma_b-\gamma_{k_l}^-)e^{-\gamma_{k_l}^- t}}{\gamma_{k_l}^- -\gamma_{k_l}^+}\nonumber\\ & {}-i\lambda_{k_l} \langle\hat{b}_{k_l}\rangle_0\frac{e^{-\gamma_{k_l}^+ t}-e^{-\gamma_{k_l}^- t}}{\gamma_{k_l}^- -\gamma_{k_l}^+}\bigg\} \theta(t), \end{eqnarray} \tag{ 12 }$$

where $\theta(t)$ is the Heaviside step function. Notice that in all cases the solution for $\langle \hat{b}_{k_l}\rangle$ is obtained from the one of $\langle \hat{a}_{k_l}\rangle$ by exchanging the labels $a\leftrightarrow b$ and the initial conditions. In this regime the poles are real and read as $2\gamma^\pm_{k_l}=\bar{\gamma}\pm\sqrt{(\gamma_a-\gamma_b)^2-4\lambda_{k_l}^2}$.

ii) Exceptional point:

$$\begin{equation} \langle \hat{a}_{k_l}\rangle=e^{-\frac{\bar{\gamma}}{2}t}\bigg\{\langle\hat{a}_{k_l}\rangle_0\bigg[1+\frac{\gamma_b-\gamma_a}{2}t\bigg]-i\lambda_{k_l} t \langle\hat{b}_{k_l}\rangle_0 \bigg\}\theta(t), \end{equation} \tag{ 13 }$$

iii) Two-frequency regime:

$$\begin{eqnarray} \langle \hat{a}_{k_l}\rangle=\,& e^{-\frac{\bar{\gamma}}{2}t}\bigg\{\langle\hat{a}_{k_l}\rangle_0\bigg[\frac{\gamma_b-\gamma_a}{2\nu_{k_l}}\sin(\nu_{k_l}t)+\cos(\nu_{k_l}t) \bigg]\nonumber\\ & {}-i\frac{\lambda_{k_l}}{\nu_{k_l}} \langle\hat{b}_{k_l}\rangle_0 \sin(\nu_{k_l}t) \bigg\}\theta(t), \end{eqnarray} \tag{ 14 }$$

Here the poles are complex, and the imaginary part is $\nu_{k_l}=\frac{1}{2}\sqrt{4\lambda_{k_l}^2-(\gamma_a-\gamma_b)^2}$.

The solution in the site basis is obtained reversing the orthogonal transformation, i.e., $\vec{t}{u}=\mathsf{P}\widetilde{\vec{t}{u}}$ (eq. (5)). Hence, in the site basis, the local dynamics is a combination of the different k-modes ones. As anticipated, the number of frequencies at which the modes can oscillate change depending on the regime. Below the first EP, the modes oscillate mono-chromatically, despite the finite coherent coupling. This region of monochromatic evolution of all units occurs in a parameter window that depends on the size of the system: in the limit of large N the decay contrast $|\gamma_a-\gamma_b|>4\lambda$ is required to overcome coupling, while for just a pair of oscillators one would get a smaller contrast $|\gamma_a-\gamma_b|>2\lambda$.

When the previous condition is not fulfilled, for instance increasing the coupling λ between elements, more k-modes enter in the two-frequency regime and, in the site basis, up to 2N frequencies are observed (fig. 1). At the EP, the first moments present an exponential-power-law decay [21]. The exponent of the polynomial corrections can be at most h  − 1, where h is the number of eigenvalues that coalesce, i.e., the order of the EP. For higher-order EPs, more frequencies merge at the EP, and the variability in the number of frequencies displayed by the system is magnified. This could occur in our system when periodically modulating losses over larger cells. Finally notice that, despite the peculiar behavior at the EP and its singular spectral character, the solutions are continuous in the transition from one regime to the other, as also observed in other physical contexts (see below).

Let us now move to the case of a system with driving, i.e., non-zero $\vec{t}{F}$. The details on the general solutions are provided in the supplementary material Supplementarymaterial.pdf (SM). Here we focus on the particular case in which we only drive one mode with an oscillatory input of frequency $\omega_n$, i.e., $\alpha_n(t)=\eta_c\alpha e^{i\Delta_n t}\theta(t)$ with $\Delta_n=\omega_0-\omega_n$. Then in the k-basis the driving terms read as $\alpha_{k_l}(t)=\eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}e^{i\Delta_n t}\theta(t)$ and $\beta_{k_l}(t)=0~\forall\,l$, where α is the driving strength and $\eta_c$ characterizes the input channel. Assuming vanishing initial conditions (SM), the solution in the s-domain is given by

$$\begin{eqnarray} \langle \hat{a}_{k_l}(s)\rangle&= \frac{(s+\gamma_b)\eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}}{(s-i\Delta_n)(s+\Omega^+_{k_l})(s+\Omega^-_{k_l})}, \end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray} \langle \hat{b}_{k_l}(s)\rangle&= \frac{-i\lambda_{k_l} \eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}}{(s-i\Delta_n)(s+\Omega^+_{k_l})(s+\Omega^-_{k_l})}. \end{eqnarray} \tag{ 16 }$$

Notice that the only pure imaginary pole is $s=i\Delta_n$. Hence, in the transient dynamics towards the steady state, the different dynamical regimes are manifested in the dynamical response of the system to the oscillatory drive through $\Omega_{k_l}^\pm$. In fact, in the steady state, the intrinsic dynamics of the system appears through the residue of $s=i\Delta_n$, which leads to the following stationary amplitudes for the k-modes:

$$\begin{eqnarray} \langle \hat{a}_{k_l}\rangle_{ss}&= \frac{(i\Delta_n+\gamma_b)\eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}e^{i\Delta_n t}}{(i\Delta_n+\Omega^+_{k_l})(i\Delta_n+\Omega^-_{k_l})}, \end{eqnarray} \tag{ 17 }$$

$$\begin{eqnarray} \langle \hat{b}_{k_l}\rangle_{ss}&= \frac{-i\lambda_{k_l}\eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}e^{i\Delta_n t}}{(i\Delta_n+\Omega^+_{k_l})(i\Delta_n+\Omega^-_{k_l})}. \end{eqnarray} \tag{ 18 }$$

It is useful for the discussion that follows to particularize the expression for $\langle \hat{a}_{k_l}\rangle_{ss}$ in the two dynamical regimes, making use of the explicit expression of $\Omega^\pm_{k_l}$. In the single-frequency regime we have

$$\begin{equation} \langle \hat{a}_{k_l}\rangle_{ss}=\frac{\eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}e^{i\Delta_n t}}{\gamma^+_{k_l}-\gamma^-_{k_l}}\bigg[\frac{\gamma^+_{k_l}-\gamma_b}{i\Delta_n+\gamma^+_{k_l}} -\frac{\gamma^-_{k_l}-\gamma_b}{i\Delta_n+\gamma^-_{k_l}} \bigg], \end{equation} \tag{ 19 }$$

while in the two-frequency regime we have

$$\begin{eqnarray} \langle \hat{a}_{k_l}\rangle_{ss}=\,& \frac{\eta_c\alpha\mathcal{S}_{n,k_l}^{(a)}e^{i\Delta_n t}}{2i\nu_{k_l}}\bigg[\frac{\frac{\gamma_a-\gamma_b}{2}+i\nu_{k_l}}{i(\Delta_n+\nu_{k_l})+\frac{\bar{\gamma}}{2}}\nonumber\\ & {}-\frac{\frac{\gamma_a-\gamma_b}{2}-i\nu_{k_l}}{i(\Delta_n-\nu_{k_l})+\frac{\bar{\gamma}}{2}} \bigg]. \end{eqnarray} \tag{ 20 }$$

We can find signatures of these dynamical regimes in the transmission spectrum of the driven mode. Within the assumption that the system is driven through the environment described in eq. (2), the input-output relation for the driven mode is just $\hat{a}^{\textrm{out}}_n=\hat{a}^{\textrm{in}}_n-\eta_c\hat{a}_n$, where $\langle \hat{a}^{\textrm{in}}_n\rangle =\alpha e^{i\Delta_n t}\theta(t)$ and $\eta_c=\sqrt{\gamma_a}$ [47]. The transmission coefficient is defined as $T=\langle\hat{a}^{\textrm{out}}_n\rangle/\langle\hat{a}^{\textrm{in}}_n\rangle$, which in the steady state reads

$$\begin{equation} T_{ss}=1 -\sqrt{\gamma_a}\sum_{l=1}^N \frac{ \mathcal{S}_{n,k_l}^{(a)}\langle \hat{a}_{k_l}\rangle_{ss}}{\langle\hat{a}^{\textrm{in}}_n\rangle}. \end{equation} \tag{ 21 }$$

It is clear that the different regimes described explicitly by eqs. (19) and (20) translate into different behaviors of the transmission spectrum. This can be observed in fig. 2, in which we plot eq. (21) (red solid lines) for three different cases: in (a) and (b) $\lambda/\omega_0=0.02$ all k-modes are in the single-frequency regime, in (c) $\lambda/\omega_0=0.1$ (only the modes of k₃ are in the single-frequency regime), while in (d) $\lambda/\omega_0=0.35$ (all k-modes are well into the two-frequency regime).

Zoom In Zoom Out Reset image size

A surprising feature of this spectrum is that, in cases (a), (b), (c), the system displays a transparency window centered at $\omega_0$, despite the presence of k-modes of this frequency that naively should resonate perfectly with the driving. This can be understood by analyzing the stationary amplitude of the k-modes in the single-frequency regime. Indeed, from eq. (19) we see that these amplitudes have two contributions that interfere destructively, both centered at the same frequency but with different decay rate. In fact, for weak coupling $\gamma_a\gg\lambda_{k_l}$, large disparity $\gamma_a\gg\gamma_b$, and near the resonant driving $\gamma_a\gg|\Delta_n|$, we can approximate the transmission spectrum as

$$\begin{equation} T_{ss}\approx \sum_{l=1}^N \frac{(\lambda_{k_l}^2/\gamma_a) (\mathcal{S}_{n,k_l}^{(a)})^2}{i\Delta_n+\gamma_b+\lambda_{k_l}^2/\gamma_a}, \end{equation} \tag{ 22 }$$

where we have made the replacements $\gamma_{k_l}^+\approx\gamma_a$, $\gamma_{k_l}^-\approx \gamma_b+\lambda_{k_l}^2/\gamma_a$ and $\gamma_a+i\Delta_n\approx \gamma_a$, accordingly to the above approximations. In fig. 2(a) and (b) we compare this approximate expression (blue dashed lines) with the exact one, finding good agreement in the expected region. Thus, eq. (22) enables us to clearly visualize the Lorentzian transparency window that emerges in the single-frequency regime. Notice that the height and width of this window increase with the coupling strength, as is clear from this expression and panels (a) and (c). The interaction-induced transparency windows have been already observed in different systems and configurations, as, for instance, in optical resonators [27,28] or optomechanical systems [29–31], and can be ultimately traced back to the rather general phenomena of Fano interference and EIT [45,48].

In the two-frequency regime we find the more intuitive result that each k-mode contributes to the transmission spectrum with a Lorentzian absorption line centered at its frequency (see eq. (20) and fig. 2(d)). In this case it is no surprise that the system is fully transparent at $\Delta_1=0$ as there are no k-modes resonant to this frequency. Furthermore, a distinctive feature of our system is that it can work in both regimes simultaneously, as in case of fig. 2(c), in which there is the transparency window due to the k₃-modes, but there are also two double-peak Lorentzian absorption lines due to the other k-modes. Finally we observe that the transmission coefficient is continuous along the different dynamical regimes, as follows from eqs. (17) and (21) and the fact that $\Omega^\pm_{k_l}$ are continuous. This is in spite of the singular character of the EPs and in accordance with what we have commented previously.

The system of differential equations describing the dynamics of the second moments is divided into two uncoupled blocks: a homogeneous set of equations for the dynamics of all possible moments of the kind $\langle \hat{a}^\dagger_n \hat{a}^\dagger_{n'}\rangle$, $\langle \hat{b}^\dagger_n \hat{b}^\dagger_{n'}\rangle$, $\langle \hat{a}^\dagger_n \hat{b}^\dagger_{n'}\rangle$ and their Hermitian conjugates, and an inhomogeneous system of equations set for the dynamics of all terms of type $\langle \hat{a}^\dagger_n \hat{a}_{n'}\rangle$, $\langle \hat{b}^\dagger_n \hat{b}_{n'}\rangle$, $\langle \hat{a}^\dagger_n \hat{b}_{n'}\rangle$ together with their conjugates. The inhomogeneous block of equations (the only relevant part in our analysis) can be readily written from eqs. (1) and (2) (as done in the SM) and its structure is that of a square lattice with four different sites per unit cell, fig. 3(a). This suggests that the system of equations can be brought to a block diagonal form by means of the orthogonal transformation $\mathsf{P}$ generalized to a 2D lattice. Thus, two wave vectors are needed, k_x and k_y, and each block reads

$$\begin{eqnarray} \dot{A}_{k_x,k_y}=\,& i(\lambda_{k_x}D_{k_x,k_y}-\lambda_{k_y}C_{k_x,k_y})\nonumber\\ & {}-2\gamma_a A_{k_x,k_y}+2\gamma_a^h\delta_{k_x,k_y}, \end{eqnarray} \tag{ 23 }$$

$$\begin{eqnarray} \dot{B}_{k_x,k_y}=\,& i(\lambda_{k_x}C_{k_x,k_y}-\lambda_{k_y}D_{k_x,k_y})\nonumber\\ & {}-2\gamma_b B_{k_x,k_y}+2\gamma_b^h\delta_{k_x,k_y}, \end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray} \dot{C}_{k_x,k_y}&= i(\lambda_{k_x}B_{k_x,k_y}-\lambda_{k_y}A_{k_x,k_y}) -\bar{\gamma} C_{k_x,k_y}, \end{eqnarray} \tag{ 25 }$$

$$\begin{eqnarray} \dot{D}_{k_x,k_y}&= i(\lambda_{k_x}A_{k_x,k_y}-\lambda_{k_y}B_{k_x,k_y})-\bar{\gamma} D_{k_x,k_y}, \end{eqnarray} \tag{ 26 }$$

where we have defined $A_{k_x,k_y}=\langle\hat{a}^\dagger_{k_x}\hat{a}_{k_y} \rangle$, $B_{k_x,k_y}=\langle\hat{b}^\dagger_{k_x}\hat{b}_{k_y} \rangle$, $C_{k_x,k_y}=\langle\hat{a}^\dagger_{k_x}\hat{b}_{k_y} \rangle$, and $D_{k_x,k_y}=\langle\hat{b}^\dagger_{k_x}\hat{a}_{k_y} \rangle$ and we have used the dot in place of the time derivative. Hence, we have transformed a general set of 4N coupled differential equations into a reduced one of only 4 coupled equations, which can readily be solved. In the following we calculate explicitly the general solution of eqs. (23)–(26) in the cases in which EPs are observed and the stationary state of the system.

Zoom In Zoom Out Reset image size

Let us take an initial condition that is uniform along each sublattice, $\langle \hat{a}^\dagger_n \hat{a}_{m} \rangle_0=n_a\delta_{n,m}$, $\langle \hat{b}^\dagger_n \hat{b}_{m} \rangle_0=n_b\delta_{n,m}$ and $\langle \hat{b}^\dagger_n \hat{a}_{m} \rangle_0=\langle \hat{a}^\dagger_n \hat{b}_{m} \rangle_0=0$ $\forall\, n,m$, and analyze the evolution of the population of each mode. Using the orthogonality properties of the canonical transformation defined in eq. (5), we can show that in the k-basis the initial condition reads as $\langle \hat{a}^\dagger_{k_x} \hat{a}_{k_y} \rangle_0=n_a\delta_{k_x,k_y}$, $\langle \hat{b}^\dagger_{k_x} \hat{b}_{k_y} \rangle_0=n_b\delta_{k_x,k_y}$ and the rest of the second moments zero. Then all sectors of eqs. (23)–(26) with $k_x\neq k_y$ vanish identically at all times. The solutions in the k-basis and in s-domain are

$$\begin{eqnarray} A_{k_l,k_l}(s)&= \frac{[2\lambda^2_{k_l}\!+\!(s\!+\!2\gamma_b)(s\!+\!\bar{\gamma})]A(s)\!+\!2\lambda^2_{k_l}B(s)} {(s\!+\!\bar{\gamma})(s\!+\!2\Omega^+_{k_l})(s\!+\!2\Omega^-_{k_l})},\qquad \end{eqnarray} \tag{ 27 }$$

$$\begin{eqnarray} B_{k_l,k_l}(s)&= \frac{[2\lambda^2_{k_l}\!+\!(s\!+\!2\gamma_a)(s\!+\!\bar{\gamma})]B(s)\!+\!2\lambda^2_{k_l}A(s)} {(s\!+\!\bar{\gamma})(s\!+\!2\Omega^+_{k_l})(s\!+\!2\Omega^-_{k_l})}, \end{eqnarray} \tag{ 28 }$$

$$\begin{eqnarray} C_{k_l,k_l}(s)&= \frac{i\lambda_{k_l}[(s+2\gamma_a)B(s)-(s+2\gamma_b)A(s)]} {(s+\bar{\gamma})(s+2\Omega^+_{k_l})(s+2\Omega^-_{k_l})}, \end{eqnarray} \tag{ 29 }$$

with $D_{k_l,k_l}(s)=C^*_{k_l,k_l}(s)$, $A(s)=n_a+2\gamma^h_a/s$ and $B(s)=n_b+2\gamma^h_b/s$. As before, the dynamics of the second moments in the site basis can be obtained transforming back to time domain and writing $\langle\hat{a}^\dagger_n\hat{a}_m\rangle=\sum_{l=1}^N A_{k_l,k_l}S_{m,k_l}^{(a)}S_{n,k_l}^{(a)}$, $\langle\hat{b}^\dagger_n\hat{b}_m\rangle=\sum_{l=1}^N B_{k_l,k_l}S_{n,k_l}^{(b)}S_{m,k_l}^{(b)}$ and $\langle\hat{a}^\dagger_n\hat{b}_m\rangle=\sum_{l=1}^N C_{k_l,k_l}S_{n,k_l}^{(a)}S_{m,k_l}^{(b)}$. From eqs. (27)–(29) it is clear that $A_{k_l,k_l}$, $B_{k_l,k_l}$ and $C_{k_l,k_l}$, display the same dynamical regimes as the first moments, and, thus, the time evolution of the second moments presents also the characteristic features of each regime: only one frequency ($\omega_0$) and multiple decay rates in the single-frequency regime, polynomial corrections in the decay just at the EP, and multiple frequencies with the same decay rate in the two-frequency regime. Indeed, here at the EP three poles coalesce, which makes the polynomial corrections to be quadratic in time. Examples of the emergent dynamics in the site basis are plotted in fig. 3(b).

Two key observations are needed to find the stationary state: only the equations for the blocks with $k_x=k_y$ have sources, and all blocks described by a homogeneous system of equations go to zero in the stationary state. Then, setting the time derivatives of eqs. (23)–(26) to zero and solving the algebraic equations, we can obtain the non-zero moments in the stationary state. In this way the obtained stationary second moments recombined in the site basis read as

$$\begin{eqnarray} \langle \hat{a}^\dagger_n \hat{a}_m\rangle_{ss}&= \sum_{l=1}^N\bigg[ \frac{(\gamma_b\bar{\gamma}\!+\!\lambda_{k_l}^2)\gamma_a^h\!+\!\lambda_{k_l}^2\gamma_b^h}{\bar{\gamma}\lambda_{k_l}^2 \!+\!\gamma_a\gamma_b\bar{\gamma}}\bigg]\mathcal{S}^{(a)}_{n,k_l}\mathcal{S}^{(a)}_{m,k_l},\qquad \end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray} \langle \hat{b}^\dagger_n \hat{b}_m\rangle_{ss}&= \sum_{l=1}^N\bigg[\frac{(\gamma_a\bar{\gamma}\!+\!\lambda_{k_l}^2)\gamma_b^h \!+\!\lambda_{k_l}^2\gamma_a^h}{\bar{\gamma}\lambda_{k_l}^2\!+\!\gamma_a\gamma_b\bar{\gamma}} \bigg]\mathcal{S}^{(b)}_{n,k_l}\mathcal{S}^{(b)}_{m,k_l}, \end{eqnarray} \tag{ 31 }$$

$$\begin{eqnarray} \langle \hat{a}^\dagger_n \hat{b}_m\rangle_{ss}&= \sum_{l=1}^N\bigg[\frac{i\lambda_{k_l}(\gamma_a\gamma_b^h-\gamma_b\gamma_a^h)} {\bar{\gamma}\lambda_{k_l}^2+\gamma_a\gamma_b\bar{\gamma}}\bigg]\mathcal{S}^{(a)}_{n,k_l}\mathcal{S}^{(b)}_{m,k_l}, \end{eqnarray} \tag{ 32 }$$

and $\langle \hat{b}^\dagger_n \hat{a}_m\rangle_{ss}$ follows from eq. (32) exchanging the labels $a\leftrightarrow b$. In case that sublattice effective temperatures are the same, $n_{0}=\gamma^h_a/\gamma_a=\gamma^h_b/\gamma_b$, we recover the familiar scenario, induced by the presence of local losses in the Lindblad equation (2) in which the only non-zero moments are $\langle \hat{a}^\dagger_n \hat{a}_n\rangle_{ss}=\langle \hat{b}^\dagger_n \hat{b}_n\rangle_{ss}=n_{0}$ $\forall\, n$. However, in the general case off-diagonal terms are also present, which can lead to the presence of correlations in the stationary state of the system. Indeed, as the master equation (2) is bilinear in the mode operators, if the initial state is Gaussian the stationary state is Gaussian too, and we can fully describe it constructing the covariance matrix, σ, from the above solutions [49], which displays the characteristic block diagonal structure of our system: $\sigma=\bigoplus_{l=1}^N \sigma_{k_l}$. Notice that even though the difference in sublattice effective temperature can induce non-local correlations in the stationary state, such correlations are not related to the transient dynamics and to the presence of EPs and will not be reported here.