Spatial and Spectral Mode-Selection Effects in Topological Lasers with Frequency-Dependent Gain

Frequency-Dependent Gain Matteo Seclì,1, a) Tomoki Ozawa,2 Massimo Capone,1, 3 and Iacopo Carusotto4 International School for Advanced Studies (SISSA), Via Bonomea 265, I-34136 Trieste, Italy Advanced Institute for Materials Research, Tohoku University, Sendai 980-8577, Japan CNR-IOM Democritos, Via Bonomea 265, I-34136 Trieste, Italy INO-CNR BEC Center and Dipartimento di Fisica, Università di Trento, I-38123 Povo, Italy


I. INTRODUCTION
Topological lasers (in short, topolasers) are one of the most promising applications of topological photonics. Such devices are obtained by including a suitable gain material in a topological system so to induce laser oscillation in a topological edge state. 1,2 Stimulated by pioneering theoretical proposals, 3-7 experimental realizations were first reported for the zero-dimensional edge states of one-dimensional arrays. [8][9][10][11] Extension to nanolasers based on corner states of two-dimensional lattices was reported in Ref. 12.
Scaling up in dimension, the crucial advantages for optoelectronic applications offered by topological lasing into the one-dimensional edge modes of a two-dimensional lattice have been theoretically highlighted: [3][4][5]13 the topological protection of chirally propagating one-dimensional edge modes appears as a promising strategy towards an efficient phase-locking of the laser oscillation at the different sites. In this way, gain can be distributed over a large number of sites, while maintaining a globally stable single-mode coherent emission, which is very promising to realize high-power coherent sources. Experiments along these lines were reported shortly afterwards using photonic crystals under a strong magnetic field 14 and arrays of coupled ring microcavities, 15 followed by more recent valley-Hall quantum cascade 16 and telecomwavelength 17 lasers.
These experimental advances have stimulated an active theoretical research to characterize the peculiar properties of the novel devices. Whereas the experiments in Refs. 14 and 15 have shown a clean single-mode emission from topolasers, the possibility of secondary instabilities as a result of the interplay of optical nonlinearities and slow carrier dynamics has been theoretically pointed out in Ref. 18. A semiclassical study of the novel features a) matteo.secli@sissa.it introduced by the chirality of the lasing state was reported in Refs. 19 and 20. Finally, extensive theoretical studies based on a stochastic approach have anticipated the robustness of the long-time coherence against static disorder by including quantum and thermal fluctuations into an idealized model of topolaser. 13 In this work, we take a different point of view and we investigate the various mode-selection mechanisms that determine whether a topolaser device is going to lase in an edge or in a bulk state. Rather than dealing with the complex nonlinear dynamics of the lasing state as done in Refs. 13,18, and 20, we focus our attention on identifying the mode that is responsible for the first instability of the vacuum state. This is a common strategy in laser physics 21 and typically provides a good intuition on the system behaviour. For instance, if different modes of the laser resonator have different spatial profiles, a specific mode can be selected just by increasing its spatial overlap with the gain material. Analogously, if modes are well separated in frequency, it is enough to tune a narrowband gain material in the spectral vicinity of the desired mode to suppress competing modes. Here, we take inspiration from the recent topolaser experiment in Ref. 14 and we investigate more subtle schemes where a suitable combination of spectral and spatial mode-selection mechanisms conspires to stabilize laser oscillation into a chiral edge state.
The structure of the work is the following. In Sec. II we introduce the physical system under investigation, namely a photonic Harper-Hofstadter lattice embedding population-inverted two-level atoms as gain medium, and we develop the theoretical model based on a Bloch-Harper-Hofstadter set of equations. In Sec. III, we show how the use of a narrow-band gain stabilizes the edge mode lasing even when the gain material is uniformly distributed across the whole system. In Sec. IV we show how a suitable combination of spectral and spatial selection mechanisms is able to stabilize the edge mode lasing under weak conditions on the gain lineshape and its spatial localization. Some experimental implications of our results are discussed in Sec. V. Conclusions are finally drawn in Sec. VI. Additional details on the derivation of the theoretical models, on the topological lasing features, and on our spatial-spectral mode selection mechanism are given in the supplementary material.

II. THE THEORETICAL FRAMEWORK
In this Section we summarize the theoretical model used for our calculations. As an archetypal model, we consider a photonic Harper-Hofstadter lattice where optical gain is introduced by including population-inverted two-level atoms at each site. For a complete derivation of the equations of this Bloch-Harper-Hofstadter model, we refer the interested reader to Sec. S.1 of the supplementary material.

