<figure><img src="image_1.png" /></figure>

HAL
open science

Impact of Cavity Geometry on Microlaser Dynamics

Kyungduk Kim, Stefan Bittner, Yuhao Jin, Yongquan Zeng, Qi Jie Wang, Hui
Cao

► To cite this version:

Kyungduk Kim, Stefan Bittner, Yuhao Jin, Yongquan Zeng, Qi Jie Wang, et al.. Impact of Cavity Geometry on Microlaser Dynamics. Physical Review Letters, 2023, 131 (15), pp.153801. 10.1103/Phys-RevLett.131.153801. hal-04239711

HAL Id: hal-04239711

https://hal.science/hal-04239711v1

Submitted on 12 Oct 2023

**HAL** is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Impact of cavity geometry on microlaser dynamics

Kyungduk Kim*,¹ Stefan Bittner*,²,³ Yuhao Jin,⁴ Yongquan Zeng,⁴ Qi Jie Wang,⁴ and Hui Cao†¹

¹Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USA

²Université de Lorraine, CentraleSupélec, LMOPS, 2 rue Edouard Belin, Metz 57070, France

³Chair in Photonics, CentraleSupélec, LMOPS, 2 rue Edouard Belin, Metz 57070, France

⁴Center for OptoElectronics and Biophotonics, School of Electrical and Electronic Engineering, School of Physical and Mathematical Science, and Photonics Institute, Nanyang Technological University, 639798, Singapore

We experimentally investigate spatio-temporal lasing dynamics in semiconductor microcavities with various geometries, featuring integrable or chaotic ray dynamics. The classical ray dynamics directly impacts the lasing dynamics, which is primarily determined by the local directionality of long-lived ray trajectories. The directionality of optical propagation dictates the characteristic length scales of intensity variations, which play a pivotal role in nonlinear light-matter interactions. While wavelength-scale intensity variations tend to stabilize lasing dynamics, modulation on much longer scales causes spatial filamentation and irregular pulsation. Our results will pave the way to control the lasing dynamics by engineering the cavity geometry and ray dynamical properties.

Controlling nonlinear dynamics of complex systems is crucial in, e.g., nonlinear optics, hydrodynamics, and laser physics [1, 2]. It remains, however, a challenge to control semiconductor laser dynamics due to extremely fast inherent time scales [3]. Instead of relying on external feedback [4–8] or optical injection [9–12], we propose a more direct and compact approach based on modifying the intrinsic light-matter interaction inside the laser by tailoring the cavity geometry [13–15].

Broad-area Fabry-Perot cavities, commonly used for semiconductor lasers, often result in lasing instabilities [16–22]. An optical lensing effect in high-intensity regions, caused by spatial hole-burning and carrier-induced refractive index changes, leads to self-focusing of light and the formation of filaments, which are unstable and induce irregular pulsations. We recently showed that modifying the resonator shape can suppress spatio-temporal instabilities, e.g., in a D-shaped microcavity that features fully chaotic ray dynamics [13], or in a stable cavity with optimized mirror curvature that exhibits integrable ray dynamics [14, 15]. The stability of lasing dynamics depends critically on the characteristic length scale of optical intensity variations in the cavity. If high-intensity regions, formed by constructive interference of propagating waves, are too small to induce a lensing effect, filamentation and lasing instability will be prevented. However, it is not clear how the cavity geometry affects the characteristic size of high-intensity regions, and whether it is possible to predict the lasing dynamics based on classical ray dynamics.

To address these questions, we experimentally study how the classical ray dynamics is related to the spatio-temporal dynamics of semiconductor microcavity lasers. The ray dynamics of two-dimensional (2D) optical microcavities is entirely determined by the cavity shape and the boundary conditions. Hence, microlasers based on dielectric resonators correspond to billiards with (partial) ray escape according to Fresnel’s laws [23, 24]. Previously the ray-wave correspondence was studied to reveal the relation between long-lived trajectories and resonances

with high quality (Q) factors [24–28]. Most studies have concentrated on static properties like emission spectra and far-field distributions [29–34], except a few numerical [35, 36] and experimental studies [13, 37–40] of the dynamic properties of deformed microcavity lasers.

Here we extend the paradigm of ray-wave correspondence to engineering the nonlinear lasing dynamics by utilizing classical ray dynamics. We judiciously choose several cavity shapes not only with distinct types of ray dynamics: chaotic vs. integrable, but more importantly, with varying structural sizes and degrees of spatial localization of their lasing modes. We find the strength and occurrence of irregular pulsations of the laser emission are strongly correlated with the local structure size of the lasing modes, which are normally the most long-lived passive cavity resonances. This characteristic size is determined by the local directionality of optical propagation inside the cavity, which can be predicted by the classical ray dynamics.

Figure 1 shows five cavity geometries, D-cavity, stadium, Limaçon, ellipse, and square. These five shapes include cavities with chaotic (D-cavity, stadium, Limaçon) and integrable (ellipse, square) ray dynamics. Furthermore, D-cavity, stadium, and square have spatially extended modes featuring small structure size due to low local directionality of wave propagation, while Limaçon and ellipse have spatially localized whispering gallery modes (WGMs) exhibiting large structure size due to directional wave propagation.

D-cavity and stadium have fully chaotic ray dynamics [41, 42], and long-lived rays explore the bulk of the cavities [Figs. 1(a₁, a₂)]. Due to ray-wave correspondence [23, 43–46], the high-Q modes are spatially extended and irregularly structured [Figs. 1(b₁, b₂)]. The Limaçon-shaped cavity features predominantly chaotic ray dynamics [47–49], but has long-lived trajectories concentrated at the cavity boundary [Fig. 1(a₃)], resulting in irregular WGMs [Fig. 1(b₃)]. The ellipse features integrable ray dynamics [50] with long-lived trajectories along the cavity boundary that are confined by to-
<figure><img src="image_3.png" /><figcaption>FIG. 1. Five cavity geometries with different classical ray dynamics. (a) Typical long-lived ray trajectories in D-shaped ($a_1$), stadium ($a_2$), Limaçon ($a_3$), ellipse ($a_4$), and square ($a_5$) dielectric cavities. (b) Numerically calculated mode profile with high Q-factor in a cavity with smooth boundary. The scale bars are 5 µm long.</figcaption></figure>

tal internal reflection [Fig. 1($a_4$)], resulting in regularly-structured WGMs [Fig. 1($b_4$)]. The square also features integrable ray dynamics, but its long-lived trajectories explore the whole cavity [Fig. 1($a_5$)], and hence its high-Q resonances [51–53] are spatially extended with regular, fine features [Fig. 1($b_5$)].

We fabricate edge-emitting semiconductor microlasers with these five resonator shapes [54]. The devices are fabricated on a commercial laser diode wafer with a GaAs/AlGaAs quantum well (Q-Photonics QEWLD-808) by photolithography and inductively coupled plasma etching, followed by deposition of a top metal contact for electric current injection. Multiple cavities for each of the five geometries, with identical cavity areas ($2.53 \times 10^4$ µm²), are fabricated on the same wafer.

We study their spatio-temporal lasing dynamics experimentally using a streak camera (Hamamatsu C5680/M5676) to record the time-resolved near-field emission intensity profiles. We measure 10 ns long time windows with a resolution of ~30 ps. Lasing occurs in all cavities with electrical pumping at room temperature [54]. Figure 2(a) shows exemplary streak images. The spatial profiles of emission intensity agree with the output patterns of high-Q resonances [40, 55]. The irregular pulsations with typical periods of sub-nanoseconds are caused by unstable lasing dynamics.

We analyze the laser intensity fluctuations using the short-term radio-frequency (RF) spectra and their temporal evolution. We calculate the temporal Fourier transform of every 10-ns-long streak image and average its magnitude squared in space. Then a broadband continuous signal from the spatio-temporal beating of lasing modes and amplified spontaneous emission is subtracted from the RF spectra, which highlights discrete peaks from irregular pulsations due to lasing instabilities [54]. The subtracted RF-spectra in Figs. 2($b_1-b_5$) clearly show different degrees of stability for the five geometries.

For a more quantitative statistical analysis, we characterize the RF-spectra with two measures in Fig. 3(a). The first one (left axis) is the total RF power $S_{tot}$, obtained by integrating the short-term RF spectra [Fig. 2(b)] in both frequency and time. The second one (right axis) measures the frequency of occurrence of RF peaks by the participation ratio of the RF spectra [54]. To account for the cavity-to-cavity variations, we average these quantities over five different lasers per geometry. Both the overall fluctuation power and the occurrence of intensity pulsations vary by several orders of magnitude for different cavity geometries. The D-cavity microlasers have the weakest and rarest pulsations, followed by the stadia. Limaçon cavities have stronger and more frequent pulsations than D-cavities and stadia. The ellipse lasers have the strongest and most frequent pulsations. The square cavities, in contrast, are much more stable than ellipses and close to stadia.

The spatial structure of the lasing modes strongly influences the nonlinear interaction between the optical field and gain material, which in turn affects the strength and occurrence of irregular pulsations. To reveal the underlying mechanism, we numerically characterize the fine structure of the lasing modes. We calculate the passive resonances of 2D cavities with boundary roughness to account for fabrication defects [54]. The cavity dimensions are 10 times smaller than the actual ones to reduce computational load. Furthermore, we employ steady-state ab-initio lasing theory with single-pole approximation (SPA-SALT) to determine which cavity resonances will lase and to calculate their lasing intensities [54, 56–59].

The mode profiles in Fig. 4(b) show significant differences in their fine structure for the different geometries. The D-cavities and stadia [Figs. 4($b_1,b_2$)] show random and isotropic intensity variations on the scale of the intermediate wavelength. For the Limaçon [Fig. 4($b_3$)], in contrast, the fine structure is anisotropic with elongated high-intensity grains. This anisotropy is even more pronounced for the ellipses [Fig. 4($b_4$)]. In contrast, the squares [Fig. 4($b_5$)] feature a more regular structure with a feature size similar to that of D-cavities or stadia.

To characterize the typical size of the fine structure, we compute the spatial intensity correlation functions of high-Q resonances in a local area [54]. The contour lines at half-maximum of the intensity correlation functions [blue dashed lines in Fig. 4(c)] are fitted by an ellipse [red solid lines in Fig. 4(c)]. While the length of the minor axis is consistent for all cavity shapes, the major axis of the ellipse varies significantly with the cavity geometry, and it is defined as the local structure size $s(\vec{r})$.

