1. Introduction

Most physical systems are subject to inherent losses, either because of dissipation or openness, and are described mathematically by a non-Hermitian Hamiltonian. This generally implies that the eigenstates of the Hamiltonian are complex as well as non-orthogonal. In such systems, the interrelation between pairs of eigenstates under the variation of a set of external parameters is essentially dictated by the existence of exceptional points (EPs). At an EP not only eigenvalues but also eigenstates degenerate, i.e. eigenvectors have the same amplitudes and a constant phase shift. The eigenvalue surfaces display a singular topology in the presence of an EP [1] while the eigenstates accumulate a residual geometric phase upon encircling the EP in parameter space [2,3]. Since their introduction by Kato in 1966 [4], EPs have been associated with numerous often counterintuitive physical effects [5–8] that include mode hybridization [9], quantum phase transitions [10], reversal of the pump dependence of a laser [11], Parity-Time (PT) symmetry breaking [12,13], increased sensitivity to perturbations [14,15] or spectrally broad coherent perfect absorption [16]. They have been observed experimentally in a variety of different systems such as microwave billiards [17], chaotic optical microcavities [18], optomechanical arrangements [19], Fabry-Pérot resonators [20], two level atoms in high-Q cavities [8], coupled photonic-crystals [21], coupled acoustic cavities [22], optical fibres [23] or electronic circuits [24].

Open random media constitute a special class of non-Hermitian systems. Here, modal confinement depends almost entirely on the degree of scattering. In the regime of strong scattering, modes can experience Anderson-localization—a multi-scattering interference phenomenon where the spatial extension of the modes becomes smaller than the system size resulting in transport inhibition [25]. These disorder-induced localized states, for example, provide a natural optical cavity in random lasers [26–28] and showed to be good candidates for cavity quantum electrodynamics (QED) [29,30], with the main advantage of being inherently disorder-robust. By locally tuning the disorder, e.g. by gradually adjusting the dielectric permittivity, these states can be engineered to spatially and spectrally overlap to form so-called necklace states [31–33], which represent open channels in a nominally localized system [34,35]. The formation of necklace states is considered a key mechanism in the transition from localization to the diffusive regime [36]. While EPs and, in general, PT-symmetry have been studied in the context of disordered media [21,28,37–41], coalescing localized modes and their associated degeneracies have not been studied so far.

In this work, the coalescence of two Anderson-localized optical modes in a two-dimensional (2D) dielectric random system is invoked by manipulating the dielectric permittivity at two different locations in the system. We first propose a general analytical analysis to follow the spectral and spatial evolution of modes in 2D dielectric open media. Then, this approach is applied to the specific case of Anderson-localized modes to identify the position of an EP in the parameter space. The prediction is tested numerically with finite element method (FEM) simulations. Remarkably, in this open disordered system a conventional two-mode model, prominently used in a variety of EP-related studies, would fail badly. The modification of a few single scatterers affects a multitude of modes and eventually influences the position of the degeneracies. This emphasizes the complexity of eigenstate interactions in disordered media. We believe that our approach opens the way to a controlled local manipulation of the permittivity and the possibility to engineer hybridized and coupled localized modes. Furthermore, we think this approach can be easily extended to other kinds of resonant systems, e.g. coupled arrays of cavities [42,43].

2. Theory and applications

2.1 Generalized eigenvalue problem

First, we consider the general case of a finite-size dielectric medium in 2D space with an inhomogeneous dielectric constant distribution $\epsilon (\mathbf {r})$. The speed of light is set to $c$ = 1 such that frequencies are hence measured in multiples of $c$ having a unit of $\text {m}^{-1}$. The electric field is assumed transverse and thus its polarization is neglected such that it satisfies the scalar Helmholtz equation

Now, consider two locations $\mathbf {R}_1$ and $\mathbf {R}_2$ where the permittivity is varied

The permittivity distribution $\tilde {\epsilon }(\mathbf {r})$ describes a new disordered system associated with a new set of eigenstates $(\tilde {\Omega }_i, |{\tilde {\Psi }_i}\rangle )_{i\in \mathbb {N}}$. The electric field of the modified system written in terms of the basis of the original random system reads

Projecting $\langle {\Psi _j^{*}}|$ onto Eq. (8) using Eq. (2) and the biorthogonal product in Eq. (4) leads to

If we consider a finite set of $N$ modes, Eq. (9) can be written conveniently in the form of a generalized eigenvalue problem