A. The Bloch-Harper-Hofstadter model
We consider a two-dimensional Harper-Hofstadter lattice 1,22 where neighboring sites are connected through a hopping Hamiltonian with a non-trivial hopping phase. In the Landau gauge, this can be written as where a † m,n (a m,n ) is the operator that creates (annihilates) a photon at the (m, n) site. All sites are assumed to have a bare photon frequency ω cav , the real and positive parameter J quantifies the hopping strength, and ϑ = 1/4 is the flux per plaquette of the synthetic gauge field that is responsible for the topological properties.
As it is sketched in Fig. 1, each site is modeled as a photonic resonator of frequency ω cav and decay rate Γ, and in which a photon is created (annihilated) by a † (a). Each resonator is provided with a frequency-dependent gain medium, which is modeled here as a collection of N population-inverted two-level atoms (TLAs). At each site, the dynamics of the atoms forming the gain medium is then described by the following Hamiltonian: where ω eg = ω e − ω g is the energy difference between the atomic levels, the light-atom coupling g is assumed to be equal for all j = 1, . . . , N atoms and, for each atom, σ + = |e g| and σ − = (σ + ) † = |g e| are the raising and lowering operators between the ground |g and excited |e states. Analogously, the atomic population difference on each atom is quantified by the σ z = |e e| − |g g| operator. Each TLA is incoherently pumped from the ground to the excited state at a pumping rate γ g , while the reverse spontaneous decay from |e to |g occurs at a rate γ e .
Under a mean-field approximation, we replace the photon field operators a m,n with their classical C-number expectation values α m,n = a m,n and we assume the atomic density matrix to have a factorized form. The one-atom Hamiltonian terms H (j) at,1 ≡ H at,1 and the oneatom density matrices ρ at are then equal for all atoms at a given site and their dynamics is captured by a Lindblad master equation of the form where the first term gives the coherent evolution induced by the atom-field dynamics in (2) and L e = |g e| = σ − and L g = |e g| = σ + are the jump operators for the decay and pumping processes.
Projecting (2) and (3) onto the atomic ground and excited states then recovers the Bloch equations of the semiclassical theory of lasers. 23,24 Together with the field dynamics determined by the hopping Hamiltonian (1), these equations constitute the full set of equations of our Bloch-Harper-Hofstadter model. Measuring all energies and times in units of J and J −1 respectively, and setting the frequency zero at the empty cavity frequency ω cav , these equations have the form:                     ρ m,n ee = γ g ρ m,n gg − γ e ρ m,n ee + i α m,n ρ m,n ge − α * m,n ρ m,n eg ρ m,n gg = γ e ρ m,n ee − γ g ρ m,n gg − i α m,n ρ m,n ge − α * m,n ρ m,n eg ρ m,n eg = −i ω eg − iγ ρ m,n eg − iα m,n ρ m,n ee − ρ m,n gg α m,n = − Γα m,n + iGρ m,n eg + i α m+1,n + α m−1,n + e −i2πϑm α m,n+1 + e +i2πϑm α m,n−1 (4) where ρ m,n ee (ρ m,n gg = 1 − ρ m,n ee ) is the average atomic population of the excited (ground) state of the atoms located at site (m, n) and ρ m,n eg is the corresponding coherence. In the following, we will assume for simplicity that no additional decay channel is acting on the atomic coherence ρ eg in addition to the unavoidable ones coming from pumping and decay, γ = γ eg = (γ g + γ e )/2.
The efficiency of the gain process enters via the G g 2 N coupling strength, proportional to the number of atoms N per site and to the square of the elementary light-atom coupling g. Indicating with V the effective field volume at each site, the light-atom coupling scales as usual as V −1/2 (see Sec. S.1 of the supplementary material), which makes G proportional to the atomic density N/V . If a real atomic gas is used as the gain material, the gain strength can be tuned by changing the density of the gas. In a solid state photonic crystal where the The energy difference between the two atomic levels is ωeg; pumping of the atoms from |g to |e occurs at a rate γg, while their spontaneous decay from |e to |g occurs at a rate γe. All cavities decay at an equal loss rate Γ. The synthetic magnetic field is included as a non-trivial hopping phase: a photon that hops around a plaquette (orange arrows) picks up an extra phase 2πϑ due the synthetic magnetic flux, in our case ϑ = 1/4.
TLAs are used to model a more complex electromic dynamics, the same effect can be achieved by varying the filling factor within the unit cell and/or the overlap of the Bloch mode with the gain material. To keep our analysis as simple as possible, in the following we will use G as the parameter controlling the strength of the gain in the different regions of space. This choice automatically includes the possibility of having different atomic densities in the different regions; compared to the pumping rate γ g , it also allows to simplify the presentation by avoiding the complications due to the simultaneous dependence of several other parameters on γ g .

B. Frequency-dependent gain
The peculiar features of the frequency-dependent gain can be simply understood by looking at the laser operation in a single site geometry. Indicating with ω L the lasing frequency, the steady-state is determined by the single-site limit of the Bloch equations (4). This gives α(t) =α e −iω L t with the amplitudeα satisfying the gain/loss balance equation with the effective pump strength and the saturation coefficient Equation (5) is formally analogous to the one of a broadband saturable gain considered in Refs. 5 and 20, with the key difference that the parameters P, β are here frequency-dependent. In particular the effective pump strength P involves a Lorentzian factor γ/[(ω eg − ω L ) 2 + γ 2 ] accounting for the non-trivial gain spectrum: this is centered at ω eg and has a HWHM set by the atomic decoherence rate γ.
This frequency-dependent gain directly reflects into an analogous dependence of the laser threshold. In a singlesite geometry this is immediately obtained from Eqs. (5) and (6) as the lowest value of G for which (unsaturated) gain exceed losses P ≥ Γ. This leads to the threshold condition where is the single-cavity lasing threshold exactly on resonance, that is for ω L = ω eg . As expected, the threshold is minimum when laser operation occurs on resonance with the pupolation-inverted atoms. Then, it increases quadratically with the detuning ω L −ω eg : the faster the atomic decoherence rate γ, the weaker this increase. In the following of the work, we will exploit this frequency-dependence of the threshold as a way to select the desired mode for lasing.