We calculate the average structure size $\langle s \rangle$ by first averaging over all lasing modes weighted by their intensities, and then averaging over all spatial locations weighted by local intensity [54]. Figure 3 shows that $\langle s \rangle$ [Fig. 3(b), left axis] is strongly correlated with the lasing instabilities in the experiment [Fig. 3(a)]. The irregular pulsations of broad-area semiconductor lasers originate from carrier-
<figure><img src="image_4.png" /><figcaption>FIG. 2. Experimentally measured spatio-temporal lasing dynamics. (a₁-a₅) Streak images of lasing emission from five different cavity geometries. The intensity is normalized by its maximum. The pump current is 500 mA for all cavities, well above their lasing thresholds. Intensity pulsations indicate unstable lasing dynamics. (b₁-b₅) Time-resolved RF-spectra $\hat{S}(f, t_d)$ of laser emission intensities obtained from the streak images. Insets in (b₁,b₂) are the magnifications of the white boxes. The cyan dashed lines indicate the times of the images in (a). RF peaks (bright spots) correspond to intensity pulsations.</figcaption></figure>

induced modulational instability. In a GaAs quantum well, high optical intensity depletes the local gain by spatial hole burning, which increases the refractive index locally. The resulting optical lensing effect and self-focusing lead to the formation of spatial filaments, which are inherently unstable and cause pulsations [17–19]. For cavities with $\langle s \rangle \sim \lambda$, intensity variation on a wavelength scale causes a refractive index change on the same scale, which is too small to focus light [13], thus preventing filamentation and instability. Conversely, a large feature size $\langle s \rangle \gg \lambda$ is more likely to create a lensing effect, which leads to stronger and more frequent pulsations. Therefore, the structure size appears to be a good predictor for the level of lasing stability.

<figure><img src="image_5.png" /><figcaption>FIG. 3. Relation between lasing instabilities and spatial structure of lasing modes. (a) Experimentally measured total RF-power $S_{\text{tot}}$ of intensity pulsations (black circles, left axis) and the frequency of their occurrence (red triangles, right axis) for different cavity shapes. The results of five lasers per geometry are averaged on a logarithmic scale. (b) Numerically calculated local structure size $\langle s \rangle$ in five cavities with boundary roughness, averaged over the cavity area (black squares, left axis). It is determined by the local directionality of light propagation $\langle D_W \rangle$ and $\langle D_R \rangle$ (red symbols, right axis), obtained from wave (+) and ray simulations (×).</figcaption></figure>

The question is what determines the local structure size of lasing modes. The granular structure of lasing modes is formed by the interference of waves propagating in different directions. Thus, the distribution of their directions plays a significant role, which can be unraveled by the spatial Fourier transform of the field profiles. Since it is the *local* directionality that determines the structure size of the intensity distributions, we determine the wave propagation directions in small regions by performing a wavelet transform [60], which can be considered as a local Fourier transform [green solid lines in Fig. 4(c)] [54]. The local directionality is defined as the inverse of the angular spread of wavelet distributions [Fig. 4(a)]. The D-cavities and stadia [Figs. 4(c₁-c₂)] show wave propagation along almost all directions. The interference of these wave components yields a small structure size. For the Limaçon with WGMs [Fig. 4(c₃)], the distribution of wave propagation is more directional parallel to the cavity boundary, leading to a larger structure size perpendicular to the boundary. The ellipses [Fig. 4(c₄)] exhibit a smooth and highly directional distribution, which explains the significantly elongated fine structure of their intensity distributions. Lastly, the square [Fig. 4(c₅)] features four double-peaked lobes along the diagonals, which correspond to the eight plane-wave components of high-Q modes [51, 52]. Despite the low number of plane-wave components, their interference produces nearly isotropic, wavelength-scale intensity grains because the propagation directions are roughly orthogonal.

To quantify how directional the wave propagation is in a local area, we compute the local directionality $D_W(\vec{r})$ as the inverse of the standard deviation of the wavelet distribution squared [54]. Figure 4(d) shows the spatially resolved directionality $D_W(\vec{r})$, averaged over
<figure><img src="image_6.png" /><figcaption>FIG. 4. Numerically simulated local directionality of optical propagation. (a₁-a₃) Illustration of local structure size and local directionality. (a₁) An exemplary high-Q mode in an elliptic cavity features highly directional optical propagation, which yields elongated bright regions. The red arrows denote the structure size. (a₂) The half-maximum contour line of the spatial intensity correlation function C(Δr) (blue dotted lines) for the white box centered at r₀ in (a₁). It is fitted by an ellipse (solid red line) whose major axis is equal to the elongated feature size in (a₁) and defines the local structure size of 0.6 µm (red arrows). (a₃) Polar plot of the amplitude of local Fourier component |w<sub>μ</sub>(r̄₀, θ)| (green), obtained by wavelet transform at the center r̄₀ of the white box in (a₁). The standard deviation of its amplitude squared is 6.7° (black arrows), whose inverse yields the local directionality. (b₁-b₅) Typical high-Q modes of five cavity shapes with surface roughness. The scale bars are 5 µm long. (Insets) Magnifications of the white boxes. (c₁-c₅) The analysis illustrated in (a) is performed with the modes in (b). (d₁-d₅) Maps of local directionality D<sub>W</sub>(r̄) of the simulated lasing modes. The color denotes the local directionality (red: high, blue: low), and the brightness denotes the local intensity of lasing modes. (e₁-e₅) Map of local directionality D<sub>R</sub>(r̄) obtained from ray tracing, showing excellent agreement with wave simulations in (d).</figcaption></figure>

the lasing modes. The D-shaped and stadium cavities [Figs. 4(d₁,d₂)] have lower $D_W(\vec{r})$ than Limaçon resonators [Fig. 4(d₃)], which demonstrates the difference between spatially extended and whispering-gallery modes. The local directionality of the ellipses [Fig. 4(d₄)] is even higher than for the Limaçon, probably due to the integrable ray dynamics of the ellipse, which limits the propagation directions of whispering-gallery trajectories more than for the chaotic trajectories of the Limaçon. The squares [Fig. 4(d₅)] have low and almost uniform $D_W(\vec{r})$ over the cavity area.

These different degrees of local directionality originate from the classical ray dynamics. We hence perform ray tracing simulations for cavities with smooth boundaries and define the local directionality $D_R(\vec{r})$ analogously by sampling the long-lived ray trajectories [54]. Figure 4(e) shows an excellent agreement between the ray and wave simulations. We average the local directionality over the entire cavity area weighted by the ray intensity, and the results are summarized in Fig. 3(b) (right axis). The good agreement of local directionality between ray $\langle D_R \rangle$ and wave simulations $\langle D_W \rangle$ indicates that a cavity with

sufficiently small boundary roughness can be efficiently simulated by the ray tracing of smooth cavities. More importantly, the strong correlation between $\langle D_R \rangle$ and $\langle s \rangle$ confirms that the structure sizes of lasing modes are determined by the local directionality of optical propagation. Hence, our results demonstrate that ray dynamics can be an efficient tool to qualitatively predict the spatio-temporal lasing dynamics.

Apart from the structure size of lasing modes, spatial localization of the modes can also promote the nonlinear processes in the gain medium. Cavities with WGMs like Limaçon or ellipse feature high local optical intensities, and this may facilitate the self-focusing effect and result in stronger output intensity pulsations. Even though simulations with a detailed model of semiconductor carrier dynamics for asymmetric cavities are desirable for a full understanding, our analysis of experiments and passive cavity modes already yields important insights into the relation between the spatial structure of cavity modes and nonlinear lasing dynamics.

To conclude, we establish the resonator geometry as a powerful design parameter to control the spatio-temporal
dynamics of semiconductor microlasers. The lasing dynamics is related to the local directionality of wave propagation, which directly corresponds to the ray dynamics. Our findings enable us to engineer the lasing dynamics by designing the cavity shape based on ray-dynamical principles. In contrast to the design of chaotic microlasers by tailoring temporal oscillation frequencies [39, 61], our approach is based on tailoring the spatial frequencies of lasing modes, which provides a huge and unexplored parameter space. From a practical perspective, customizing the cavity shape enables compact devices to be easily integrated on-chip, in contrast to optical injection and time-delayed feedback [3, 62, 63]. Potential applications are the development of high-power broad-area lasers with stable dynamics and compact lasers for chaos-based applications [63, 64]. Furthermore, the principle of controlling the nonlinear dynamics via its geometry can also find application in other types of lasers such as broad-area VCSELs [65–67] or random lasers [68], as well as in other nonlinear dynamic systems in aerodynamics, fluid dynamics, and plasma physics.

ACKNOWLEDGMENTS

The authors thank Roland Ketzmerick, Jan Wiersig, Ortwin Hess, Stefano Guazzotti, Takahisa Harayama, and Douglas Stone for fruitful discussions.

H.C. and K.K. acknowledge the computational resources provided by the Yale High Performance Computing Cluster (Yale HPC).

The work done at Yale is supported partly by the National Science Foundation under Grant No. ECCS-1953959, and the Office of Naval Research under the Grant No. N00014-221-1-2026.

S.B. acknowledges funding for the Chair in Photonics by Ministère de l'Enseignement supérieur, de la Recherche et de l'Innovation; GDI Simulation; Région Grand-Est; Département Moselle; European Regional Development Fund (ERDF); CentraleSupélec; Fondation CentraleSupélec; and Metz Metropole.

* These authors contributed equally.

† hui.cao@yale.edu

[1] E. Ott, C. Grebogi, and J. A. Yorke, Controlling chaos, Phys. Rev. Lett. **64**, 1196 (1990).

[2] R. Roy, T. Murphy Jr, T. Maier, Z. Gills, and E. Hunt, Dynamical control of a chaotic laser: Experimental stabilization of a globally coupled system, Phys. Rev. Lett. **68**, 1259 (1992).

[3] J. Ohtsubo, *Semiconductor Lasers - Stability, Instability and Chaos*, 3rd ed. (Springer, Heidelberg, 2013).

[4] J. Mørk, J. Mark, and B. Tromborg, Route to chaos and competition between relaxation oscillations for a semiconductor laser with optical feedback, Phys. Rev. Lett. **65**, 1999 (1990).

[5] I. Fischer, O. Hess, W. Elsäßer, and E. Göbel, High-dimensional chaotic dynamics of an external cavity semiconductor laser, Phys. Rev. Lett. **73**, 2188 (1994).

[6] I. Fischer, G. Van Tartwijk, A. Levine, W. Elsäßer, E. Göbel, and D. Lenstra, Fast pulsing and chaotic itinerancy with a drift in the coherence collapse of semiconductor lasers, Phys. Rev. Lett. **76**, 220 (1996).