Equation (11) extends this result to any number of interacting modes $N>2$ essentially describing a network of linearly coupled oscillators.

2.2 Identifying an EP between two Anderson-localized states

The formalism presented above is now applied to a 2D random collection of 896 circular dielectric scatterers (radius $60~\text {nm}$) with dielectric permittivity $\epsilon = 4$ embedded in a host material of index $\epsilon _{\text {mat}} = 1$, with a filling fraction of $40\%$ [Fig. 1(a)]. The system dimensions are $L\times L = 5.3\,\times 5.3 \,\mu \text {m}^{2}$. The interference effects, originating from multiple scattering leads to spatial confinement by the disorder known as Anderson localization and characterized by a localization length $\xi$. In the spectral range considered in the following, $\xi$ is computed through an averaged correlation between modes [49] and estimated to be around $1 \mu \text {m} \ll L$. Thus, the energy of the modes is mostly confined in the system: The modes are localized. The two circular regions with $340~\text {nm}$ diameter centered at $\mathbf {R}_1$ and $\mathbf {R}_2$ are shown in Fig. 1(a) as red and green circles, respectively. Within these regions the dielectric permittivity of the scatterers is varied from $\epsilon$ to $\epsilon + \Delta \epsilon _1$ and $\epsilon + \Delta \epsilon _2$, respectively.

 

Fig. 1. (a) 2D random medium: 896 scatterers with dielectric permittivity $\epsilon = 4$ are embedded in vacuum $\epsilon _{\text {mat}}=1$. The system is open at its boundaries. The permittivity is modified in the two colored regions, red and green circles, centered at positions $\mathbf {R}_1$ and $\mathbf {R}_2$, respectively. (b) Unperturbed eigenvalues $\Omega _i$, $i \in 1,\dots,90$, computed by FEM and sorted in the complex plane according to the distance $d(1,i)$ (see text). The insert points out eigenvalues of interest (namely $\Omega _1$ and $\Omega _2$). (c) and (d) Spatial distribution of the absolute value of the amplitude for the eigenvectors $|{\Psi _1}\rangle$ and $|{\Psi _2}\rangle$, respectively.

Download Full Size | PDF

The initial modes $(\Omega _i,\, |{\Psi _i}\rangle )$, $i \in 1,\dots,N$, which are the only input requested by Eq. (11), are computed using FEM simulation [50,51] with absorbing boundary conditions that are placed $0.4\, {\mu }\text {m}$ away from each side of the system. A large number of modes ($N = 90$) are computed for the initial system [Fig. 1(b)] in a narrow spectral range and we check that none of them are degenerate. Two localized states $|{\Psi _1}\rangle$ and $|{\Psi _2}\rangle$, corresponding to $\Omega _1$ and $\Omega _2$ respectively, are selected for being spectrally close [Fig. 1(b)] but spatially distinct [Fig. 1(c),(d)]. The spectral distance between mode $i$ and mode $1$ is defined as $d(1,i) = |\Omega _1 - \Omega _i|$ and is visualized in the color-coding of Fig. 1(b). Here, mode $2$ has the largest spectral overlap with mode $1$ but we will see that the influence of other nearby modes cannot be neglected in the modal interaction. The initial modes obtained via a FEM are spatially defined in a finite spatial domain $V$ while the biorthogonal product of Eq. (4) is an integral over all $\mathbb {R}^{2}$. However, it can be split into an integral over $V$ and a second term at the boundary of $V$ which replaces outside propagation [52]. Because the modes are localized their amplitudes along the boundary of $V$ are very small. Neglecting the edge term in the biorthogonal product leads to an inaccuracy of around $0.8\%$ in the ensuing calculation of the EP position. Boundary terms can therefore be safely neglected in the biorthogonal product.