III. NARROWBAND GAIN
In the previous Section we have restricted our attention to the single-site case. This provides us the conceptual building blocks to understand laser operation in a topological lattice. As a first step, in this Section we consider the simplest case where the narrowband gain spectrum is concentrated within a topological gap. Under this conditions, a stable topological lasing can be obtained even under a spatially uniform pumping. Still, because of mode-pulling effects, a non-trivial relation is found between the lasing frequency and the bare frequencies of the discrete set of edge modes.
This narrowband configuration can be obtained for instance by tuning the atomic frequency ω eg in the middle of the topological bandgap with a gain linewidth γ much smaller than the gap width, as sketched in Fig. 2(c). In this way, the frequency-dependence of gain strongly increases the effective threshold for laser operation in the off-resonant bulk band states while the one for edge state lasing remains almost unaffected.
Laser operation in this regime is illustrated in Fig. 2(ab). Emission into the edge state is stable and monochromatic and remains so up to high pump strengths well above the laser threshold. 25 Only at very high pump strength, also the bulk modes go above threshold and the dynamics recovers the chaotic behaviour found in Ref. 20 for a broadband gain distributed in the whole system.
In contrast to the broadband gain case where the laser frequency ω L is typically locked to the bare mode frequency ω 0 , for a narrowband gain a sizable mode pulling effect can occur on the laser frequency. 24 The lasing frequency results then from a weighted average of the atomic resonance ω eg and the bare mode frequency ω 0 via the mode pulling formula where S = Γ/γ is the so-called the stabilization factor. For equal Γ = γ, the stabilization factor is S = 1 and the mode pulling effect becomes a simple average. Let us explore the impact of this effect in our case of topological laser operation for a narrowband gain centered inside the topological bandgap, thus perfectly overlapping with the edge state dispersion. The numerical calculations illustrated in Fig. 2(d) show an interesting interplay of mode pulling with the intrinsic discreteness of the edge mode frequencies discussed in Ref. 20 and due to the intrinsic periodic boundary conditions around the system. The steady-state lasing frequency is plotted as a function of the atomic resonance position ω eg as blue dots: as expected, the overall behavior roughly follows ω eg (solid red line), but also displays a series of slant steps.
When mode pulling effects are negligible, for instance in a broadband gain case, one expects the system to lase at the frequency of the edge mode that is closest to the resonance ω eg for which gain is strongest. In our plot, this would correspond to a staircase of flat steps separated by a spacing ∆ω = 2πv g /L determined by the overall length L of the system edge and the group velocity v g of the edge mode. In the presence of the narrowband gain, laser operation still occurs in the discrete mode that is closest to ω eg . This gives again a staircase structure, but mode pulling effects make the steps to have now a finite slope equal to S 1+S . In the figure we have taken γ = Γ, so the stabilization factor is S = 1 and the slope is 1/2 (black dashed lines).

IV. BROADBAND GAIN
In the previous Section, we have seen an efficient scheme to stabilize topolaser operation with a uniformly distributed gain, by spectrally concentrating the gain spectrum in the topological gap. While conceptually interesting, this scheme is hardly useful in practical semiconductor systems, where the gain linewidth is typically comparable if not larger than the width of the topological band gaps. In this Section we will explore a more sophisticated scheme that is able to stabilize topological lasing in a much wider range of parameters of potential technological relevance. As a practically relevant scenario, we consider a broadband gain with a linewidth much larger than the topological gap.
The configuration we consider is inspired by the experiment Ref. 14 and is sketched in Fig. 3(d): we consider a central region, which has a narrow topological gap, surrounded by a region with a much wider and topologically trivial gap. Chiral boundary modes are localized at the interface between the two regions. We assume that the gain is stronger in the trivial region than in the topological region, and also assume that the width of the gain is much larger than the narrow topological gap but comparable to the large trivial gap. Even though the distribution of the gain material is homogeneous in each respective region, we can expect lasing in the topological boundary states because the boundary states at the interface penetrates also into the trivial region with stronger gain. In this Section, we explain in detail how this idea works.
We first explain how we prepare a narrow topological gap in the central region. We start from the ϑ = 1/4 Harper-Hofstadter model, which contains multiple topological band gaps with topological invariants adding up to zero. We want to isolate one topological band gap from these multiple gaps. To this end, we add a checkerboardshaped on-site frequency detuning ±∆: the frequencies of the (m, n) sites are thus alternated and equal to ω cav + ∆ · (−1) m+n . The photonic bands of such a bipartite Harper-Hofstadter model are shown in panel (a) of Fig. 3: because of the checkerboard detuning, the Dirac touching point between the middle two bands opens into a trivial gap with size ∼ 2∆. The two small topological gaps between the lower two bands and the upper two bands persist together with the corresponding edge states. We focus on the topological band gap of the two upper bands; the gain spectrum is centered around the frequency of the two upper bands. The two lower bands are, instead, off resonant, and would not be relevant in the laser operation and the discussion below.
Next, we explain how we prepare the surrounding region with a wide trivial gap. We again start from the ϑ = 1/4 Harper-Hofstadter model and add a checkerboardshaped detuning ∆ trivial , which is larger than ±∆ in the topological region. We add a global shift of all site frequencies by ω trivial so that the large topologically trivial gap between the middle two bands is centered around the two upper bands of the topological region. The corresponding photonic bands are shown in panel (b). Although the two upper and two lower bands have narrow topological gaps, they are pushed away by the large ∆ trivial , and thus we can focus on the effect of the wide trivial gap between the middle two bands. We call this surrounding region a "trivial" region in this sense. The gain spectrum is indicated by the yellow shading in panel (b), which is centered at the middle of the wide trivial gap and completely encompasses the topological gap in the central region.
The gain strength in the surrounding trivial region, G trivial , can be reinforced either by locally increasing the pumping strength as done in Ref. 15 or by keeping a spatially uniform pumping but increasing the density of gain material with respect to the central region, as discussed in Sec. II A. Focusing on this latter case, we can write G trivial = G·d, where d can be interpreted as the effective density of gain material in the surrounding region relative to the central region, and treat G/G res,0 as a global measure of the uniformly distributed pumping strength in units of the single-resonator resonant threshold.
The results of the numerical simulations are summarized in Fig. 3, where we show a phase diagram of the different regimes of laser operation as a function of the relative effective density of gain material in the surrounding region d = G trivial /G and of the pumping strength in the central topological region in units of the resonant, single-site threshold, G/G res,0 .
When the surrounding trivial region is purely passive  and does not display any gain (d = 0, Fig. 3(i) and (j)), the system is almost equivalent to a bipartite 15 × 15 lattice without the surrounding region. We therefore expect the system to only lase above the resonant single-site lasing threshold, G/G res,0 = 1, as shown by the red region at the bottom of the phase diagram. Since the gain is effectively broadband with respect to the upper pair of photonic bands in the central region, both bulk and boundary modes equally participate in the lasing process, forbidding a stable topological laser operation (panel (j)). Compared to the bulk states, boundary states are even slightly disfavoured by the worse spectral overlap with the gain spectrum and by their evanescent tail that penetrates into the surrounding trivial region and reduces the spatial overlap with the gain region.
We can induce lasing from the topological edge modes at the boundary by making the gain in the surrounding trivial region to be stronger than the one in the central topological region, that is d > 1. In this case, a region appears in the parameter space where the system displays the desired topological laser behaviour (panel (g)). Thanks to their evanescent tail overlapping with the stronger amplifying surrounding region, the effective threshold of the topological boundary modes is in fact pushed well below the one G/G res,0 = 1 of the bulk modes (thick solid black line in the phase diagram), opening a window where only these modes can lase (blue region). For G/G res,0 > 1, the bulk modes also go above threshold and the laser emission acquires a complex multi-mode character (blue-to-white-to-red region). Still, thanks to the stronger gain of the surrounding region, the intensity of boundary mode lasing can remain significantly stronger than the one of the central bulk modes even at values of G above the single-cavity lasing threshold (panel (h)).
For even higher values of the surrounding density d above the dash-dotted black line in the phase diagram, we reach a point where the spectral selection is no longer sufficient to suppress bulk lasing in the surrounding region and topological lasing is no longer possible. In this phase (yellow area in the phase diagram), the much stronger gain of the surrounding region makes the laser emission to be concentrated in this region (panels (e) and (f)).
A quantitative analytical discussion of the location of the transition lines in the phase diagram is given in Sec. S.4 of the supplementary material. As expected, the area of the topological lasing region in parameter space can be increased by either increasing the trivial bandgap in the surrounding region or by using a narrower gain spectrum. This trend is confirmed by the additional numerical simulations with different values of the parameters that are shown in Sec. S.3 of the supplementary material.