[7] J. Martin-Regalado, G. Van Tartwijk, S. Balle, and M. San Miguel, Mode control and pattern stabilization in broad-area lasers by optical feedback, Phys. Rev. A **54**, 5386 (1996).

[8] S. K. Mandre, I. Fischer, and W. Elsäßer, Spatiotemporal emission dynamics of a broad-area semiconductor laser in an external cavity: stabilization and feedback-induced instabilities, Opt. Commun. **244**, 355 (2005).

[9] S. Wieczorek, B. Krauskopf, and D. Lenstra, A unifying view of bifurcations in a semiconductor laser subject to optical injection, Opt. Commun. **172**, 279 (1999).

[10] S. Hwang and J. Liu, Dynamical characteristics of an optically injected semiconductor laser, Opt. Commun. **183**, 195 (2000).

[11] B. Krauskopf, S. Wieczorek, and D. Lenstra, Different types of chaos in an optically injected semiconductor laser, Appl. Phys. Lett. **77**, 1611 (2000).

[12] S. Takimoto, T. Tachikawa, R. Shogenji, and J. Ohtsubo, Control of spatio-temporal dynamics of broad-area semiconductor lasers by strong optical injection, IEEE Photon. Technol. Lett. **21**, 1051 (2009).

[13] S. Bittner, S. Guazzotti, Y. Zeng, X. Hu, H. Yılmaz, K. Kim, S. S. Oh, Q. J. Wang, O. Hess, and H. Cao, Suppressing spatio-temporal lasing instabilities with wavechaotic microcavities, Science **361**, 1225 (2018).

[14] K. Kim, S. Bittner, Y. Zeng, S. Guazzotti, O. Hess, Q. J. Wang, and H. Cao, Massively parallel ultrafast random bit generation with a chip-scale laser, Science **371**, 948 (2021).

[15] K. Kim, S. Bittner, Y. Jin, Y. Zeng, S. Guazzotti, O. Hess, Q. J. Wang, and H. Cao, Sensitive control of broad-area semiconductor lasers by cavity shape, APL Photonics **7**, 056106 (2022).

[16] I. Fischer, O. Hess, W. Elsäßer, and E. Göbel, Complex spatio-temporal dynamics in the near-field of a broad-area semiconductor laser, Europhys. Lett. **35**, 579 (1996).

[17] O. Hess and T. Kuhn, Maxwell-bloch equations for spatially inhomogeneous semiconductor lasers. ii. spatiotemporal dynamics, Phys. Rev. A **54**, 3360 (1996).

[18] J. Marciante and G. Agrawal, Spatio-temporal characteristics of filamentation in broad-area semiconductor lasers, IEEE J. Quant. Electron. **33**, 1174 (1997).

[19] J. R. Marciante and G. P. Agrawal, Spatio-temporal characteristics of filamentation in broad-area semiconductor lasers: experimental results, IEEE Phot. Tech. Lett. **10**, 54 (1998).

[20] A. Klaedtke and O. Hess, Ultrafast nonlinear dynamics of whispering-gallery mode micro-cavity lasers, Opt. Express **14**, 2744 (2006).

[21] D. Scholz, H. Braun, U. T. Schwarz, S. Brüninghoff, D. Queren, A. Lell, and U. Strauss, Measurement and simulation of filamentation in (Al,In)GaN laser diodes, Opt. Express **16**, 6846 (2008).

[22] M. Arahata and A. Uchida, Inphase and antiphase dynamics of spatially-resolved light intensities emitted by
a chaotic broad-area semiconductor, IEEE J. Sel. Top. Quant. Electron. 21, 1800609 (2015).

[23] E. G. Altmann, J. S. E. Portela, and T. Tél, Leaking chaotic systems, Rev. Mod. Phys. 85, 869 (2013).

[24] H. Cao and J. Wiersig, Dielectric microcavities: Model systems for wave chaos and non-Hermitian physics, Rev. Mod. Phys. 87, 61 (2015).

[25] R. K. Chang and A. J. Campillo, *Optical processes in microcavities*, Vol. 3 (World scientific, 1996).

[26] A. D. Stone, Wave-chaotic optical resonators and lasers, Phys. Scr. 2001, 248 (2001).

[27] H. E. Türeci, H. G. L. Schwefel, P. Jacquod, and A. D. Stone, Modes of wave-chaotic dielectric resonators, Prog. Opt. 47, 75 (2005).

[28] T. Harayama and S. Shinohara, Two-dimensional microcavity lasers, Laser Photonics Rev. 5, 247 (2011).

[29] A. Levi, R. Slusher, S. McCall, J. Glass, S. Pearton, and R. Logan, Directional light coupling from microdisk lasers, Appl. Phys. Lett. 62, 561 (1993).

[30] J. U. Nöckel, A. D. Stone, G. Chen, H. L. Grossman, and R. K. Chang, Directional emission from asymmetric resonant cavities, Opt. Lett. 21, 1609 (1996).

[31] G. Chern, H. Tureci, A. D. Stone, R. Chang, M. Kneissl, and N. Johnson, Unidirectional lasing from InGaN multiple-quantum-well spiral-shaped micropillars, Appl. Phys. Lett. 83, 1710 (2003).

[32] H. G. L. Schwefel, H. E. Tureci, A. D. Stone, and R. K. Chang, Progress in asymmetric resonant cavities: Using shape as a design parameter in dielectric microcavity lasers, Optical Processes in Microcavities, World Scientific, Singapore (2003).

[33] H. G. L. Schwefel, N. B. Rex, H. E. Tureci, R. K. Chang, A. D. Stone, T. Ben-Messaoud, and J. Zyss, Dramatic shape sensitivity of directional emission patterns from similarly deformed cylindrical polymer lasers, J. Opt. Soc. Am. B 21, 923 (2004).

[34] T. Fukushima, T. Tanaka, and T. Harayama, Unidirectional beam emission from strained InGaAsP multiple-quantum-well quasistadium laser diodes, Appl. Phys. Lett. 86, 171103 (2005).

[35] S. Sunada, T. Harayama and K. S. Ikeda, Multimode lasing in two-dimensional fully chaotic cavity lasers, Phys. Rev. E 71, 046209 (2005).

[36] T. Harayama, T. Fukushima, S. Sunada, and K. S. Ikeda, Chaos and multi-attractors in fully chaotic 2D microcavity lasers, Prog. Theor. Phys. Suppl. 166, 104 (2007).

[37] M. Choi, T. Fukushima, and T. Harayama, Alternate oscillations in quasistadium laser diodes, Phys. Rev. A 77, 063814 (2008).

[38] S. Shinohara, T. Fukushima, S. Sunada, T. Harayama, K. Arai, and K. Yoshimura, Anticorrelated bidirectional output of quasistadium-shaped semiconductor microlasers, Opt. Rev. 21, 113 (2014).

[39] C.-G. Ma, J.-L. Xiao, Z.-X. Xiao, Y.-D. Yang, and Y.-Z. Huang, Chaotic microlasers caused by internal mode interaction for random number generation, Light: Sci. Appl. 11, 187 (2022).

[40] K. Kim, S. Bittner, Y. Jin, Y. Zeng, Q. Wang, and H. Cao, Spatiotemporal lasing dynamics in a Limaçon-shaped microcavity, Opt. Lett. 48, 574 (2023).

[41] L. A. Bunimovich, On the ergodic properties of nowhere dispersing billiards, Commun. Math. Phys. 65, 295 (1979).

[42] S. Ree and L. E. Reichl, Classical and quantum chaos in

a circular billiard with a straight cut, Phys. Rev. E 60, 1607 (1999).

[43] M. V. Berry, Regular and irregular semiclassical wavefunctions, J. Phys. A 10, 2083 (1977).

[44] A. I. Shnirel'man, Ergodic properties of eigenfunctions, Uspekhi Mat. Nauk 29, 181 (1974).

[45] R. Ketzmerick, K. Clauß, F. Fritzsch, and A. Bäcker, Chaotic resonance modes in dielectric cavities: Product of conditionally invariant measure and universal fluctuations, Phys. Rev. Lett. 129, 193901 (2022).

[46] M. You, D. Sakakibara, K. Makino, Y. Morishita, K. Matsumura, Y. Kawashima, M. Yoshikawa, M. Tonosaki, K. Kanno, A. Uchida, et al., Universal single-mode lasing in fully chaotic billiard lasers, Entropy 24, 1648 (2022).

[47] M. Robnik, Classical dynamics of a family of billiards with analytic boundaries, J. Phys. A 16, 3971 (1983).

[48] J. Wiersig and M. Hentschel, Combining directional light output and ultralow loss in deformed microdisks, Phys. Rev. Lett. 100, 033901 (2008).

[49] H. R. Dullin and A. Bäcker, About ergodicity in the family of limaçon billiards, Nonlinearity 14, 1673 (2001).

[50] M. V. Berry, Regularity and chaos in classical mechanics, illustrated by three deformations of a circular 'billiard', Eur. J. Phys. 2, 91 (1981).

[51] S. Bittner, E. Bogomolny, B. Dietz, M. Miski-Oglu, and A. Richter, Experimental observation of localized modes in a dielectric square resonator, Phys. Rev. E 88, 062906 (2013).

[52] S. Bittner, E. Bogomolny, B. Dietz, M. Miski-Oglu, and A. Richter, Dielectric square resonator investigated with microwave experiments, Phys. Rev. E 90, 052909 (2014).

[53] Y.-D. Yang and Y.-Z. Huang, Mode characteristics and directional emission for square microcavity lasers, J. Phys. D 49, 253001 (2016).

[54] See Supplemental Material [url] for a detailed description of laser characterization, measurement of lasing dynamics, wave simulations, calculation of local directionality and structure size, and ray tracing algorithm, which includes Refs. [69, 70].

[55] S. Bittner, K. Kim, Y. Zeng, Q. J. Wang, and H. Cao, Spatial structure of lasing modes in wave-chaotic semiconductor microcavities, N. J. Phys. 22, 083002 (2020).

[56] L. Ge, Y. D. Chong, and A. D. Stone, Steady-state ab initio laser theory: Generalizations and analytic results, Phys. Rev. A 82, 063824 (2010).

[57] S. F. Liew, L. Ge, B. Redding, G. S. Solomon, and H. Cao, Pump-controlled modal interactions in microdisk lasers, Phys. Rev. A 91, 043828 (2015).

