Complex-frequency dispersion analysis

The dispersion and loss curves for the transverse magnetic (transverse electric)-polarized modes of the planar waveguide were calculated semi-analytically using a modified transfer matrix method35, which allows for extraction of the complex-ω modes that remain valid solutions at the SL points. Let us consider a planar metal-dielectric stack structure of N layers of thickness d i (i=1…N) that are stacked up in y direction and extend infinitely in the x–z plane. For isotropic media, the field component ? z , orthogonal to the propagation direction x, can be written as

where m i is the layer permittivity (permeability) for ? z =H z (? z =E z ). Momentum conservation dictates that the y component of the wavevector depends on the x component k (also called propagation constant). The boundary conditions between individual layers translate into

where the matrices M i are given by

The bound mode profile leads to the initial condition in the superstrate, while in the substrate the condition can be expressed as

where define the decay constants in the cladding layers.

For plasmonic and active layer structures, the characteristic function F(ω, k2) is a complex-valued function. Its roots lie on a hyperplane and can be found efficiently using a combination of the argument principle and the Newton–Raphson methods.

In practice, when solving the characteristic equation one considers one of the arguments (ω or k2) as a real parameter and the other as complex eigenvalue. For real ω, the complex-k eigenvalue describes a mode with spatially decaying (or increasing) amplitude, a solution that resembles those obtained by cw mode injection into the waveguide. These solutions do not exist at the SL points, as the corresponding mode cannot propagate. Instead, one needs to consider complex-ω modes (with a fixed, real k). As these modes decay in time (rather than space), they do not show the typical back bending and mode splitting of the complex-k dispersion curves at the SL points.

The roots of the dispersion function can be classified according to their behaviour in the cladding layers. The bound modes, characterized by an exponential decay in the super- and substrate (Re[γ 1 ], Re[γ N ] > 0), are associated with a discrete spectrum of eigenvalues. In contrast, the radiative continuum is described by modes with a purely oscillatory behaviour (Re[γ 1 ],Re[γ N ]=0)36. The third class of modes, so-called leaky modes, show a divergent behaviour in the cladding. These modes, despite not being part of the complete basis set, are observable as resonances in the reflection spectrum. For the (weakly leaky) structure under investigation, leaky modes accurately capture the changes to the modal dispersion and the additional radiative loss caused by the coupling of the waveguide mode to the continuum of free space modes.

Maxwell-Bloch Langevin time-domain simulations

The time-resolved simulations were performed using the finite-difference time-domain method in combination with auxiliary differential equations that describe the material dispersion of the metal cladding and the nonlinear, spatially resolved polarization response of the gain media. Using the Maxwell-Bloch formalism we modelled the gain medium as a four-level system with pump (0→3) and emission (2→1) transitions37 that we associate with polarization densities P a (r, t) and P e (r, t). The dynamics of the polarization densities can be mapped on a set of differential equations

which describe the dynamics of light–matter interaction via the dipole coupling of the electric field to the electronic transitions. The resonance frequencies ω i (and transition frequencies ), the spectral half-widths γ i , the phenomenological coupling constants (that are connected to the emission and absorption cross-sections σ 0,i ), as well as the real part of the host refractive index are determined by the choice of gain material.

The temporal evolution of the occupation densities N 0 to N 3 , which in turn determine the inversions ΔN e (r,t)=N 2 (r,t)−N 1 (r,t) and for emission and absorption, is given by the set of equations

The lifetimes τ jk account for non-radiative relaxation processes between the levels. When considering a spatially homogeneous and constant local pump rate , the term connected to the polarization density of the absorption transition, (ℏ ω r,a )−1(∂P a /∂t+Γ a P a )E, can be replaced by .

Dissipation-induced fluctuations of the polarizations and carrier densities are included in the Maxwell-Bloch equations via Langevin noise terms37. The noise terms are given by F e (r, t), F a (r, t) and F 00 (r, t) to F 33 (r, t), where N cell is the number of four-level systems per unit volume cell of the finite-difference time-domain grid. Below the lasing threshold the noise acts as source for amplified spontaneous emission giving rise to the recorded background spectra (Fig. 3b). Above the lasing threshold, the noise acts as a seed aiding the transient build-up of the coherent lasing fields through feedback and stimulated emission processes in the gain medium.

Structural and material parameters

The planar metal-dielectric stack structures comprises an optically thick metal substrate, a thin dielectric core layer and a second (thin) metal cladding layer. Further control over the modal dispersion is achieved by introducing a dielectric layer on top of the stack structure when necessary26. In our analysis, we assume a constant permittivity of ε s =11.68 for the dielectric layers, a value that is typical for III–V semiconductors (such as InGaAsP). The dispersive permittivity of the metal layers, in this case a transparent conductive oxide, is described using a Drude model. In accordance with experimental data38, we use ω p =3.13·1015 rads−1 for the plasma frequency of indium tin oxide. The collision rate γ is assumed to be reduced by a factor of 10 to 1.07·1013rads−1 as could be achieved by fabricating films of very high quality and cooling them to low temperatures5. In comparison to metals such as silver or gold, the plasma frequency for, e.g., indium tin oxide can be controlled through doping and can be shifted into the visible wavelength range to around 600 nm making this material a suitable option for SL lasing in the near-infrared.

Genetic optimization of broadband dispersionless modes

The principle of sub-wavelength SL lasing relies on the formation of a flat band in the modal dispersion that maximizes the range of wavevectors contained within the frequency bandwidth of the gain. The flatness of the band can be achieved by introducing two or more SL points, which enforce a monotonous variation of the dispersion in between them. We employ a genetic algorithm technique to isolate structures of a given design that exhibit two SL points and the smallest possible frequency variation in a selected mode.

Every planar metal-dielectric stack structure is uniquely characterized by the number of layers, the layer thicknesses and their permittivities, forming an ordered parameter set, the so-called chromosome of the genetic algorithm. On each iteration of the algorithm, we hold a population of N chromosomes, which are evaluated with respect to a fitness function that returns a single numerical value representing suitability of the structure. The population is sorted to this order and a small percentage of the fittest chromosomes are passed to the next generation unchanged. The remaining places in the new population are filled by ‘tournament selection’, where each tournament is a subset of the current population whose members are chosen at random. The two fittest chromosomes from each tournament are then mutated by randomly perturbing parameters of the chromosome and optionally crossed over by swapping sections of the parameter list between the two tournament winners. The two new chromosomes are added to the next population, and the process of tournament selection is repeated until the population set is complete. The progress of each iteration can be monitored by considering the fittest chromosome in the population and the algorithm is terminated once this chromosome surpasses a threshold fitness.

For the fitness function, we first evaluate the dispersion relation of a given chromosome (that is, structure) and perform a root find for two SL points. If the calculation fails during this step, the function returns to zero, otherwise we determine the average band velocity and return its reciprocal value

The structure with the highest value of f is deemed fittest. It is common to initialize the population with chromosomes representing structures already known to exhibit two SL points, which may in turn be found by optimization using a different, appropriate fitness function.