V. DISCUSSION
In the previous Sections we have concentrated our attention on a Harper-Hofstadter model which provides a relatively straightforward insight into the basic effects, but our conclusions extend to any combination of lattices with suitable spectral and topological properties. In particular, we expect that our physical conclusions extend even outside the tight-binding approximation that has been made in all theoretical studies of topological lasing so far.
As a most intriguing example, the results of our calculations are compatible with some key observations of the pioneering experiment in Ref. 14 that, to the best of our knowledge, remain so far unexplained. In particular, topological lasing was observed in this experiment without the need to concentrate gain along the edge separating the topological and trivial regions as it was instead the case in other experiments. 5,16 A key difference between the devices used in these works consists in that the topological system used in Ref. 5 is surrounded by empty space, while in Ref. 14 the central topological system is surrounded by a topologically trivial region where the field can penetrate with a significant evanescent tail. Most importantly for our purposes, the outer region displays a larger filling factor of the unit cell (compare Figs. 2A and 2B in Ref. 14). For an equal level of optical pumping, we can thus reasonably expect the gain to be stronger in the outer region, which corresponds to d > 1 in our model. As a result, the overlap of the edge state with this stronger amplifying region favours topological lasing with respect to bulk lasing in the central region. At the same time, the much wider extension of the trivial photonic band gap of the outer region forbids laser operation in the outer region thanks to the natural frequency-dependence of gain in the used semiconductor quantum well material.
While these arguments provide a suggestive interpretation of experimental observations, they are of course not yet completely sufficient to rule out other possible explanations. For instance, in analogy to the arguments put forward in Ref. 17 for a different geometry, another potentially relevant mechanism for stabilizing the edge mode lasing could originate from the weaker losses of the edge mode compared to the ones of bulk modes. 26 In the specific system of Ref. 14, reduced radiative losses may in fact originate from the evanescent tail in the outer trivial region where bulk modes in the vicinity of the trivial gap are below the light cone. In our model theory, the reduced radiative losses of the trivial region could be explicitly included via a reduced Γ of the outer sites, but we expect their effect to be similar to the one of the increased gain G trivial considered in our calculations. On this basis we are confident that the qualitative conclusions of our theory directly apply to the experiment. However, a firm and definitive unraveling of these questions requires accurate experimental measurements and comprehensive microscopic calculations of the band structure and of the radiative and non-radiative decay rates of the different modes, 26 which go beyond our work.
As a final point, it is worth briefly mentioning some straightforward experiments that may serve to shed light on the possible interpretations of the experimental observations even in the absence of a direct measurement of the Q factor of the different modes. In the IR spectral region of the experiment, magnetic effects are quite weak as signalled by the smallness of the topological gap. This implies that the magnetic field is crucial to induce the topological edge state, but has a minor effect on the bulk regions. As a result, according to our theory in the absence of any magnetic field no laser operation should be observed up to powers well above the topological laser threshold. Some evidence in this direction is found by comparing Figs. 3(b) and 3(c) of Ref. 14. Further experimental insight could be obtained by keeping the magnetic field on and ramping up the pump intensity well above the topological laser threshold. According to our model, as discussed in Sec. IV, going up in gain strength G should move the system from the blue topological lasing region into the red/yellow ones of multimode bulk lasing. In particular, we expect that the threshold for bulk lasing at high gain strengths should be almost insensitive to the applied magnetic field.