[58] A. Cerjan, B. Redding, L. Ge, S. F. Liew, H. Cao, and A. D. Stone, Controlling mode competition by tailoring the spatial pump distribution in a laser: a resonance-based approach, Opt. Express 24, 26006 (2016).

[59] A. Cerjan, S. Bittner, M. Constantin, M. Guy, Y. Zeng, Q. J. Wang, H. Cao, and A. D. Stone, Multimode lasing in wave-chaotic semiconductor microlasers, Phys. Rev. A 100, 063814 (2019).

[60] I. Daubechies, *Ten Lectures on Wavelets* (Society for Industrial and Applied Mathematics, 1992).

[61] J.-C. Li, J.-L. Xiao, Y.-D. Yang, and Y.-Z. Huang, Random bit generation based on self-chaotic microlasers with enhanced chaotic bandwidth, arXiv preprint arXiv:2301.00111 (2022).

[62] M. C. Soriano, J. García-Ojalvo, C. R. Mirasso, and
I. Fischer, Complex photonics: Dynamics and applica-
tions of delay-coupled semiconductors lasers, Rev. Mod.
Phys. **85**, 421 (2013).

[63] M. Sciamanna and K. A. Shore, Physics and applications of laser diode, Nat. Photonics **9**, 151 (2015).

[64] X.-Q. Qi and J.-M. Liu, Photonic microwave applications of the dynamics of semiconductor lasers, IEEE J. Sel. Top. Quant. Electron. **17**, 1198 (2011).

[65] A. Brejnak, M. Gebski, A. K. Sokól, M. Marciniak, M. Wasiak, J. Muszalski, J. A. Lott, I. Fischer, and T. Czyszanowski, Boosting the output power of large-aperture lasers by breaking their circular symmetry, Optica **8**, 1167 (2021).

[66] S. Bittner and M. Sciamanna, Complex nonlinear dynam-

ics of polarization and transverse modes in a broad-area
VCSEL, APL Photonics **7**, 126108 (2022).

[67] O. Alkhazragi, M. Dong, L. Chen, D. Liang, T. K. Ng, J. Zhang, H. Bagci, and B. S. Ooi, Modifying the coherence of vertical-cavity surface-emitting lasers using chaotic cavities, Optica **10**, 191 (2023).

[68] S. Bittner, S. Knitter, S. F. Liew, and H. Cao, Random-laser dynamics with temporally modulated pump, Phys. Rev. A **99**, 013812 (2019).

[69] V. F. Lazutkin, The existence of caustics for a billiard problem in a convex domain, Math. USSR Izv. **7**, 185 (1973).

[70] J. B. Keller and S. I. Rubinow, Asymptotic solution of eigenvalue problems, Ann. Phys. **9**, 24 (1960).
# Supplementary Material: Impact of cavity geometry on microlaser dynamics

Kyungduk Kim*,<fnref n="1" /> Stefan Bittner*,<fnref n="2" />,<fnref n="3" /> Yuhao Jin,<fnref n="4" /> Yongquan Zeng,<fnref n="4" /> Qi Jie Wang,<fnref n="4" /> and Hui Cao†<fnref n="1" />

## A. Laser characterization resonators.

### 1. Device fabrication

We fabricate edge-emitting semiconductor microcavity lasers on a commercial epwafer with a GaAs/AlGaAs quantum well (Q-Photonics QEWLD-808). First, the bottom metal contact is deposited on the backside of the wafer. Next, lateral boundaries of different cavity shapes are defined via UV-lithography and inductively coupled plasma (ICP) dry etching. The etch depth is 3.5 µm, reaching the bottom cladding layer to ensure strong light confinement within the cavity by a high refractive index difference at the sidewalls. The sidewalls are almost vertical and have a small but non-negligible surface roughness (see Ref. [1]). Finally, the top metal (Ti/Au) contacts for electrical pumping are deposited. The top contacts are withdrawn by 6 µm from the cavity edges in order to avoid misalignment during photolithography, causing the metal hanging down over the sidewall and blocking the edge emission (see Ref. [2] for details of the fabrication process).

Multiple devices for each of the five different geometries (D-shaped, stadium, Limaçon, ellipse, and square) are fabricated concurrently on the same epwafer (see Fig. S1). The D-cavity is a circle with a straight cut, where the radius is $R_D = 100$ µm, and the distance from the circle center to the cut is set to $R_D/2$ to render the ray dynamics as chaotic as possible. The stadium consists of a square with side length $a = 119$ µm between two semicircles of radius $a/2$. The boundary of the Limaçon cavity [3–5] is defined in polar coordinates by $ρ(φ) = ρ_0(1 + ϵ \cos φ)$ where $φ$ is the azimuthal angle, $ρ_0 = 86$ µm the mean radius and $ϵ = 0.42$ the deformation parameter which is chosen to obtain unidirectional emission [4]. The ellipse cavities have an aspect ratio of $b/a = 2$ with a minor (major) diameter of $a = 127$ µm ($b = 254$ µm). Lastly, the square microlasers have a side length of 159 µm. These different cavity geometries result in approximately equal cavity areas ($2.53 × 10^4$ µm<fnref n="2" />).

The guiding layer of the quantum well wafer supports a single vertical excitation mode, and thus the optical field propagates in the plane of the resonators with a phase velocity $c/n_{eff}$, where c is the speed of light in vacuum and $n_{eff} = 3.37$ the effective refractive index. The laser emission has transverse electric (TE) polarization; that is, the electric field is parallel to the plane of the

### 2. Device testing

We perform experiments on five microlasers for each cavity shape. The microlasers are contacted with a Tungsten needle (Quater Research H-20242) and pumped electrically at room temperature by a diode driver (DEI Scientific PCX-7401). We use 2 µs long current pulses with a repetition rate lower than 1 Hz to reduce heating. All the spatio-temporal measurements are done at the same current of 500 mA (corresponding to about 2.0 kA/cm<fnref n="2" />) for all cavity shapes.

The optical spectra of lasing emission, with a central wavelength of about 800 nm, are measured with an imaging spectrometer (Acton SP300i) equipped with an intensified CCD camera (ICCD, Andor iStar DH312T-18U-73). For measurements of time-resolved optical spectra, we time-gate the ICCD camera with a temporal resolution of 50 ns.

To investigate the spatio-temporal dynamics of the microlasers, one facet of a microcavity (marked by cyan dashed lines in Fig. S1) is imaged onto the entrance slit of a streak camera using a microscope objective (20×, NA = 0.4) and a tube lens. Then the near-field laser emission is measured by a streak camera (Hamamatsu C5680) with a fast single-sweep unit (M5676) (see Ref. [6] for further details of the setup). We measure 10 ns long streak images with a temporal resolution of about 30 ps.

### 3. Lasing thresholds

Figure S2 shows the lasing thresholds for different resonator geometries. The LI curves of all the measured cavities exhibit a sharp transition in the slope, where one example for a square cavity is shown in Fig. S2(a). Lasing thresholds for different cavity shapes, each averaged over five devices, are presented in Fig. S2(b). The fluctuations of the threshold amongst cavities of the same shape are very small, indicating the generally good reproducibility and consistency of fabrication and measurements. The typical lasing thresholds range from 80 to 150 mA (from 0.32 to 0.60 kA/cm<fnref n="2" />), strongly depending on the cavity geometry. The square, ellipse, and Limaçon cavities

<footnote n="1">¹Department of Applied Physics, Yale University, New Haven, Connecticut 06520, USA</footnote>

<footnote n="2">²Université de Lorraine, CentraleSupélec, LMOPS, 2 rue Edouard Belin, Metz 57070, France</footnote>

<footnote n="3">³Chair in Photonics, CentraleSupélec, LMOPS, 2 rue Edouard Belin, Metz 57070, France</footnote>

<footnote n="4">⁴Center for OptoElectronics and Biophotonics, School of Electrical and Electronic Engineering, School of Physical and Mathematical Science, and Photonics Institute, Nanyang Technological University, 639798, Singapore</footnote>
<figure><img src="image_7.png" /><figcaption>FIG. S1. Scanning electron microscope images of GaAs microlasers with D-shaped (a), stadium (b), Limaçon (c), ellipse (d), and square (e) dielectric cavities in perspective view. The scale bars are 50 µm long. The yellow area represents the top metal contact for electrical pumping. The cyan dashed lines mark the facets imaged onto the streak camera.</figcaption></figure>

have the low thresholds, followed by the stadium, and
the D-cavities have the highest thresholds.

B. Lasing dynamics

1. *Optical spectrum*

Figure S3 shows the spectrochronogram (time-resolved spectrum) of a square microlaser which is measured with a tomographic technique: the ICCD connected to the spectrometer is gated on during a 50 ns long interval to measure spectra with a high temporal resolution, and the time delay time $t_d$ of the gate interval with respect to the start of the pump pulse is successively increased from one pulse to the next (cf. Ref. [6]). The laser remains multimode at any point in time, though the instantaneous number of lasing peaks is smaller than that in the time-integrated spectra (the right panel of Fig. S3). Similar spectrochronograms are obtained for the other microlaser geometries (cf. Refs. [5–7]). It should be emphasized that multimode lasing is observed for all cavity shapes, and the interaction of multiple lasing modes is an important

ingredient in the development of lasing instabilities [8, 9].

The spectrum changes continuously over the course of the current pulse with a timescale of the order of 100 ns, which is much longer than the intrinsic time scales of the nonlinear lasing dynamics discussed in Section B 3. Due to Joule heating, the gain spectrum redshifts, thus the lasing peaks at shorter wavelength disappear, and new peaks emerge at longer wavelengths [6]. Consequently, the lasing state gradually evolves over the course of the pulse. This enables us to investigate an ensemble of different possible dynamical states by recording streak images at different times during the pulse, where the time window (10 ns) of individual streak images is considerably shorter than the time scale of the heating-induced drift of the lasing spectrum.

2. *Spatio-temporal dynamics*

We measure the evolution of spatio-temporal lasing dy-
namics during a pump current pulse of 2 µs length. We
recorded single-sweep streak images of 10 ns length (tem-
poral resolution about 30 ps). While repeatedly pumping
the lasers under identical conditions, we gradually shift

