Room-temperature lasing from nanophotonic topological cavities

The study of topological phases of light underpins a promising paradigm for engineering disorder-immune compact photonic devices with unusual properties. Combined with an optical gain, topological photonic structures provide a novel platform for micro- and nanoscale lasers, which could benefit from nontrivial band topology and spatially localized gap states. Here, we propose and demonstrate experimentally active nanophotonic topological cavities incorporating III–V semiconductor quantum wells as a gain medium in the structure. We observe room-temperature lasing with a narrow spectrum, high coherence, and threshold behaviour. The emitted beam hosts a singularity encoded by a triade cavity mode that resides in the bandgap of two interfaced valley-Hall periodic photonic lattices with opposite parity breaking. Our findings make a step towards topologically controlled ultrasmall light sources with nontrivial radiation characteristics.

Equations (S2) can be treated as a matrix equation ω[ψ 1 , ψ 2 ] =Ĥ(k)[ψ 1 , ψ 2 ], whereĤ(k) conforms to the two-band model of the general form,Ĥ whereσ = (σ x , σ y , σ z ) is a vector consisting of the three Pauli matrices, and the vector h(k) = (h x (k), h y (k), h z (k)) captures topology and provides a convenient way for its visualization. The vector h(k) determines pseudospin textures for the eigenstates of the model. Specifically for staggered graphene, Expanding h(k) in the vicinity of the Dirac points k = K ± + δk, where δk is a small detuning, we obtain the effective Hamiltonian at two valleyŝ where v D = √ 3 2 ta is the velocity parameter, and µ = ta 2 8 . Keeping only the first order in δk, we get where δU = [M, −M ] is the staggered onsite potential.

B. Analytical solution for edge states at the domain wall
Next, we consider a honeycomb lattice with a staggered on-site potential in the stripe geometry. The stripe consists of two domains with the zigzag interface, where the sign of the parity breaking is flipped, as shown in Fig. S1(a). The tight-binding model for this lattice is written as a system of discrete equations: where k ≡ k x is the momentum vector along the x-aligned interface, ρ = a/2, ψ s,j (n) is the field at site n, j located in the domain s = I, II, index j = A, B denotes sublattice, and integer argument n is a dimer number. At the domain wall, the TBM equations are: The stripe is considered finite, which formally implies the following boundary conditions at the external boundaries of the stripe Here, we focus on the edge states localized at the domain wall. They are concentrated at sites (0,B) in domain I and sites (0,A) in domain II. The solutions for the edge states assume the form: The parameter κ −1 characterizes the decay length away from the interface. Due to parity symmetry of the system, the amplitudes at the interface are connected as where the sign E = ±1 corresponds to symmetric/antisymmetric solutions with respect to the interface. Utilizing the ansatz (S10) with the boundary conditions, we get the continuity condition for the components a I = a II , b I = b II . Substituting Eq. (S10) into Eq. (S8) (for domain I) results in a system of two equations We then apply the condition b I = Ea I ω = −tE − 2tE cos(kρ)e κ and derive the dispersion relation of edge states in the form which can be rewritten as Four valley edge states with different parities are found, among which two are localized at the interface made of a dimer of negative mass −|M | and the other pair is confined to the positive mass |M | dimer. Near the Dirac point at k x = 4π/3a, we can employ the expansion cos(kρ) = cos(k x a/2) ≈ − 1 2 − √ 3 4 δk x a, and for the antisymmetric edge state the spectrum is simplified to We find out that Eq. (S15) for |M | |t| transforms to ω ≈ −v D δk x , while for |t| |M | it is approximated by . The decay rate depends on the propagation constant as follows κ √ 3 2 a = ln M 2 + 4t 2 cos 2 kρ + |M | |2t cos kρ| .
(S16) At |M | |t| this expression turns to κ = |M |/v D , that suggests a constant decay factor for the valley edge states near the Dirac points.
Alternatively, spectrum of edge states can be calculated numerically by solving the eigenvalue problem for a zigzag ribbon composed of two domains. We first compose N × N matrix for one of the domainŝ The second matrixĤ II is obtained by replacement M → (−M ). The final matrixĤ F for a finite stripe of size 2N ×2N is then constructed by stitching the two matrices and incorporating a bond between them in the matrix elements The outer zigzag cut of such finite ribbon supports edge states with a flat band dispersion ω = M . The analytically derived dispersion of the edge modes perfectly agrees with the numerical tight-binding calculations, see Fig. S2. Imposing periodic boundary conditions H F (1, 2N ) = −t, H F (2N, 1) = −t corresponds to the geometry with two domain walls, whose spectrum includes all 4 solutions given by Eq. (S14).