VI. CONCLUSIONS
In this work, we have developed a general semiclassical theory of topological laser operation that is able to include the peculiar structure of the photonic modes of the underlying topological lattice and the frequencydependence of a realistic gain material. As a specific example of application of our theory, we have investigated the lasing threshold in a configuration that displays a subtle interplay between the spatial overlap of the modes with the gain medium and the position and width of the frequency gaps in the different regions. Based on our theory, an interpretation of the recent experimental observations of topolaser operation under a spatially uniform pumping 14 is proposed.
A natural next step will be to include our theory of frequency-dependent gain into the Bogoliubov description of collective excitations around a topologically lasing state 27 so to characterize the stability of realistic models of topological laser operation in the different regimes of gain parameters. This will be of great interest as a new tool to tame all those instability mechanisms that may hinder a clean and efficient single-mode topological laser emission in practical semiconductor devices.

SUPPLEMENTARY MATERIAL
See supplementary material for full derivations of the theoretical models, additional simulations at different gain widths and an extended discussion on the topological lasing features.

ACKNOWLEDGMENTS
We warmly acknowledge continuous stimulating exchanges with Boubacar Kanté. We are grateful to Aurelian Loirette-Pelous and Ivan Amelio for continuous discussions on topolaser dynamics and we thank Maxim Gorlach for discussions on our preliminary results. We acknowledge financial support from the European Union H2020-FETFLAG-2018-2020 project "Pho-QuS" (n.820392) and from the Provincia Autonoma di Trento. TO acknowledges support from JSPS KAKENHI Grant Number JP20H01845, JST PRESTO Grant Number JPMJPR19L2, and JST CREST Grant Number JP-MJCR19T1.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

. Model and Equations
Lamb's semiclassical theory of lasing is usually derived by coupling a semiclassical Hamilonian for the light-atom interaction, where the electric field is classical, with Maxwell's equations. The same theory can also be obtained under a mean-field treatment of the photon field, by starting from a fully quantum Hamiltonian, e.g. the Jaynes-Cummings Hamiltonian; while in the former case one has to make the assumption of slowly-varying fields, in the quantum treatment we work under the dipole and the rotating wave approximations.
The starting point is an Hamiltonian description of N identical two-level atoms (TLAs) inside a cavity; the choice of working with TLAs comes from the fact that they're the most simple -yet effective -model of an emitter. The multi-atom version of the Jaynes-Cummings Hamiltonian, also called Travis-Cummings or Dicke model, reads as: 1 where (r) stands for "radiation" and (ar) for "atom-radiation" interaction. Here ω cav is the natural frequency of the cavity, ω eg = ω e − ω g is the energy difference between the atomic levels, g = g eg = g ge = − degE0 is the light-atom ge is the matrix element e | e r | g of the dipole moment -usually called the transition dipole moment, E 0 = ωcav 2ε0V and V is the volume of the cavity. In the Hamiltonian, a † is an operator that creates a cavity photon and, for every atom n, σ + = |e g| and σ − = (σ + ) † = |g e| are operators that describe transitions between the atomic levels, and σ z = |e e| − |g g|. Additionally |e e| + |g g| = I, with I the identity operator.
The central assumption we now make is for the total density matrix ρ to be factorizable as a tensor product between a density matrix for the atoms and a density matrix for the photons, i.e. ρ = ρ N −at ⊗ ρ ph , where the atomic part can be obtained by tracing out the photonic part: ρ N −at = Tr ph ρ. With this consideration one can similarly obtain a time-evolution equation for the atomic density matrix ρ N −at by starting from the evolution of the full density matrix, dρ dt = − i [H, ρ], and tracing out the photonic part. The result is an equation very similar to the one for the full a) matteo.secli@sissa.it density matrix, but in which H is replaced by a mean-field Hamiltonian: where in H MF all the photon operators have been replaced by their expectation values: If we further assume that the atomic density matrix can be factorized as a tensor product of the density matrices of all the atoms, i.e.
where ρ (S5) If we trace away all the other atoms except for the n-th one, we get which is exactly the same time-evolution equation one would get from a single-atom Jaynes-Cummings Hamiltonian. This means that all the atoms are equivalent and have exactly the same evolution, so we can drop the (n) superscript. Finally, we can provide every atom with a decay + pumping mechanism at the master equation level, by adding Lindblad terms of the following form: where the dots . . . indicate the rhs of (S6), L † e = |e g| = σ + (the population of |e is increased) and L e = |g e| = σ − (the population of |e is decreased), and conversely L † g = |g e| = σ − and L g = |e g| = σ + . With these ingredients we can write down the equations for the matrix elements of the atomic density matrix of a generic atom in the cavity: where γ = γ eg +γ ph , γ eg = γe+γg 2 and γ ph , in the more general treatment (see Ref. 2), comes from dephasing processes; for the sake of our discussion though, γ ph = 0. The fact that the evolution of the population of one of the levels is the negative of the evolution of the other one is due to the fact that the Lindblad losses preserve the trace of the density matrix (which is 1), and indeed from ρ gg + ρ ee = 1 we getρ gg = −ρ ee .
We still miss an equation for the cavity field a, or better, for its expectation value. The Heisenberg equation for the field operator is By noticing that σ − is an atom-only operator, so that we can then take the expectation value of (S9) and use the fact the the atoms are all equivalent to write Finally, by redefining a → a for simplicity and by adding a phenomenological cavity loss coefficient Γ, we have a set of coupled semi-classical equations for the elements of the atomic density matrix and the expectation value of the photon field: Interestingly the number of atoms in the cavity, N , appears only in the equation for the average photon field, in front of the off-diagonal atomic matrix element. Indeed, if you compare with the semiclassical theory, 2,3 you find that the time-evolution for the macroscopic electric field iṡ where E + is the positive-frequency component of the electric field and P + is the positive-frequency component of the (macroscopic) polarization. Each atom here has a tiny dipole moment that can be shown to be e r = d eg (ρ eg +ρ ge ) = d eg (ρ eg +ρ * eg ), so one can assume that the macroscopic polarization vector is just the sum of these tiny dipole moments: When you then identify P + ≡ N d eg ρ eg , you geṫ 2ε0 , which has the same form of (S11).