<figure><img src="image_8.png" /><figcaption>FIG. S2. Experimentally measured lasing thresholds of asym-
metric cavities with deformed geometry. (a) An exemplary LI
curve of a square microcavity laser: emission power obtained
from the time-integrated lasing spectra as a function of in-
jected current. The lasing threshold is 80 mA (blue arrow).
(b) Lasing thresholds for the five cavity geometries. The er-
rorbars represent the variation between five different devices
of the same geometry.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>FIG. S3. Measured spectrochronogram of a square micro-
cavity laser at 500 mA pump current. The horizontal axis
indicates the time delay $t_d$ of the gate interval with respect
to the start of the 2-µs-long current pulse. The white scale
bar indicates the time resolution of 50 ns. The right panel
shows the emission spectrum integrated over the entire pulse,
exhibiting multiple overlapping lasing peaks.</figcaption></figure>
the gated time window of measurement to scan the entire emission pulse. More specifically, the delay time $t_d$ between the streak image and the beginning of the current pulse is varied in the range of 0.3 – 1.9 µs [6]. Each streak image provides the spatio-temporal intensity pattern of lasing emission $I(x, t_d, t_d + t)$, where $t = 0–10$ ns is the time during one image. For each microlaser, 161 consecutive images covering a total time of 1.61 µs are measured, omitting only the transient dynamics at the start and end of the 2-µs-long pulses. These images are concatenated in the time domain to obtain microsecond-long traces. To account for cavity-to-cavity variations, we obtain streak images for five different devices for each cavity geometry. The examples for each cavity geometry are shown in Fig. 2(a) of the main text.

Several dynamic processes in our microlasers contribute to the spatio-temporal fluctuations observed in the streak images [5]. First, there are relatively strong intensity fluctuations with a rather narrow radio-frequency (RF) band, which appear as peaks in the short-term RF-spectra for some delay times $t_d$ during the current pulses and can be spatially localized as well. Second, we observe a ubiquitous spatio-temporal speckle caused by the beating of many lasing modes with different optical frequencies and spatial patterns [2, 10] and/or amplified spontaneous emission with a broad spectrum. This spatio-temporal speckle is common for multimode lasers like the ones investigated here but is not the focus of the current study. Hence we will separate the first type of fluctuations, resulting from filamentation and lasing instabilities, from the second type due to spatio-temporal interference.

### 3. Radio-frequency spectrum

Here we separate intensity fluctuations of different origins in the measured streak images. First, each streak image $I(x, t_d, t_d + t)$ is normalized such that $\langle I(x, t_d, t_d + t) \rangle_{x,t} = 1$ to enable a quantitative comparison of the fluctuation strength of different measurements. Next, we calculate the spatially-resolved Fourier transform (FT) of the streak image,

$$ \tilde{I}(x, f, t_d) = \int_0^T dt I(x, t_d, t_d + t)e^{-2\pi i f t}, \quad (S1) $$

where $f$ is the frequency in the range up to 12 GHz, which is limited by the temporal resolution of the streak image, and $T = 10$ ns is the length of a single streak image. Then we spatially average it to obtain the RF power spectrum

$$ S(f, t_d) = \langle |\tilde{I}(x, f, t_d)|^2 \rangle_x \quad (S2) $$

of the streak image.

The spatio-temporal beating of lasing modes and/or ASE produces a very broad range of RF components. The highest possible frequency, given by the width of

the optical spectrum of lasing emission, is of the order of hundreds of GHz. Hence, the measured spatio-temporal speckle is limited by the temporal resolution of our streak camera. The Lorentzian-shaped point spread function for the temporal response of our streak camera [10] results in a continuous RF spectrum that decays exponentially with frequency. The lasing instabilities manifest in relatively narrow peaks on top of it.

We develop a process to separate the RF peaks from the broadband contribution of spatio-temporal speckle. Figure S4(a) shows an exemplary streak image of a square microcavity laser, which contains intensity fluctuations from lasing instabilities as well as spatio-temporal beating. Temporal Fourier transform gives the spatially-resolved RF spectrum in Fig. S4(b). Since the RF power of the spatio-temporal speckle has a comparable magnitude to that of the lasing pulsations, it is not easy to visually separate the two. The spatially-averaged RF-spectrum [blue curve in Fig. S4(c)] displays a peak at frequency 2.2 GHz due to almost regular intensity oscillations.

To determine the RF-power level of the spatio-temporal speckle, we plot the spatially-averaged RF-spectra $S(f, t_d)$ at every time $t_d$ together as shown in Fig. S4(c). All spectra feature a continuous signal decaying exponentially with frequency, which originates from the resolution-limited spatio-temporal speckle [10]. We perform a linear fit of the exponential decay in a semilogarithmic scale. It yields the broadband continuous RF-signal $S_b(f)$, which we consider as the RF-power level of spatio-temporal speckle.

To separate the RF peaks above $S_b(f)$, we subtract this fitted exponentially-decaying power spectrum by

$$ \hat{S}(f, t_d) = S(f, t_d) - S_b(f), \quad (S3) $$

and setting the negative values of $\hat{S}(f, t_d)$ to 0. The subtracted RF-spectrum $\hat{S}(f, t_d)$ is presented in Fig. S4(e), in which only the oscillatory peaks caused by lasing instabilities are left and can thus be well discerned. Figure 2(b) in the main text shows the subtracted RF-spectra $\hat{S}(f, t_d)$ for each cavity geometry.

### 4. Lasing instabilities

For a more quantitative evaluation of the RF-spectra, we consider two quantities characterizing different aspects in the following. First, we look at the integrated RF power,

$$ S_{\text{tot}} = \sum_{f, t_d} \hat{S}(f, t_d). \quad (S4) $$

Second, we calculate the frequency of occurrence of instabilities via the participation ratio of the RF spectra,

$$ S_{\text{PR}} = \frac{\langle \hat{S}(f, t_d) \rangle_{f, t_d}^2}{\langle [\hat{S}(f, t_d)]^2 \rangle_{f, t_d}}. \quad (S5) $$
<figure><img src="image_10.png" /><figcaption>FIG. S4. Subtraction of RF-spectrum to highlight frequency components from lasing instabilities. (a) Streak image of lasing emission from a square microcavity for a time window of 10 ns length, measured at the delay time $t_d = 1.20~\mu s$ from the beginning of a 2-µs-long pulse. (b) Spatially resolved RF spectrum $|\tilde{I}(x, f, t_d)|^2$ of the emission intensity, obtained from temporal Fourier transform of (a). (c) Spatially-averaged RF-spectra $S(f, t_d)$ at all delay times $t_d$ are plotted together (gray curves). The blue curve highlights $S(f, t_d)$ from (b), which exhibits a peak standing out at 2.2 GHz. The red dashed line indicates the fit $S_b(f)$ of the exponential decay of the RF-power from the resolution-limited spatio-temporal speckle. The RF peaks above the red line result from lasing instabilities. (d) Time-resolved RF-spectrum $S(f, t_d)$ before subtracting $S_b(f)$. (e) Modified RF-spectrum $\hat{S}(f, t_d)$ obtained by subtracting $S_b(f)$ from $S(f, t_d)$. The peaks due to lasing instabilities are clearly visible once the contributions from spatio-temporal beating are removed.</figcaption></figure>

These two quantities are extracted from measurements with five different microlasers per geometry, and the results averaged in logarithmic scale are presented in Fig. 3(a) of the main text.

### C. Wave simulations

#### 1. Resonances

We model the laser cavities as two-dimensional (2D) dielectric resonators, that is, the passive cavity modes are the solutions of the 2D scalar Helmholtz equation

$$ \{\Delta + n^2(\vec{r})k^2\}\Psi(\vec{r}) = 0, \quad (S6) $$

with outgoing-wave boundary conditions [11, 12], where $\vec{r} = (x, y)$ are the coordinates in the plane of the cavity, $n(\vec{r})$ is the effective refractive index structure of the cavity, $n(\vec{r}) = n_{\text{eff}} = 3.37$ inside and $n(\vec{r}) = 1$ outside of the cavity, and $\Psi$ corresponds to the z-component of the magnetic field, $H_z$, for TE-polarized modes. The solutions are computed numerically using the COMSOL eigenfrequency solver module, where the outgoing-

wave boundary conditions are implemented via perfectly matched layers at the boundaries of the computational domain. A real part of the resonance wave numbers $k$ corresponds to the resonance frequency and the imaginary part to the decay rate (inverse of the lifetime), hence the quality factors are given by $Q = -\text{Re}(k)/[2\text{Im}(k)]$.

The wave simulations are performed for cavities with ten times smaller linear dimensions (area of about 253 µm²) due to computational constraints. Since both the simulated and the fabricated cavities are well within the semiclassical regime $kR \gg 1$, where $R$ is a typical linear dimension (e.g., the radius of the D-cavity), the simulation results of smaller cavities can nonetheless be considered representative of the actual, larger cavities. Exemplary high-Q resonances of five different cavity geometries with smooth boundaries are shown in Fig. 1(b₁-b₅) of the main text.

#### 2. Surface roughness

Since surface roughness is inevitable in fabricated microwaves, we calculate the resonances of passive res-
onators with boundary roughness. The roughness model we use is similar to the one in Refs. [6, 13]. The cavity boundaries are perturbed by a superposition of high-order harmonics with random phase and amplitude,

$$ \Delta R(s) = C \sum_{m=m_1}^{m_2} a_m \cos(2\pi ms + \theta_m), \quad (S7) $$

where $s \in [0, 1)$ is the normalized arc length coordinate along the cavity boundary. At every location $s$, the perturbation is added in the direction perpendicular to the boundary. The random variables $a_m$ determine the amplitude of harmonic perturbation, which follow a uniform distribution in the range of $a_m \in [-A, A]$ nm, where $A$ is set to 10 nm. The random phases $\theta_m$ are uniformly distributed in the range $[0, 2\pi)$. The normalization constant $C$ is given by $\sqrt{6/(m_2 - m_1 + 1)}$ so the root-mean-square of the fluctuation $\Delta R(s)$ is $A$.

The lowest-order harmonics $m_1$ (with the longest period) is determined by $s_0/L_{max}$, where $s_0$ is the perimeter of the cavity boundary and $L_{max}$ is the maximal length scale of boundary fluctuations parallel to the cavity boundary. Based on high-resolution SEM images of experimental cavities, we set $L_{max} = 1$ µm. The highest harmonics $m_2$ (with the shortest period) is determined by $s_0/L_{min}$, where $L_{min}$ is the shortest scale of boundary modulation considered. We set its value to half of the in-medium wavelength, $L_{min} = \lambda/(2n_{eff}) = 0.12$ µm, as any perturbation on length scales much finer than the wavelength cannot be resolved by light. Given that the perimeters $s_0$ of the simulated cavities of all five different shapes are roughly 60 µm, the approximate range of the harmonics $m$ is from $m_1 = 60$ to $m_2 = 500$, where the exact values of $m_1$ and $m_2$ depend on the cavity shape.