C. Corner states
In the limit of weak coupling |t| → 0, |t| |M |, we can find the frequencies of edge states from two coupled equations for a dimer a II (0) As above, we assume (−M ) interface in the geometry of Fig. S1. Splitting of the level results in frequencies ω 1,2 = −M ± t with symmetric and antisymmetric eigenfunctions (S19) Similarly, in the limit |t| |M |, for the states localized at one corner depicted in Fig. S1(b), we write approximate solutions based on the coupled equations for three coupled sites (a corner trimer) which yields three levels ω 1 = −M , ω 2,3 = −M ± √ 2t, and corresponding eigenvectors: These simple considerations give us basic intuition. With increasing t, the frequency of the corner state gradually moves across the band gap from −M to M , and finally disappears merging into the upper bulk band.
In the schematic Fig. S1, we distinguish two sublattices and dimers positioned at (la 1 + ha 2 ) by two integer Miller indices l, h. Here we choose the basis vectors a 1 = a(1/2, − √ 3/2), a 2 = a(−1/2, − √ 3/2). In the next order in t/|M |, we analytically derive and the corresponding distribution of the excited elements near the corner When t ≈ 2M , the frequency of the lowest-frequency staggered mode V 2 is approaching ω ≈ M . Writing equation for three corner sites (for localized solutions the nearest neighbor a 00 remains negligible, |a 00 | |b 00 |), and taking into account a mirror symmetry of the structure, a −10 = a 0−1 , we deduce a −10 = a 0−1 ≈ −b 00 /2. Remarkably, these relations are found in agreement with the numerically calculated distribution of the out-of-plane magnetic field component H z near the domain-wall corner in the triade cavity mode of our open metasurface. The field is seen primarily concentrated around these three holes. Next, for sake of simplicity, we assume that a reduced scattering potential V (r) associated with the corner mode is dominated by this distribution and can be roughly described as delta functions with a weight V 0 centered at the air holes: where δ n are offsets of elements with respect to the corner hexagon center. The Fourier transformṼ (k) is obtained by integrating the potential Being multiplied by the radiation pattern of the elements (simplistically treated as point magnetic dipoles), it essentially determines the far-field diagram in the reciprocal space. The norm of this transform given by hosts a zero-order singularity and exhibits elliptical deformation, as shown in Fig. S3 calculated for approximate parameters of the structure studied experimentally. Fitting numerically calculated band structures of the photonic metasurface (with varying structural parameters as described in the main text), we estimate effective parameters of the tight-binding model, t/M ∼ 1.8 -2.2, a = 1. Figure S4 shows an exemplary spectrum of a small-scale triangular cavity depicted in Fig. 1(b,c) of the main text. In the tight-binding calculations, the cavity is embedded in a large rectangular domain with a side length of 30 dimers. The band gap contains 8 modes, the triade cavity mode composed of coupled in-phase corner states is marked with a red circle. Some discrepancies between the tight-binding and photonic models can be associated with a finite thickness of the dielectric slab and electromagnetic interactions beyond nearest neighbors. For instance, in the tight-binding model including only nearest-neighbor interaction, the spectrum of bulk dimerized graphene is symmetric about zero energy. The simulated photonic band structure does not exhibit this property. Similarly, incorporating next-nearest-neighbor interactions in the tight-binding model lifts the symmetry of the spectrum. However, as we have shown, the existence of the triade cavity mode can be inferred within the nearest-neighbor tight-binding approximation.
Remarkably, circularly polarized sources operating at frequencies within the minigap can be employed to selectively excite corners of the cavity, see Fig. S5. This directional coupling stems from the chirality of the evanescent edge states [see Figs. S5(a,b)]. In turn, individual corner modes emit beams with singularities in the far field.
For a large triangular cavity where the corners are weakly coupled, the dependence of the triade mode frequency on TBM parameters t/M is shown in Fig. S6. It has two asymptotics: given by expression Eq. (S22) for small t, and for t approaching 2M , where ∆t = 2M − t, and coefficient C ∼ 0.5. The range covered by the edge state dispersion ω = t − M 2 + 4t 2 cos 2 (k x a/2) at fixed t/M is defined by the minimum ω min = t − √ M 2 + 4t 2 and maximum ω max = t − |M | values also plotted in Fig. S6.

II. Continuum model of valley-Hall insulators
By substituting δk x,y = −i∂ x,y , into Eq. (S7), the evolution equations for the two-component wavefunction are transformed to (2 + 1)-dimensional Dirac model For propagating bulk modes e −iωt+iδkxx+iδkyy in a homogeneous medium with a constant mass, the reduced two-band model (S7) gives the (hyperbolic) dispersion relations for two bands: where δk 2 = δk 2 x + δk 2 y .