S.1.2. Stationary State
We can get additional insights and build up some intuition by solving the equations (S12) for the stationary state.
In order to do that, let us first perform a change of variable. We define the overall oscillation frequency of the field as ω L -which is not a priori guaranteed to be equal to ω cav , so that it is convenient to define a new variableã as where nowã(t) is a slowly-varying quantity, i.e. |ȧ/ã| ω L . We also defineρ eg = ρ eg e +iω L t andρ αα ≡ ρ αα (α = e, g). If then one solves for the population difference at the steady state, plugs it back into the equation for the coherencẽ ρ eg and finally combines with the equation for the field, the latter takes the following form: where (S20) Here {ω} indicates the dependence on the frequencies (ω cav , ω eg and ω L ) and {γ} indicates the dependence on the microscopic loss coefficients (γ e , γ g , γ ph ).ρ The key point here is that the P coefficient, which if the zero-field population difference is positive then acts as a gain, is frequency-dependent; you get most out of your cavity system when the lasing frequency ω L is in resonance with the transition frequency between the two atomic levels ω eg . If instead your lasing frequency gets farther apart from the transition frequency, then the gain decays by following a Lorentzian behavior -with the linewidth of the Lorentzian controlled by γ. Additionally, the lasing frequency experiences an intensity-dependent frequency shift that decays as the square of the detuning for large detunings but is exactly zero at resonance; this shift is responsible for the so-called mode pulling.

S.1.3. Active Harper-Hofstadter Lattice
Multiple active cavities like the one studied in Sec. S.1 can be further arranged into a lattice with a suitably engineered hopping Hamiltonian. In our case, we are interested in a Harper-Hofstadter model with a hopping term of the form: where m = 1, . . . , N x and n = 1, . . . , N y are the site indices in the x and y direction of the lattice plane. The introduction of a hopping Hamiltonian amounts to an additional term in the evolution of the fields (third equation in (S12)), i.e. the photon field at each site, before taking the expectation value, now evolves aṡ where the dots . . . indicate the rhs of (S9). It is also useful to perform a change of variable so to measure all frequencies with respect to the cavity frequency ω cav . Note that this is different from the transformation done in (S16); while the latter is useful for a theoretical discussion, is of little help in the numerical simulations since the lasing frequency ω L is unknown. In practice, we rewrite a(t) in the following way: where nowα(t) is not a slowly varying quantity, since it oscillates at the detuning between the field and the cavity. We further defineρ eg = ρ eg e iωcavt ,ρ nn ≡ ρ nn (n = e, g),ω eg = ω eg − ω cav andω L = ω L − ω cav , with the consideration that equation (S22) holds for the transformed fieldα as well.
For numerical purposes, it is also useful to perform a rescaling to dimensionless quantities. This is done by performing the re-definitions This parameter, which is proportional to the volumetric density N V of TLAs and to the natural frequency of the cavity, acts as a gain parameter.

(S25)
In all the simulations, we've set the on-site losses to a reasonable value Γ = 0.1.

S.1.4. Lasing with Frequency-Dependent Gain
We now seek the lasing condition of a single lattice cavity, bearing in mind that as we have shown in Ref. 4 the lasing condition for a lattice of identical cavities has to be greater or equal than the threshold of a single cavity.
In actual devices the amplification is controlled by γ g , the rate at which the atoms are excited from their ground state |g to a higher-energy state |e . However, this has also the additional side effect of possibly modifying both the height and the linewidth of the Lorentzian profile of the gain (see Sec. S.1.2), depending on the relative value of γ ph , as well as contributing to frequency-dependent shifts of the lasing frequency itself and of the saturation factor of the amplification term. In order to skim these side effects off the model and obtain a simpler, yet still meaningful description, we'll instead use G as the primary amplification parameter, and we'll write the lasing condition for G itself. Increasing G instead of ρ eg (through γ g ) has the same effect on the amplification term in equation forã, but helps avoiding the effects described above and, if the range of probed G values is not too large, it still gives a faithful description.
The lasing condition in terms of G is easily obtained from the steady-state analysis done in Sec. S.1.2, which shows we must have P 1+β|ã| 2 > Γ. By expanding this expression, we have the threshold condition where ∆ det = ω eg − ω L measures the detuning of the lasing frequency ω L from the atomic transition frequency ω eg . This relation shows that the threshold for a given mode has a quadratic dependence on its frequency, and that the mode that lases first -i.e. for which the threshold is as small as possible -is the one at resonance, ω L = ω eg . In this resonant case, the threshold condition becomes where we have defined G res,0 as the single-cavity resonant lasing threshold.
In the numerical simulations, for the atomic parameters we have set γ ph = 0 and we have taken γ g = 4 3 γ and γ g = 2 3 γ, with 2γ = γ g + γ e being the FWHM of the Lorentzian profile of the gain. With this choice of parameters, we then have G res,0 = 3γΓ. The remaining parameters are changed throughout the text.