### 3. Lasing modes

An important aspect of microcavity lasers is the mode competition that determines the number of lasing modes. While simulating the dynamics of our asymmetric microcavity lasers is extremely computationally expensive [6], the Steady-state ab-initio Lasing Theory (SALT) [14, 15] has proven useful to estimate the number of lasing modes of such lasers, although it is strictly speaking applicable only to lasers that have reached a steady-state. As our microcavities have relatively high Q-factors, we employ the single-pole approximation to SALT (SPA-SALT), which assumes that each lasing mode is represented by a single passive cavity mode (see Refs. [1, 13–15] for details).

We perform SPA-SALT simulations for three different realizations of rough boundaries for each cavity shape. All high-Q resonances within the spectral range of 797–803 nm are considered. The first lasing threshold is inversely proportional to the quality factor $Q_{max}$ of the most long-lived mode. Figure S5(a) shows the computed lasing thresholds for the five cavity shapes. Limaçon,

ellipse, and square cavities have the lowest thresholds, whereas those of D-cavity and stadium are significantly higher. This trend is consistent with the experimental data [Fig. S2(b)], and demonstrates how strongly the cavity shape affects the lasing threshold.

The onset of further lasing modes is determined by a combination of their Q-factors compared to the first lasing mode, which saturates the gain, and the spatial overlap of their intensity distributions. The actual number of lasing modes and the powers of lasing modes depend on spatial gain-competition effects, which are taken into account by the SPA-SALT simulations. We perform SPA-SALT simulations at 10 times the lasing threshold even though we pump less strongly in experiments. This is because the microcavities in simulation have much smaller sizes and hence a smaller number of lasing modes than those in experiments. Figure S5(b) shows the computed number of lasing modes for the five cavity shapes. All cavity shapes exhibit multimode lasing, with a slightly different number of lasing modes at 10 times their lasing thresholds. Note that individual lasing modes have different powers, given by the SPA-SALT simulations.

<figure><img src="image_11.png" /><figcaption>FIG. S5. Lasing behaviors simulated by SPA-SALT. (a) Lasing thresholds and (b) the number of lasing modes at ten times the lasing threshold for five different cavity geometries. The errorbars denote the variation from three different realizations of surface roughness for each geometry.</figcaption></figure>

### 4. Spatial Fourier transform

The main objective of wave simulations is to investigate how the cavity shape affects the optical propagation directions inside the cavity. A straightforward approach is performing the Fourier transform of a passive cavity resonance. The spatial Fourier transform of its field distribution $\Psi(\vec{r})$ is

$$ \tilde{\Psi}(\vec{k}_M) = \int d\vec{r} \Psi(\vec{r}) \exp[-i(\vec{k}_M \cdot \vec{r})], \quad (S8) $$

where $\vec{k}_M = (k_x, k_y)$ is the in-medium wave vector, and the integral is over the interior of the cavity, corresponding to an expansion of the field distribution in plane-wave components.
<figure><img src="image_12.png" /><figcaption>FIG. S6. Spatial Fourier transform of resonances in (a) D-cavity, (b) stadium, (c) Limaçon, (d) ellipse and (e) square micro-cavities with surface roughness. The spatial distributions of field amplitude of the modes are shown in Figs. 4(b₁-b₅) of the main text. The chaotic cavities (a-c) have uniformly distributed plane-wave components in all directions, while ellipse (d) and square (e) resonators feature dominant propagation in vertical and diagonal directions, respectively. The Fourier transform, however, does not include information about spatially-resolved propagation directions.</figcaption></figure>

The wave-vector distributions of the resonances shown in Figs. 4(b₁-b₅) of the main text are presented in Fig. S6. In general, all plane-wave components have the same $|\vec{k}_M| = 2\pi n/\lambda$, where n is the refractive index and $\lambda$ is the resonant wavelength, and hence all the wave-vector distributions in Fig. S6 are localized on a circle with that radius. For the D-cavity, stadium, and Limaçon resonators [Figs. S6(a-c)], the circle is covered quite homogeneously, indicating that there is no preferred propagation direction of the underlying ray trajectories. This is not surprising in view of their chaotic ray dynamics [cf. Figs. 1(a₁-a₃)]. For the ellipse [Fig. S6(d)], the components propagating in the vertical direction (parallel to the major axis) are stronger. This is due to the vertically elongated shape of the ellipse: the whispering gallery trajectories spend more time propagating parallel to the approximately vertical boundaries since these are longer compared to the approximately horizontal boundary parts near the vertices [see Fig. 4(a₄)]. For the square [Fig. S6(e)], the wave-vector components are stronger in the diagonal directions $\theta = \pm 45^\circ$ and $\pm 135^\circ$ since trajectories traveling in these directions are most strongly confined by TIR. For a smooth square, the high-Q modes consist of exactly 8 plane-wave components with directions close to the diagonals [16, 17], but the surface roughness creates additional propagation directions via scattering and thus forms the four clusters of peaks around the diagonals that we observe in Fig. S6(e).

The main drawback of the spatial Fourier transform, in particular for analyzing resonances of wave-chaotic cavities, is that the information about the propagation directions comes without spatial resolution. Hence, in order to obtain the propagation directions locally in the cavity, we apply a wavelet transform.

## D. Local directionality and structure size

### 1. Wavelet transform

We calculate the wavelet transform of spatial field distribution $\Psi_\mu$ of the lasing mode $\mu$, which can be considered a local Fourier transform [18], as