Now, the permittivity distribution $\epsilon (\mathbf {r})$ is altered by $\left (\Delta \epsilon _1,\Delta \epsilon _2\right )$ and the corresponding eigensystem is calculated from Eq. (11) for $N=60$ interacting modes for each pair $\left (\Delta \epsilon _1,\Delta \epsilon _2\right )$ in this two-dimensional parameter space, from which a new set of modes $(\tilde {\Omega }_i,\, |{\tilde {\Psi }_i}\rangle )$ is obtained. The spectral distance $\tilde {d}(1,2) = |\tilde {\Omega }_1 - \tilde {\Omega }_2|$ is shown in Fig. 2(a) for all values of the pair $\left (\Delta \epsilon _1,\Delta \epsilon _2\right )$. The sharp drop of $\tilde {d}(1,2)$ to zero observed at $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}} = (0.939,0.90)$ reveals the existence of an EP, where both eigenvalues become degenerate. The evolution in the parameter space of the real and imaginary part of the eigenvalues, $\tilde {\Omega }_1$ and $\tilde {\Omega }_2$, is represented separately in Fig. 2(c),(d). The intricate topology of self-intersecting Riemann sheets at precisely the position $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}}$, corroborates the existence of an EP, where both the real and imaginary part of the eigenvalues cross. At this singular position, the corresponding eigenfunctions become collinear. This is demonstrated in Fig. 2(b), where the function $(1-|\langle {\tilde {\Psi }_1}|{\tilde {\Psi }_2}\rangle |)^{-1}$ is found to diverge.

 

Fig. 2. (a) The eigenvalue difference $\tilde {d}(1,2)_{\text {dB}} = |\tilde {\Omega }_1 - \tilde {\Omega }_2|_{\text {dB}}$ (arbitrary units) vanishes at the position $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}}$ in parameter space. (b) The function $(1-|\langle {\tilde {\Psi }_1}|{\tilde {\Psi }_2}\rangle |)^{-1}$ measures the collinearity of the eigenvectors and has a clear divergent maximum at $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}}$ where the two modes coalesce. (c) and (d) Real and imaginary part of the eigenvalues $\tilde {\Omega }_{1}$ and $\tilde {\Omega }_{2}$. The eigenvalue surfaces display the typical structure of intersecting Riemann sheets with a singular point at $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}}$. (e) Trajectories $\text {Im}(\tilde {\Omega }_1)= \text {Im}(\tilde {\Omega }_2)$ and $\text {Re}(\tilde {\Omega }_1)= \text {Re}(\tilde {\Omega }_2)$ marked by red and blue dots, respectively: Both curves join continuously at the EP. The insets show the real and imaginary part of modes 1 and 2 near the EP: the modes become collinear and satisfy $\tilde {\Psi }_1 = -i \tilde {\Psi }_2$ ($\text {Re}(\tilde {\Psi }_1)=\text {Im}(\tilde {\Psi }_2$) and $\text {Im}(\tilde {\Psi }_1)=-\text {Re}(\tilde {\Psi }_2)$).

Download Full Size | PDF

2.3 Testing the theoretical prediction

To substantiate the analytical prediction, the first two eigenvectors ($|{\tilde {\Psi }_1}\rangle$, $|{\tilde {\Psi }_2}\rangle$) of the modified permittivity distribution $\tilde {\epsilon }(\mathbf {r})$ are calculated via numerical FEM simulations for each point $(\Delta \epsilon _1,\Delta \epsilon _2)$ on a grid in a range that encloses $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}}$. The sampling is set to $0.04$ along the two parameters $\Delta \epsilon _1$ and $\Delta \epsilon _2$. The numerically computed eigenvalues merge at $(\Delta \epsilon _1,\Delta \epsilon _2) = (0.92,0.88) \pm (0.04,0.04)$, such that the predicted value from Eq. (11) is within the error bars. The collinearity of the eigenvectors is also numerically confirmed (not shown). Finally, the coalescing mechanism of modes 1 and 2 is described in Fig. 2(e) where modes with identical real (imaginary) parts are plotted in the parameter space, $(\Delta \epsilon _1,\Delta \epsilon _2)$, as blue (red) dots. At the EP where real and imaginary parts are equal simultaneously, the red and blue trajectories meet. Near this position, we compute the real (imaginary) part of the two eigenvectors, $\tilde {\Psi }_1(x,y)$ and $\tilde {\Psi }_2(x,y)$. Their distribution is shown in the insets of Fig. 2(e)]. Starting from the unperturbed modes $\Psi _1(x,y)$ and $\Psi _2(x,y)$ [shown in Fig. 1(c),(d)], both eigenvectors are progressively modified when approaching the EP. In the vicinity of the EP, $\text {Re}(\tilde {\Psi }_1)$ converges to $\text {Im}(\tilde {\Psi }_2)$ and $\text {Im}(\tilde {\Psi }_1)$ to $-\text {Re}(\tilde {\Psi }_2)$. Therefore, the eigenvectors display the same amplitude and a phase shift of $\pi /2$, $\tilde {\Psi }_1 = -i \tilde {\Psi }_2$, i.e. they become collinear at the EP. The amplitude of the degenerate eigenstate at the EP, $|{\tilde {\Psi }_{\text {EP}(1,2)}}\rangle$, is shown in the inset of Fig. 3(a). It forms a beaded chain which connects both ends of the system and is similar to necklace states studied in [31].

 