S.2. MODE PULLING AND FREQUENCY DISCRETIZATION
It is a general feature of laser physics that a narrowband gain causes an additional mode pulling effect. For a single cavity this means that the lasing frequency does not correspond to the bare cavity frequency, as it is pulled towards the center of the Lorentzian profile of the gain, 2 i.e. towards ω eg : where S = Γ γ is called the stabilization factor. For equal Γ = γ, one has S = 1 and the mode pulling effect becomes a simple average.
In the case of a lattice of cavities, we have to replace ω cav with a frequency ω 0 that corresponds to the specific mode of the underlying lattice that the system is selecting for lasing. This effect has to be considered for a reliable calculation of frequency-dependent lasing thresholds and to correctly relate the lasing frequencies that we measure through a power spectral density analysis to the underlying lattice eigenmodes.
Let us explore the impact of this effect on topolaser operation for ω eg located in the topological bandgap. In this case, one may naively expect that mode pulling effects are irrelevant. The system is in fact expected to lase in the topological edge mode resonant with the gain for which the threshold is lowest. Since it is resonant, this mode does not experience any mode pulling -setting ω 0 = ω eg yields in fact ω L = ω0+Sωeg 1+S = ω eg . In contrast to this expectation, the numerical experiment shown in Fig. S1 shows that we actually have a richer physics due to the intrinsic discreteness of the modes.
In the left panel, reproduced here from Fig. 2(d) of the main text, we have simulated a lattice with narrowband gain for multiple values of ω eg and we have measured the steady-state frequency at which the system was emitting via a power spectral density analysis. While the measured frequencies obey an overall approximate relation ω L ∼ ω eg , the reality is that ω L is bound to assume only discretized values due to the periodic boundary conditions, as discussed in Ref. 4. In a broadband gain regime where mode pulling effects are negligible (S 1), this would result in ω L following a stair-like behavior as a function of ω eg with uniform spacing ∆ω = v g · 2π L determined by the length L of the topological edge. Within each step, the emission frequency is locked to the one of the mode that is closest to resonance and for which gain is strongest. For a narrowband gain, the behavior is instead that of a slant stair and the slope of each of these slant steps is exactly S 1+S (dashed lines). This results from the fact that ω 0 cannot continuously follow ω eg , but it has to do so in steps, so it is actually constant in a certain frequency interval; then ω L = ω0+Sωeg 1+S yields a behavior of type y = a · x + b with a ≡ S 1+S and b ≡ ω0 1+S . In the particular case S = 1, we get a = 1/2 and b = ω 0 /2.
This also gives us the chance to quantitatively verify the relation between the spacing of the lasing frequencies and the group velocity of the topological edge mode in the presence of mode pulling (Fig. S1, right panel). We can calculate the spacing between successive modes by just calculating the intersection frequencies ω int between the fits of type y = 1 2 · x + b (dashed lines) and the line y = x (red line) and then taking the difference ∆ω inc between successive frequencies. We assign each value ∆ω int to the frequency at the center of its corresponding frequency interval; this first order approximation becomes exact for L → ∞. From our periodic boundary condition argument, we should also have that the spacing between successive modes is ∆ω ∼ v g · 2π L ; as it's not clear how to treat the corner sites in the calculation of L, we show values ranging from L = 4 × 14 to L = 4 × 16 as a blue shaded area in the plot, with the solid blue line marking the value L = 4 × 15.
The data shows a good agreement with respect to our predictions, highlighting the importance of including modepulling effects as well as the presence of mode discretization even in the narrowband case. We show here additional simulations of the bipartite Harper-Hofstadter laser with the central topological region and the surrounding trivial region discussed in Sec. IV of the main text. In particular, we have pointed out in the main text that changing the gain linewidth reflects in a change of the area of the topological lasing region in the phase diagram. In Fig. S2 we show simulations with an increased (top row) or reduced (bottom row) gain linewidth when compared to the value used in the main text.
For the increased (reduced) gain linewidth we take 2γ = 8.0 (2γ = 4.0), to be compared with the value 2γ = 5.2 used in the main text. With these settings the gain linewidth is around ×12.5 (×6.3) wider than the linewidth of the topological bandgap of the central region, so it is roughly 96 (48) times wider with respect to the topological bandgap than the gain linewidth shown in the narrowband case. The gain linewidth is also around 57% (29%) the linewidth of the trivial bandgap in the surrounding region.
When using this increased/reduced gain linewidth, the area of the region of parameters that yields a topological laser shrinks/widens; in particular, the upper transition line (which describes the lasing threshold of the surrounding region) is extremely sensitive to the linewidth of the gain. It is possible to predict in advance the effects of such a reduced/increased gain linewidth on the phase diagram by a direct calculation of the transition lines -see Sec. S.4. It is also instructive to explore the phase diagram, with respect to the emitted intensity, along both a vertical and a horizontal cut; we do that in Fig. S3 for the bottom right phase diagram of Fig. S2.
If we fix the effective density of gain material d of the surrounding region at some value d = 9.80 and we plot the emitted intensity vs G/G res,0 (left panel of Fig. S3) we notice that the data can be grouped in four sets, each one fitting a distinct linear branch. The fitted linear branches highlight the presence of three thresholds: the first one for the topological edge, the second one for the central bulk and the third one for the surrounding region.
We can also explore the emitted intensity along a vertical cut of the phase diagram, i.e. as a function of the density d of the surrounding region for a fixed value of the gain strength G/G res,0 = 1.06; for the sake of clarity, it is best to separate the emission from the topological edge region (blue), from the central bulk (red) and from the surrounding region (yellow). First of all we have to make a crucial observation: the scale of the emission from the surrounding region (right vertical axis) is two orders of magnitude bigger than the scale of the topological edge and of the central bulk (left vertical axis). This means that if we look at real-space snapshots of the system when the surrounding region is lasing (e.g. Fig. 3(e)-(f) in the main text) we get the wrong impression that the rest of the system (i.e. the central bulk and topological edge) is not lasing; the rest of the system is indeed lasing, but the value of the emitted intensity is negligibly small when compared to the one emitted from the surrounding region.
The other crucial observation is that the plot shows a clear evidence that the topological edge is indeed lasing by exploiting the unused amplification in the surrounding region. First of all, after the surrounding region is above threshold its emitted intensity has a non-linear ramp-up as a function of the density. This is due to the fact that, as one increases the density of the surrounding region, more and more modes from the surrounding region reach the lasing threshold. The first ones to lase are localized in the bulk of the surrounding region; however, at higher densities, the modes from the surrounding region localized on the border reach the lasing threshold as well. These modes, due to their sizable shared overlap, are in direct competition with the topological edge mode that is lasing in the central region; thus the emitted intensity of the latter starts decreasing (marked by a vertical dashed line in the plot).