$$w_{\mu}(\vec{r}, \theta) = \int d\vec{r}' \Psi_{\mu}^{*}(\vec{r}') \hat{\Psi}(\vec{r}' - \vec{r}, \theta), \quad (S9)$$

where $\hat{\Psi}$ is the Morlet wavelet [Fig. S7(b)]

$$\hat{\Psi}(\vec{r}, \theta) = e^{-ik_M x_R - x_R^2/(2\sigma^2) - y_R^2/\sigma^2}. \quad (S10)$$

The rotated coordinate frame $(x_R, y_R)$ is given by

$$\begin{aligned} x_R &= x \cos \theta + y \sin \theta \\ y_R &= -x \sin \theta + y \cos \theta \end{aligned} \quad (S11)$$

where $\theta$ is the azimuthal angle with respect to the x-axis and $k_M = 2\pi n/\lambda$ is the in-medium wave number of the mode $\Psi_\mu$. The width $\sigma$ of the Morlet wavelet is $\sigma = 4\lambda/n = 0.95$ µm, which is four times the in-medium wavelength of resonant modes. Therefore, $|w_\mu(\vec{r}, \theta)|$ yields the amplitude of waves propagating in direction $\theta$ in a small region of approximate diameter $2\sigma$ around $\vec{r}$.

Figure S7 exemplifies the wavelet transform using a wave function of an ellipse resonator with surface roughness. The field amplitude is shown in Fig. S7(a) and the Morlet wavelet for $\theta = 0^\circ$ in Fig. S7(b). The wavelet transform for $\theta = 0^\circ$, that is, horizontal wave propagation, is shown in Fig. S7(c). It features high amplitudes only close to the two vertices of the ellipse where the whispering gallery trajectories propagate horizontally [cf. Fig. 1(a₄)]. In contrast, the wavelet transform for vertical propagation ($\theta = 90^\circ$) in Fig. S7(d) features high intensity near the vertically extended parts of the boundaries where ray trajectories propagate in this direction. This demonstrates how the wavelet transform gives spatially
<figure><img src="image_13.png" /><figcaption>FIG. S7. Wavelet transform of a cavity resonance. (a) Field amplitude of a high-Q mode ($\lambda = 797.5$ nm, $Q = 1.56 \times 10^4$) in an ellipse resonator with rough boundary. (b) Real part of the Morlet wavelet $\Psi$ with orientation $\theta = 0^\circ$. Its standard deviation $\sigma$ is 0.95 $\mu$m, four times the in-medium wavelength. The imaginary part of $\Psi$ is similar to its real part, except for a phase shift of $\pi/2$. (c) Wavelet transform $|w_\mu(\vec{r}, \theta)|$ of the wave function for $\theta = 0^\circ$ and (d) for $\theta = 90^\circ$, where the orientation is indicated by a white arrow. The wavelet transform reveals the dominant propagation directions with spatial resolution.</figcaption></figure>

resolved information of the wave propagation directions, which corresponds well to the classical dynamics, in contrast to the spatial Fourier transform shown in Fig. S6(d).

## 2. Local directionality

We calculate the spatially-resolved local directionality of wave propagation over the entire cavity area. First, we calculate the integrated wavelet distribution for all lasing modes $\mu$,

$$W(\vec{r}, \theta) = \sum_{\mu} P_{\mu} |w_{\mu}(\vec{r}, \theta)|^2, \quad (S12)$$

and the average intensity profile

$$I(\vec{r}) = \sum_{\mu} P_{\mu} \int |w_{\mu}(\vec{r}, \theta)|^2 d\theta. \quad (S13)$$

Here $P_\mu$ is the lasing mode power from SPA-SALT, so the lasing modes with higher power contribute more to

the average. The local angular spread of propagation is then defined as

$$\sigma_W(\vec{r}) = \min_{\theta_0} \sqrt{\frac{\int_{-\frac{\pi}{2}}^{+\frac{\pi}{2}} \theta^2 [W(\vec{r}, \theta + \theta_0) + W(\vec{r}, \theta + \theta_0 - \pi)] d\theta}{\int_{-\pi}^{+\pi} W(\vec{r}, \theta) d\theta}} \quad (S14)$$

Here $\theta_0$ indicates the orientation to calculate the angular spread, and the two terms in the numerator indicate the forward and backward half of the angular domain with respect to $\theta_0$. We search for the direction $\theta_0$ in which the angular spread is minimal. It corresponds to the dominant propagation direction. Equation (S14) yields the minimal angular spread in terms of the standard deviation of the wavelet distribution as a function of the polar angle, and the local directionality is defined as its inverse,

$$D_W(\vec{r}) = \frac{1}{\sigma_W(\vec{r})}. \quad (S15)$$

Maps of the local directionality $D_W(\vec{r})$ for cavities with rough boundaries and different geometries are presented in Figs. 4(d1-d5) of the main text.

Finally, we calculate the overall local directionality of a cavity by averaging $D_W(\vec{r})$ over the position weighted by the local intensity,

$$\langle D_W \rangle = \frac{\int d\vec{r} D_W(\vec{r}) I(\vec{r})}{\int d\vec{r} I(\vec{r})}. \quad (S16)$$

The mean local directionality $\langle D_W \rangle$ for different cavities is presented in Fig. 3(b) of the main text.

## 3. Local structure size

The calculation of the local structure size of a lasing mode is illustrated by the example in Fig. S8. Figure S8(a) shows the intensity distribution of a high-Q mode of a Limaçon resonator with surface roughness, and Fig. S8(b) is a magnification of it in a local area (white box). It reveals a few standing waves in the direction $\sim 120^\circ$. The feature size in the direction transverse to propagation, often called the transverse wavelength, is obtained from the spatial intensity correlation function of mode $\mu$ around a point $\vec{r}_0 = (x_0, y_0)$,

$$C_{\mu}(\Delta\vec{r}; \vec{r}_0) = \int_{x_0-\delta x}^{x_0+\delta x} dx \int_{y_0-\delta y}^{y_0+\delta y} dy I(\vec{r}) I(\vec{r} + \Delta\vec{r}), \quad (S17)$$

where $I(\vec{r}) = |\Psi(\vec{r})|^2$ is the intensity distribution of the mode, $\Delta\vec{r} = (\Delta x, \Delta y)$, and the correlation function is normalized to $C_\mu(0) = 1$. We note that the integration is over a $2\delta x \times 2\delta y = 2\sigma \times 2\sigma$ large region around $\vec{r}_0$, where $2\sigma = 8\lambda/n = 1.9$ $\mu$m is the full width of the wavelet we used.

The spatial intensity correlation function is shown in Fig. S8(c). The central maximum of the function can
<figure><img src="image_14.png" /><figcaption>FIG. S8. Local structure size of a cavity resonance. (a) Calculated intensity distribution of a typical high-Q mode in a Limaçon cavity with rough boundary ($\lambda = 798.6$ nm, $Q = 1.05 \times 10^4$). (b) Magnification of a $2\sigma \times 2\sigma$ region [white box in (a)], where $2\sigma = 8\lambda/n = 1.9$ µm is the full width of the wavelet in Fig. S7(b). (c) Spatial intensity correlation function $C_\mu(\Delta\vec{r}; \vec{r}_0)$ in the magnified region. The blue dotted line indicates the half-maximum contour line where $C_\mu = 1/2$. (d) The half-maximum contour line (blue dotted line) is fitted by an ellipse (solid red line) with a major (minor) axis of $0.28$ µm ($0.06$ µm). The solid green line indicates the directionality diagram $|w_\mu(\vec{r}_0, \theta)|$ in a polar plot. Its maximum at $\theta = 123^\circ$ gives the local wave propagation direction, consistent with (b).</figcaption></figure>

have a circular or an elongated shape, depending on whether the local intensity structure is more isotropic or anisotropic. Its shape, size, and orientation are a good measure of the fine structure of the intensity distribution. We determine the structure size by extracting the half-maximum contour line where $C_\mu = 1/2$ around $\Delta\vec{r} = (0, 0)$ [blue dotted lines in Fig. S8(c)]. This contour line has an approximately elliptic shape, and it is plotted together with an ellipse fit (solid red line) in Fig. S8(d). The fitted ellipse has a major (minor) axis of $0.28$ µm ($0.06$ µm) oriented in the direction $\theta = 118^\circ$ ($28^\circ$). The transverse wavelength is given by the major diameter of the fitted ellipse, which is perpendicular to the local wave propagation direction as shown below.

The orientation of the ellipse is dictated by the local directionality, as shown by a polar plot of the directionality diagram $|w_\mu(\vec{r}_0, \theta)|$ (obtained from the wavelet transform) superimposed in Fig. S8(d) as a solid green line. Its maximum is at $\theta = 123^\circ$, which agrees well with

the orientation of the minor ellipse axis. Thus, it represents the local direction of wave propagation, and the major ellipse axis $s_1$ yields the transverse wavelength. If it is much larger than the in-medium wavelength, optical lensing effects may cause light to focus in the propagation direction, thus forming spatial filaments that can lead to lasing instabilities.

For every lasing mode $\mu$, we repeat the above procedure for every location $\vec{r}$ inside the resonator and obtain $\tilde{s}_\mu(\vec{r})$. The spatially-resolved structure size of a microlaser is obtained by averaging over all lasing modes, weighted by their power $P_\mu$ (obtained from the SPA-SALT simulation):

$$s(\vec{r}) = \frac{\sum_{\mu} P_{\mu} \tilde{s}_{\mu}(\vec{r})}{\sum_{\mu} P_{\mu}}, \qquad (S18)$$

Finally, we average the local structure size over the entire cavity area, again weighted by local intensity $I(\vec{r})$ [averaged intensity of all lasing modes, Eq. (S13)],

$$\langle s \rangle = \frac{\int d\vec{r} s(\vec{r}) I(\vec{r})}{\int d\vec{r} I(\vec{r})}, \qquad (S19)$$

The mean structure size $\langle s \rangle$ for each of the five cavity shapes is presented in Fig. 3(b) of the main text.

## E. Ray dynamics

### 1. Cavity geometry

In this section, we discuss in detail the ray dynamical properties of the five cavity geometries: D-shape, stadium, Limaçon, ellipse, and square.

The D-cavity is a circle of which a section has been cut off, where the distance from the circle center to the cut is $R_D/2$, where $R_D$ is the radius of the cavity. The stadium consists of a square with side length $a$ between two semi-circles of radius $a/2$. Both the D-cavity and stadium have completely chaotic ray dynamics [19, 20] with the exception of marginally unstable periodic orbits with measure zero in phase space. Their classical trajectories cover the whole billiard and are spatially extended, but due to the leakiness of the dielectric cavity, the long-lived trajectories feature central regions of reduced intensity [7].

The boundary of the Limaçon cavity [3] is defined in polar coordinates by $r(\varphi) = R_0(1 + \epsilon \cos \varphi)$, where $\varphi$ is the azimuthal angle, $R_0$ is the mean radius (averaged over $\varphi$), and $\epsilon = 0.42$ is the deformation parameter. In contrast to the D-cavity and stadium, its ray dynamics is predominantly, but not completely chaotic: since it is a convex cavity, there remain invariant surfaces with whispering-gallery type trajectories very close to its boundary [4, 21], and there is a small stable island with integrable ray dynamics around the horizontal diameter orbit [3, 4], as well as further tiny stable islands around other stable periodic orbits [22]. So technically it
is a mixed system, but we will consider it a chaotic billiard, since the integrable regions of phase space are too small to support any localized resonances for the cavity size considered here. The long-lived trajectories are concentrated at the cavity boundary and confined by total internal reflection. Even though the chaotic trajectories eventually explore the whole phase space, in an open dielectric cavity, they will quickly lose their intensity once they reach the central region of the cavity, as their incident angles onto the boundary are below the critical angle for total internal reflection.

Our elliptic cavities have the aspect ratio $b/a = 2$, where $b$ and $a$ represent the major and minor axes of the ellipse. Like the circle, the ellipse has completely integrable ray dynamics [23]. The long-lived trajectories are spatially localized near the boundary and are confined by total internal reflection. While there exists another type of modes based on so-called librator trajectories [24], these have a much shorter lifetime and thus cannot lase at the pump levels considered in this study.

Finally, the square has completely integrable ray dynamics like the ellipse, however, the structure of ray trajectories is quite different. Since the moduli of the momentum components, $|k_x|$ and $|k_y|$, are conserved quantities, there exist trajectories confined by total internal reflection at every reflection, which give rise to modes with extremely long lifetimes. The trajectories exhibit a very regular structure like for the ellipse, however, they cover the whole cavity homogeneously in contrast to the highly localized trajectories of the ellipse.

## 2. Ray tracing algorithm

We perform ray tracing for all five cavity shapes with smooth boundaries. We set the radius $R_D$ of the D-shaped cavity to unity, and the sizes of the other four cavity shapes are chosen such that they have the same cavity area ($\approx 2.53 \times R_D^2$). We consider the transverse-electric (TE) polarization of light, to be consistent with the polarization of edge emission from GaAs quantum-well lasers. The effective refractive index inside the cavity is $n_{\text{eff}} = 3.37$, and it is 1 outside.

We launch the optical rays with unit intensity and let them propagate with the in-medium speed of light, $c/n_{\text{eff}}$. Since the dielectric billiards we consider here are open systems, the ray intensity can change at reflections on the cavity boundary. Depending on the angle of incidence with respect to the boundary normal, $\chi \in (-\pi/2, \pi/2)$, the intensity either decreases due to refractive loss ($\sin|\chi| < 1/n_{\text{eff}}$) according to the Fresnel reflection coefficients or remains unchanged thanks to total internal reflection ($\sin|\chi| \ge 1/n_{\text{eff}}$), ignoring evanescent tunneling. We track the trajectories in space as well as their intensity decay as a function of the actual time of flight (instead of the number of reflections).

The aim of the ray tracing simulations is to calculate the local directionality analogously to the case of

the wave simulations discussed above in order to demonstrate the correspondence of rays and waves in this context. This requires averaging over specific ensembles of rays, which are different for chaotic and integrable billiards, respectively. The appropriate ensembles can be chosen via the initial conditions of the rays, as explained below.

For cavities with chaotic ray dynamics (D-cavity, stadium, and Limaçon), all high-Q resonances correspond to the same ensemble of long-lived trajectories [25, 26]. We launch an ensemble of $10^6$ trajectories at uniformly distributed random positions on the cavity boundary into uniformly distributed random initial directions $\sin\chi \in (-1, 1)$. In order to select the most long-lived trajectories, we sample the rays during the time interval in which the averaged ray intensity exhibits exponential decay after initial transients [7]. We choose the time intervals of [10.0, 30.0], [11.9, 35.7], and [43.0, 60.2] for D-cavity, stadium, and Limaçon, respectively, where the unit of time is $n_{\text{eff}}R_D/c$.

For integrable cavities (ellipse and square), each individual resonance corresponds to one torus in phase space [27]. We hence launch all rays at a fixed position on the cavity boundary with $10^3$ different angles of incidence $\chi$ such that $\sin\chi$ is evenly spaced in $(-1, 1)$. Then each ray propagates on a different torus in phase space and explores its torus entirely given enough time<fnref n="1" />. Here, we only sample the ray trajectories that are permanently confined by total internal reflection since the decaying trajectories correspond to low-Q modes that are irrelevant for lasing. The time interval for sampling rays is $[0, 10^4]n_{\text{eff}}R_D/c$, which is sufficiently long to homogeneously sample the tori.

A small modification of the initial conditions is needed for the square. Due to its symmetry with respect to the diagonal axes, a resonance of a square resonator exhibits eight plane-wave components with wave vectors ($k_x = \pm k_1, k_y = \pm k_2$) and ($k_x = \pm k_2, k_y = \pm k_1$) [16]. However, a single ray trajectory in a square billiard represents only four of these eight propagation directions. Therefore, we combine two trajectories with initial angles of incidence of $\chi$ and $\chi + \pi/2$ to represent a single resonance of the square cavity, where $\chi \in (-\pi/2, 0)$.

## 3. Local directionality

From all sampled ray trajectories, we sum the intensities of all rays with direction $\theta$ inside a local area centered at $\vec{r}$ to obtain $R(\vec{r}, \theta)$, where $\vec{r}$ is the position inside the cavity, and $\theta$ is the direction of ray propagation with respect to the horizontal axis. The bin size for sampling the direction $\theta$ is 1°. The radius of the local area is chosen as $(9.5 \times 10^{-2})R_D$, so that when $R_D$ is 10 µm as

<footnote n="1">¹ Except for periodic trajectories, which we exclude here.</footnote>
<figure><img src="image_15.png" /><figcaption>FIG. S9. Spatial intensity distributions of lasing modes. The top row is obtained from wave simulations for (a) D-cavity, (b) stadium, (c) Limaçon, (d) ellipse, and (e) square cavities with surface roughness. The field intensity is averaged over all lasing modes in three cavities with different surface roughness realizations that are statistically equivalent and described in subsection C 2. The scale bars are 5 µm long. The bottom row (f-j) represents the ray intensity distribution from ray tracing simulations with smooth cavity boundaries. The spatial structures inside the cavity show excellent agreement between wave and ray simulations.</figcaption></figure>

in the wave simulations, this radius becomes 0.95 µm, which is exactly the half-width of the wavelet used for wave simulations (see Fig. S7).

The local angular spread of ray propagation directions is defined as

$$ \sigma_R(\vec{r}) = \min_{\theta_0} \sqrt{\frac{\int_{-\frac{\pi}{2}}^{+\frac{\pi}{2}} \theta^2 [R(\vec{r}, \theta + \theta_0) + R(\vec{r}, \theta + \theta_0 - \pi)] d\theta}{\int_{-\pi}^{+\pi} R(\vec{r}, \theta) d\theta}}, \quad (S20) $$

which is analogous to Eq. (S14) for the wave simulations. The local directionality is defined by the inverse of the angular spread,

$$ D_R(\vec{r}) = \frac{1}{\sigma_R(\vec{r})}. \quad (S21) $$

The maps of $D_R(\vec{r})$ are presented in Figs. 4(e1-e5) of the main text, showing an excellent agreement with the local directionality maps $D_W(\vec{r})$ obtained from wave simulations.

Finally, we calculate the mean directionality, averaged over the entire cavity area

$$ \langle D_R \rangle = \frac{\int d\vec{r} D_R(\vec{r}) I(\vec{r})}{\int d\vec{r} I(\vec{r})}, \quad (S22) $$

where $I(\vec{r}) = \int R(\vec{r}, \theta) d\theta$ is the ray intensity distribution inside the cavity. The calculated $\langle D_R \rangle$ for the five different cavity shapes, as well as the comparison with $\langle D_W \rangle$, are presented in Fig. 3(b) of the main text.

### 4. Spatial intensity distributions

Finally, we confirm ray-wave correspondence by comparing the results of ray tracing with wave simulations. To this end, we present the averaged intensity distributions inside the cavity in Fig. S9.

For wave simulations of microcavities with rough boundary (top row), the lasing modes at the pumping level of 10 times the lasing threshold are calculated using SPA-SALT, and the average intensity distributions are given by $I(\vec{r}) = \sum_{\mu} P_{\mu} |\Psi_{\mu}(\vec{r})|^2$, where the lasing modes are weighted by their power $P_{\mu}$. We also averaged the intensity distributions over three surface roughness realizations.

For ray simulations (bottom row), we show $I(\vec{r}) = \int R(\vec{r}, \theta) d\theta$. Here we use a finer spatial resolution of $(1.2 \times 10^{-3}) R_D$, in order to reveal detailed features in the intensity profiles inside a cavity.

The ray and wave simulations yield qualitatively consistent intensity profiles. The D-cavity and stadium have a quite homogeneous average intensity distribution, but we observe central regions of lower intensity due to the leakage at the dielectric cavity boundaries [7]. The average intensity distributions of the Limaçon and ellipse resonators are localized at the boundary due to the prevalence of whispering gallery modes. Thus they are highly localized, though the boundary roughness in the wave simulations increases the area they cover compared to ray simulations. The ray intensity distribution in the square with a smooth boundary is perfectly homogeneous over the whole cavity. However, the rough boundary in wave simulation leads to the formation of diagonal stripes of
higher intensity, rendering the wave intensity distribution a bit less homogeneous compared to the ray intensity. In general, the ray tracing simulations in 2D cavities with smooth boundaries correctly predict all the main characteristics of the spatial structure of lasing modes in wave simulations of cavities with small boundary roughness.

Considering the huge amount of computational resources required for wave simulations, ray tracing can be an efficient tool to predict the lasing characteristics of realistic microcavities.

* These authors contributed equally.

† hui.cao@yale.edu

[1] A. Cerjan, S. Bittner, M. Constantin, M. Guy, Y. Zeng, Q. J. Wang, H. Cao, and A. D. Stone, Multimode lasing in wave-chaotic semiconductor microlasers, Phys. Rev. A **100**, 063814 (2019).

[2] K. Kim, S. Bittner, Y. Jin, Y. Zeng, S. Guazzotti, O. Hess, Q. J. Wang, and H. Cao, Sensitive control of broad-area semiconductor lasers by cavity shape, APL Photonics **7**, 056106 (2022).

[3] M. Robnik, Classical dynamics of a family of billiards with analytic boundaries, J. Phys. A **16**, 3971 (1983).

[4] J. Wiersig and M. Hentschel, Combining directional light output and ultralow loss in deformed microdisks, Phys. Rev. Lett. **100**, 033901 (2008).

[5] K. Kim, S. Bittner, Y. Jin, Y. Zeng, Q. Wang, and H. Cao, Spatiotemporal lasing dynamics in a Limaçon-shaped microcavity, Opt. Lett. **48**, 574 (2023).

[6] S. Bittner, S. Guazzotti, Y. Zeng, X. Hu, H. Yılmaz, K. Kim, S. S. Oh, Q. J. Wang, O. Hess, and H. Cao, Suppressing spatio-temporal lasing instabilities with wave-chaotic microcavities, Science **361**, 1225 (2018).

[7] S. Bittner, K. Kim, Y. Zeng, Q. J. Wang, and H. Cao, Spatial structure of lasing modes in wave-chaotic semiconductor microcavities, N. J. Phys. **22**, 083002 (2020).

[8] J. Ohtsubo, *Semiconductor Lasers - Stability, Instability and Chaos*, 3rd ed. (Springer, Heidelberg, 2013).

[9] C.-G. Ma, J.-L. Xiao, Z.-X. Xiao, Y.-D. Yang, and Y.-Z. Huang, Chaotic microlasers caused by internal mode interaction for random number generation, Light: Sci. Appl. **11**, 187 (2022).

[10] K. Kim, S. Bittner, Y. Zeng, S. Guazzotti, O. Hess, Q. J. Wang, and H. Cao, Massively parallel ultrafast random bit generation with a chip-scale laser, Science **371**, 948 (2021).

[11] H. E. Türeci, H. G. L. Schwefel, P. Jacquod, and A. D. Stone, Modes of wave-chaotic dielectric resonators, Prog. Opt. **47**, 75 (2005).

[12] H. Cao and J. Wiersig, Dielectric microcavities: Model systems for wave chaos and non-Hermitian physics, Rev. Mod. Phys. **87**, 61 (2015).

[13] S. F. Liew, L. Ge, B. Redding, G. S. Solomon, and H. Cao, Pump-controlled modal interactions in microdisk lasers, Phys. Rev. A **91**, 043828 (2015).

[14] L. Ge, Y. D. Chong, and A. D. Stone, Steady-state *ab initio* laser theory: Generalizations and analytic results, Phys. Rev. A **82**, 063824 (2010).

[15] A. Cerjan, B. Redding, L. Ge, S. F. Liew, H. Cao, and A. D. Stone, Controlling mode competition by tailoring the spatial pump distribution in a laser: a resonance-based approach, Opt. Express **24**, 26006 (2016).

[16] S. Bittner, E. Bogomolny, B. Dietz, M. Miski-Oglu, and A. Richter, Experimental observation of localized modes in a dielectric square resonator, Phys. Rev. E **88**, 062906 (2013).

[17] S. Bittner, E. Bogomolny, B. Dietz, M. Miski-Oglu, and A. Richter, Dielectric square resonator investigated with microwave experiments, Phys. Rev. E **90**, 052909 (2014).

[18] I. Daubechies, *Ten Lectures on Wavelets* (Society for Industrial and Applied Mathematics, 1992).

[19] L. A. Bunimovich, On the ergodic properties of nowhere dispersing billiards, Commun. Math. Phys. **65**, 295 (1979).

[20] S. Ree and L. E. Reichl, Classical and quantum chaos in a circular billiard with a straight cut, Phys. Rev. E **60**, 1607 (1999).

[21] V. F. Lazutkin, The existence of caustics for a billiard problem in a convex domain, Math. USSR Izv. **7**, 185 (1973).

[22] H. R. Dullin and A. Bäcker, About ergodicity in the family of limaçon billiards, Nonlinearity **14**, 1673 (2001).

[23] M. V. Berry, Regularity and chaos in classical mechanics, illustrated by three deformations of a circular 'billiard', Eur. J. Phys. **2**, 91 (1981).

[24] J. B. Keller and S. I. Rubinow, Asymptotic solution of eigenvalue problems, Ann. Phys. **9**, 24 (1960).

[25] E. G. Altmann, J. S. E. Portela, and T. Tél, Leaking chaotic systems, Rev. Mod. Phys. **85**, 869 (2013).

[26] R. Ketzmerick, K. Clauß, F. Fritzsch, and A. Bäcker, Chaotic resonance modes in dielectric cavities: Product of conditionally invariant measure and universal fluctuations, Phys. Rev. Lett. **129**, 193901 (2022).

[27] M. V. Berry, Regular and irregular semiclassical wavefunctions, J. Phys. A **10**, 2083 (1977).