A. Valley-Hall edge states in the continuum Dirac model
In the framework of the continuum Dirac model, we next derive dispersion of valley-Hall edge states at the domain walls created by mass inversion. We first determine the effective boundary conditions based on the symmetry of equations. If we consider a straight interface along x axis, owing to the translational invariance, the momentum k x is a good quantum number, whereas k y should be replaced by the real-space derivative −i∂ y . Accepting the ansatz Ψ 1,2 = ψ 1,2 (y)e −iωt+iδkxx , equations (S29) are rewritten as (−ω + M (y))ψ 1 (y) + (δk x − ∂ y )ψ 2 (y) = 0 (ω + M (y))ψ 2 (y) + (−δk x − ∂ y )ψ 1 (y) = 0 .
The correspondence between the continuum and and tight-binding models is recovered through

B. Topological cavities
We consider a topological cavity in the geometry, where a circular domain ρ < R with negative mass (−M ) is surrounded by the insulator with the positive mass M > 0.
The recovered coupled equations (S29) can be formulated in the polar coordinate system (ρ, ϕ) for two valley blocks System (S34) possesses solutions with harmonic time dependence and radial symmetry where m = 0, ±1, ±2... is the azimuthal number. The eigenproblem in the polar coordinate system reads The cavity's spectrum yielded by Eqs. (S36) is discrete and contains edge states that occur in pairs with opposite pseudospins. Radial functions u, v for the K + valley obey the coupled ordinary differential equations Plugging u from Eq. (S37b) into Eq. (S37a), we obtain a second-order differential equation for v: Eq. (S38) is simply the Bessel differential equation, whose solution for ω 2 < M 2 are modified Bessel functions (the Bessel functions for complex arguments): v = CI m+1 ( √ M 2 − ω 2 ρ) for the region containing ρ = 0, and v = AK m+1 ( √ M 2 − ω 2 ρ) outside the circle. The ratio of coefficients C/A and frequency ω are found from the continuity boundary condition at ρ = R: (S39) Thus, the dispersion relation is given by the implicit relation For the cavities of radius R much larger than the decay length of the states into the bulk insulator, we can exploit the following asymptotics in Eq. (S40): . With this approximation, we obtain In the limit κR 1, the problem is reduced to the quasi-planar interface considered above by replacing m R with a wavenumber δk x . This suggests that for the modes strongly localized in the vicinity of the boundary ρ = R we can use analytical solutions obtained above for the planar geometry. Note, in the open system the modes propagating along the domain wall are leaky. With the cavity's size and radiative losses increasing, the spectrum tends to the continuum limit.

C. Bound-state formation at the domain-wall corner
We perform the Taylor-series expansion of the edge state dispersion from Eq. (S14) in the vicinity of the wavenumberK: Accepting the monochromatic exp (iωt) process and substituting operators into the dispersion relation (S41), we then recover the equation for the slowly varying field amplitude where v g = dω dk x kx=K is the group velocity, α = − 1 2 d 2 ω dk 2 x kx=K is the second-order dispersion coefficient. For the localised solution at rest v g = 0, hence,Ka = π. In this case α = a 2 t 2 2|M | ,ω = t − |M |, and Eq. (S42) transforms into a parabolic equation, In this setting, the corner of the domain wall plays the role of a trapping potential in Eq. (S43). Next, we denote by s the coordinate along the domain wall path [see Fig. S1(b)] and rewrite Eq. (S43) introducing the potential U (s) and assuming A =Ā(s)e −i∆ωt ∂ 2Ā ∂s 2 − ∆ω αĀ + U (s)Ā = 0 , where U (s) is given by with the potential amplitude A U ≡ A U ∆t 2M being a function of the parameter ∆t 2M ; R is the radius of the path curvature. Here, we specifically consider the regime t → 2M , when the edge state dispersion curve nearly merges to the upper bulk band at k x a = π. At vanishing R → 0 the potential turns to the Dirac delta function U (s) = δ(s) 2 a A U and from Eq. (S44) we find the decay rate κ = ∆ω α = 1 a A U , and the frequency of the trapped state This formula can be used to fit the numerically calculated dispersion of corner modes choosing A U = C 2 ∆t 2M , that coincides with Eq. (S28) .

III. Effective electromagnetic Hamiltonian from the plane wave expansion method
We derive the effective Hamiltonian in the vicinity of Dirac points for a 2D honeycomb photonic crystal slab with air holes by using the plane wave expansion method. Accepting the time dependency ∼ exp(−iωt), we start with Maxwell's equations, where k 0 = ω/c is the wavenumber in free space, ω is the angular frequency, and c is the speed of light. Material response is assumed to be described by effective parameters of a two-dimensional medium, namely, a scalar dielectric permittivity ε(x, y) and the constant magnetic permeability µ = 1. Implying ∂ z = 0, from Eq. (S47) we obtain the governing equation for the H z (x, y) electric field component of TE polarization