S.4. TRANSITION LINES IN THE PHASE DIAGRAM
As pointed out in Sec. IV of the main text, the topological edge mode is able to lase thanks to an amplification leak from the surrounding region. This leak can be quantified and used to analytically determine the transition lines of the phase diagrams in Fig. S2 and in Fig. 3(c) of the main text without having to resort at all to dynamical simulations.
The first ingredient is the knowledge of the energy spectrum and of the spatial structure of the eigenmodes of the lattice Hamiltonian, which we call H bare : where E i is the energy corresponding to the mode with normalized wavefunction ψ i , and ψ i can be thought as a matrix in position space, i.e. ψ i = (ψ i;m,n ), with m and n respectively the x-and the y-index. At this point we can define some spatial masks M topo , M bulk and M trivial ; these masks are matrices in position space and contain 1 in the entries representing points belonging to, respectively, the topological edge, the bulk of the central region and the surrounding region, and zero elsewhere. These masks are the same ones used to calculate the emitted intensity overlaps in the phase diagram, and also allow us to calculate the overlap of a given wavefunction with these three spatial region of the system, e.g. to calculate the overlap S topo,i of the wavefunction ψ i with the topological edge: For a given wavefunction ψ i the overlaps with the three regions sum to 1 since the wavefunction is normalized, i.e. S topo,i + S bulk,i + S trivial,i = 1. The overlaps of each mode can be visualized as in Fig. S4. The overlaps are useful to determine what is the effective amplification G eff,i experienced by a given mode ψ i as compared to the global parameter G characterizing the gain strength. For a completely flat gain profile, if all the system regions have the same density of the bulk, then G eff,i /G = 1 because the overlaps of any given wavefunction with the three regions of the system sum to 1. However, if we change the density d of the surrounding region, this ratio will be in general G eff,i (d)/G = f i (d) = S topo,i + S bulk,i + d · S trivial,i . If d > 1, then f i (d) > 1 and it acts as an enhancement factor for the amplification of that mode.
If the gain is not flat but has e.g. a Lorentzian profile, as in our case, the expression of f i (d) in terms of wavefunction overlaps has to be modified. Namely, we have to add an extra Lorentzian factor that is equal to 1 at zero energy detuning and falls to 0 at large energy differences. Analogously to Sec. S.2, the energiesẼ i are here corrected with respect to E i in order to take into account the mode pulling, i.e.
Here, the mode pulling correction has the effect of increasing the value of L i as compared to a Lorentzian factor calculated with the non-pulled energies E i . So, the overall enhancement factor will be f i (d) = S topo,i + S bulk,i + d · S trivial,i 1 + With these ingredients, we are ready to calculate the transition line between the not-lasing region and the different lasing regions, starting with the one corresponding to a topological laser. The first mode to lase, i.e. the mode that determines the threshold, will be the mode with largest effective amplification, i.e. with greater f i (d). That mode will have an enhancement factor equal tof Since we are looking for the lasing threshold among the topological edge modes, we take the max over the modes that have an overlap S topo,i with the topological edge greater than a certain threshold ε topo , e.g. ε topo = 0.80. Similarly, to find the second transition line at which the surrounding region starts to lase, one defines Note that the mode that maximizes the enhancement factor depends on the specific d being considered. A more crude approximation not involving a maximization over the modes, but that still provides very good results, would be to define e.g.f topo (d) = f itopo (d), where i topo is the index of the topological edge mode that lases more frequently just above threshold in the time-domain simulations; i topo is then independent on the choice of d. A similar consideration holds forf trivial (d) as well.
In general, the threshold for lasing into a given mode is reduced by the same factor f i characterizing the enhanced effective amplification as compared to the d = 1 case for which the threshold is at G/G res,0 = 1. As a result, the transition lines we are looking for are simply described by the equations d = 1 f topo (d) and d = 1 f trivial (d) .
Both transition lines carry a dependency on γ, mediated by the Lorentzian factor in f i (d). However, the transition line between the not-lasing and the topo-lasing regimes has a much weaker dependence. The reason is that, by construction, the topological edge modes are much closer to the center of the Lorentzian than the bands in the surrounding region. Since this detuning contributes quadratically to the position of any transition line, the transition line involving modes farther from the gain center is therefore much more sensitive to γ itself. As a final note, this calculation also allows us to easily determine the lasing thresholds when we vary the thickness of the surrounding region. The system considered in Fig. 3 of the main text has a 5-sites-thick surrounding region; in the left panel of Fig. S5 we show again its phase diagram but we calculate the transition lines of an otherwise identical system which differs by having a 1-site-thick surrounding region. There is a discrepancy between the two, albeit extremely small, which shows that a bare minimum 1-site-thick surrounding region is already enough to open a topological lasing region in the phase diagram. Already with a 2-sites-thick surrounding region even this small discrepancy disappears, and the transition lines becomes virtually indistinguishable from the 5-sites-thick case.
The fact that a 1-site-thick surrounding region is already sufficient comes from the fact that the topological edge modes we want to enhance have an overlap with the surrounding region that is mainly concentrated in the surrounding sites closest to the edge itself, i.e. in the system region covered by a 1-site-thick surrounding region. As can be inferred by comparing the overlap fractions shown in Fig. S4 and in the right panel of Fig. S5, the positive-energy topological edge modes have an overlap with the surrounding region that's at best around 2.3% both when the surrounding region is 1-site-thick and when it's 5-sites-thick, thus contributing around the same enhancing factor in the calculation of the thresholds.