Fig. 3. (a) Norm of the biorthogonal projection of the degenerate eigenstate at the EP (coalescence of modes 1 and 2) with the first 45 initial modes $| \langle {\Psi _{i}^{*}\,}|{\tilde {\Psi }_{\text {EP}(1,2)}}\rangle |$. The absolute value of the amplitude of $|{\tilde {\Psi }_{\text {EP}(1,2)}}\rangle$ is shown in the inset. The four panels on top display the spatial distributions of initial modes 1, 2, 4 and 5, respectively. Their respective biorthogonal projection coefficient is numbered. (b) Prediction of EP position in the parameter space [$\Delta \epsilon _{1,\text {EP}}$ (red line) and $\Delta \epsilon _{2,\text {EP}}$ (blue line)] for a number of modes $N$, ranging from 2 to 60. The prediction converges for $N \ge 55$. The two error bars point out the position of the EP obtained by FEM computation.

Download Full Size | PDF

2.4 Influence of other localized modes

Here, in contrast to higher-order EPs [53], the EP under study characterizes the coalescence of only two modes. Yet, when the permittivity distribution gets modified to converge to this EP, the others modes play a role in the deformation of the two coalescing modes. To investigate how much different localized modes of the system contribute to the degenerate state $|{\tilde {\Psi }_{\text {EP}(1,2)}}\rangle$, we measure their spatial overlap with the degenerate eigenstate. The norm of their biorthogonal projection, $| \langle {\Psi _{i}^{*}\, }|{\tilde {\Psi }_{\text {EP}(1,2)}}\rangle |$, is shown in Fig. 3(a). Remarkably enough, it demonstrates the fading, though not negligible, influence of nearby modes with $i>2$. Both modes 1 and 2 contribute to $60\%$ to the degenerate EP state, while other modes have a vanishing contribution with increasing index, though modes 4 and 5, still contribute significantly for $25\%$. Their influence is also highlighted in Fig. 3(b), where the position of the EP, $(\Delta \epsilon _1,\Delta \epsilon _2)_{\text {EP}}$, is calculated using the analytical prediction, Eq. (11), for an increasing number of contributing modes $N$. Fluctuations of the EP position are large for small $N$, while the convergence is reached for value of $N$ larger than $55$. This observation emphasizes a fundamental point, namely that the spatial amplitude distribution of the degenerate state at the EP results from the coupling of a large number of modes.

3. Conclusion

In conclusion, we proposed a general analytical theory to study the evolution of modes in open media when the permittivity is varied. This non-perturbative approach relies on the linearity of $\epsilon (\mathbf {r})$ in the Helmholtz equation: For any variation $(\Delta \epsilon _1, \, \Delta \epsilon _2)$ a new set of modes can be computed from Eq. (11) simply from the knowledge of the initial system. This approach describes the system by an infinite set of modes acting as oscillators coupled via the modification of the permittivity. We considered the specific case of random media in the localization regime and show that our theory can be used to investigate the mode coupling and hybridization resulting from a local perturbation. In particular, by changing the local index at two different small areas, two modes are brought to coalescence and give rise to an EP. Remarkably, the accuracy of the theoretical prediction is shown to strongly depend on the number of modes considered. A large number of initial modes is required, which shows that the mode evolution is dictated by multimode interactions. In terms of experimental realization, such a manipulation of the disorder can be easily implemented on existing setups [27,54,55]. The permittivity landscape can be shaped reversibly, e.g. using laser illumination to switch on nonlinear index changes [56–58] or by moving one or several scatterers [59]. EPs can be calculated for any pair of modes and even generalized to hybridization of three or more eigenstates, which opens the way to the control of light-matter interaction in random media. The effective creation of necklace states, for instance, allows the formation of open transmission channels in opaque systems. Alternatively, using mode repulsion in the vicinity of an EP, the disorder can be engineered to increase the spatial confinement of the modes and consequently their Q-factor. Finally, we believe this approach can be extended to others types of complex optical systems, e.g. photonic crystal cavity arrays used for quantum simulation [42,43,60] where fabrication inaccuracy could be compensated a posteriori by such an external control.