Non-Hermitian bath model for arrays of coupled nanoresonators

Nanophotonics systems have recently been studied under the perspective of non-Hermitian physics. Given their potential for wavefront control, nonlinear optics and quantum optics, it is crucial to develop predictive tools to assist their design. We present here a simple model relying on the coupling to an effective bath consisting of a continuum of modes to describe systems of coupled resonators, and test it on dielectric nanocylinder chains accessible to experiments. The effective coupling constants, which depend non-trivially on the distance between resonators, are extracted from numerical simulations in the case of just two coupled elements. The model predicts successfully the dispersive and reactive nature of modes for configurations with multiple resonators, as validated by numerical solutions. It can be applied to larger systems, which are hardly solvable with finite-element approaches.


Non-Hermitian bath model for arrays of coupled nanoresonators: supplementary document 1. MAXWELL'S EQUATIONS RESOLUTION FOR AN ENSEMBLE OF ALUMINUM GAL-LIUM ARSENIDE NANOCYLINDERS
The proof-of-concept system for the non-Hermitian analytical quantum formalism consists of an ensemble of Al 0.18 Ga 0.82 As nanocylinders, which have recently drawn a great deal of attention because of the potential of dielectric nanostuctures on low-index substrate in terms of light control, nonlinear optics and topological physics. They can be seen as nanoresonators, whose properties are individually tailored by their geometry (Fig. S1). We have purposely selected cylindrical geometry to better convey the main message of our study. Eigenmodes of cylindrical structures are close to spherical resonators' ones, which provides a basis of mode close to Mie resonances. For equal diameter and height, in-plane dipoles (MD x and MD y ) and normal dipole (MD z ) are degenerate in energy. For a fixed height, varying the radius of the cylinder shifts the eigenfrequencies of dipolar modes, and lifts their degeneracy. Furthermore, nanocylinders are relatively easy to implement with e-beam lithography and are experimentally very accessible. Consequently, they give a large flexibility and represent a convenient toy system for more advanced characterisations, designs and future experiments. In practice, Al 0.18 Ga 0.82 As nanoantennas are supported by a low optical index semiconductor substrate like Aluminum oxide (n AlOx = 1.6). However, since such substrate only results in second-order correction to the nanocylinders eigenmodes, here we neglect it and suppose the nanostructures to be surrounded by air. Fig. S1. Effect of nanocylinder radius on the real part of the eigenvalues of the magnetic dipolar modes for a 400 nm-high single nanocylinder. Orange curve corresponds to MD z , and blue curve MD x and MD y .
To compute test data for the analytical non-Hermitian formalism (see Figs. 1 and 2 in main text, and Fig. S2), two algorithms were implemented. The former relies on near-field spatial correlations that enable us to track any given mode when the gap d is varied. At each iteration, the electric field of each computed mode from FEM simulation is interpolated inside and in the close vicinity of each nanocylinder, and stored as an array of data. Spatial correlation coefficients are calculated by comparing the interpolated field from two successive iterations. With this method, by selecting the maximum of correlation coefficient one can reconstruct the evolution of a given mode when the gap is varied. Moreover, this algorithm can be initialized with the single element simulation: the modes of an N-element chain are then tracked by groups of N hybrid modes stemming from a single-element mode.
The second algorithm relies on a direct identification of modes based on the evaluation of their complex field components in the whole integration space. For example, MD z has no E z component, and two equal E y and E x components. PML modes are discarded due to their higher field density in the PML layer. Each calculation of this algorithm is thus independent of the others. This procedure can be very versatile, assuming that the user can properly extract sorting condition from physics of Maxwell's equations. It proves more efficient with longer chains, for which eigenmode coupling to PML becomes more prominent in a finite calculation volume.

DERIVATION OF THE QUANTUM LANGEVIN EQUATION
From the Hamiltonian (Eqs. (1-4) of the main text) and the commutation relations, the equation of motion for the continuum bath modes in the Heisenberg picture writes: which can be formally solved aŝ This leads to the dynamical equation forâ j : Γ jk (τ) = Θ(τ) dη g jη g * kη e −iω η τ , where Θ(τ) is the Heaviside step function. Under the limit t 0 → −∞, Eq. (S3) can be rewritten as Finally, the commutator can be calculated as which, together with the nearest-neighbour-coupling assumption, completes the derivation of the quantum Langevin equation (5) of the main text.

DETERMINING THE COUPLING FUNCTIONS FROM THE MAXWELL EQUATION SO-LUTIONS
The eigenvalues of equation (10) in the main text are given by and therefore The simulation data of MD x for the N = 2 case are shown in Fig. S2 (d), with which we will demonstrate the determination of the coupling functions in our analytical model. For a simplified treatment, with no a priori knowledge on the frequency dependence of the reservoir functions, we expand them up to the first order in ω: This assumption makes Eq. (S10) quadratic in ω, and the resonant frequencies ω res ± = ω ± − iγ ± can be expressed analytically in terms of our model parameters: Solving this for each value of the gap d, we obtain a possible set of fit parameters plotted in Fig. S3 as function of the gap, resulting in the same curves as in Fig. S2 (d). While both sets {A 1 , A 2 , B 1 , B 2 , J} and {A 1 , A 2 , −B 1 , −B 2 , −J} give the correct resonant frequencies ω res ± , one of the solutions should be discarded in order to let the mode parity match. In particular, the FEM solutions can be identified as either "bonding" (non-zero electric field at midpoint between nanocylinders) or "antibonding" (zero electric field at midpoint) for each gap value. We then choose the solution set such that the eigenvector associated to an antibonding (resp. bonding) mode is proportional to (-1,1) (resp. (1,1)). This allows us to determine the fitted functions with no sign ambiguity, as shown in Fig. S3.

EFFECTIVENESS OF THE ANALYTICAL MODEL
To provide a measure of the efficacy of the approach, we quantify in Fig. S4 and S5 the relative prediction error of the analytical model without and with the simplified next-nearest neighbour coupling for the cases N = 3 and N = 4, benchmarked against the FEM numerical solutions, which shows the good predictability of the analytical model.