††thanks: These authors contributed equally††thanks: These authors contributed
equally††thanks: These authors contributed equally

# Direct Observation of Landau Levels in Silicon Photonic Crystals

Maria Barsukova  Department of Physics, The Pennsylvania State University,
University Park, PA, USA  Fabien Grisé  Department Astronomy & Astrophysics,
The Pennsylvania State University, University Park, PA, USA  Zeyu Zhang
Department of Physics, The Pennsylvania State University, University Park, PA,
USA  Sachin Vaidya  Department of Physics, The Pennsylvania State University,
University Park, PA, USA  Department of Physics, Massachusetts Institute of
Technology, Cambridge, MA, USA  Jonathan Guglielmon  Department of Physics,
The Pennsylvania State University, University Park, PA, USA  Michael I.
Weinstein  Department of Applied Physics and Applied Mathematics and
Department of Mathematics, Columbia University, New York, NY, USA  Li He
Department of Physics and Astronomy, University of Pennsylvania, Philadelphia,
PA, USA  Bo Zhen  Department of Physics and Astronomy, University of
Pennsylvania, Philadelphia, PA, USA  Randall McEntaffer  Department Astronomy
& Astrophysics, The Pennsylvania State University, University Park, PA, USA
Mikael C. Rechtsman  Department of Physics, The Pennsylvania State University,
University Park, PA, USA<EMAIL_ADDRESS>
###### Abstract

We experimentally observe photonic Landau levels that arise due to a strain-
induced pseudomagnetic field in a silicon photonic crystal slab. The Landau
levels are dispersive (i.e., they are not flat bands) due to the distortion of
the unit cell by the strain. We employ an additional strain which induces a
pseudoelectric potential to flatten them.

When electrons are confined to a two-dimensional plane and are subjected to an
out-of-plane magnetic field, they move in circular cyclotron orbits as a
result of the Lorentz force. In the quantum domain, this cyclotron motion is
quantized, and as a consequence, the electrons’ energy spectrum splits into
discrete, highly degenerate states called Landau levels. The integer and
fractional quantum Hall effects [1, 2] arise as a direct result; in the
fractional case, it is the high degeneracy of Landau levels (i.e., that they
are flat bands) that gives rise to effectively strong electron-electron
interactions and leads to the fractionalization of charge.

In free space, photons do not respond to external magnetic fields because they
do not carry charge; yet, when propagating in magneto-optical materials, they
may respond indirectly as a result of the material’s magnetic response.
However, this response is weak at optical frequencies. In 2012, an approach
was put forward for emulating magnetic behavior in photonic systems by
inhomogeneously straining a photonic lattice. [3]. This implementation was
based on an idea proposed for electrons in graphene, where a strain pattern
imposed on the lattice would introduce an effective gauge field at the Dirac
point, causing electrons to behave as though there were a strong field
present, even in the absence of a real magnetic field [4]. The effect was
later demonstrated by directly observing Landau levels in graphene bubbles,
where a strain corresponding to an enormous ‘pseudomagnetic’ field of $300T$
was imposed [5]. Since the original photonic experiment, Landau levels were
also proposed and observed in exciton-polariton condensates [6] and in
mechanical systems [7, 8, 9]. Moreover, there have been a number of
theoretical proposals for how Landau levels may be used in the context of
photonics that are intrinsically distinct from the electronic case [10, 11].

Here, we directly observe Landau levels in two-dimensional silicon photonic
crystal slabs in the nanophotonic domain. Moreover, we go beyond purely
pseudomagnetic effects and demonstrate that strains corresponding to
pseudoelectric fields act to flatten the Landau levels that inherit dispersion
from the form of the pseudomagnetic strain. There are several key differences
and advantages of pseudomagnetism in photonic crystals compared to previous
realizations of photonic pseudomagnetism. First, photonic crystals have been
demonstrated to enhance light-matter interaction via cavity modes and flat
bands [12, 13, 14]. This enhancement is generated as a result of the lattice.
In contrast, for systems composed of individual, isolated guiding or resonant
elements (as in Refs. [3, 6]), lattice effects are not leveraged because
strong enhancement would occur even in a single site. Second, besides having
unit cells that are an order of magnitude smaller, photonic crystals can in
practice have much larger system sizes compared to previous realizations
(millions compared to hundreds of unit cells), and can be realized with
smaller loss in the silicon platform. Since Landau level degeneracy scales
with system size and the linewidth increases with loss, photonic crystals
allow for increased degeneracy and significantly improved spectral resolution
of the levels.

Further, since photonic crystals do not have an associated tight-binding
theory, the original theoretical framework relating strain to pseudomagnetism
is not directly applicable, rendering a new understanding necessary; the
appropriate effective Hamiltonians for strain-dependent emergent parameters
for two-dimensional photonic crystals were derived in our previous theoretical
work [15], and are extended to the slab geometry here (2D slab embedded in 3D
space). Our establishment of a new analytical method of understanding and
describing aperiodicity in photonic crystals (i.e., using pseudomagnetic
fields) will be useful in their optimization for many different functions;
this has traditionally been approached by using direct numerical optimization
[16, 17, 18].

Our starting point is a photonic crystal structure consisting of rounded
triangular air holes in a silicon slab [19] that rests on a silica substrate.
The holes form an underlying honeycomb pattern with $C_{6v}$-symmetry. As a
result, this lattice hosts Dirac points at the $\mathbf{K}$ and
$\mathbf{K^{\prime}}$ points in the Brillouin zone [20, 21, 22]. As these
Dirac points lie below the light line of vacuum, they are not detectable via
free-space excitation. To allow radiative coupling from outside the slab, we
introduce a small period-doubling perturbation by changing the size of some of
the holes (more details can be found in Supplementary Information Section 2).
This makes the unit cell of the lattice rectangular, and the band structure is
folded such that the Dirac cone resides along the $k_{x}$ axis and lies above
the light line of vacuum. A scanning electron microscope image of the
structure is shown in Fig. 1(a); the period-doubled unit cell is shaded in
purple.

Figure 1: Unstrained and strained photonic crystals and their calculated band
structures. (a) Scanning electron microscope image of the unstrained periodic
photonic crystal. The period-doubled unit cell is highlighted in purple; $a$,
lattice constant without period doubling, is indicated in white. (b) Simulated
band structure for the unstrained structure, showing a Dirac point for the TE-
like modes. Color indicates quality factors of the states. (c) SEM image of
the strained pattern, in which the strain causes the loss of periodicity in
the $x$ direction. The $y$-periodic unit cell of the strained pattern is
highlighted in purple. (d) Simulated band structure for the strained pattern
($\kappa=0.0632a^{-1}$), showing the emergence of Landau levels in the
vicinity of the Dirac frequency. Only the first eleven Landau levels
$(n=-5\text{ to }5)$ are shown.

We numerically compute the band structure (in the transverse electric
polarization) using the guided mode expansion method as implemented in the
open-source software package Legume [23]. Fig. 1(b) shows the linearly-
dispersing transverse electric (TE)-like bands that exhibit a Dirac point at
$\mathbf{K}$, with a frequency $\omega_{\rm D}=0.318\ [2\pi ca^{-1}]$. Here
$a$ is the lattice constant of the underlying hexagonal lattice structure and
$c$ is the speed of light. The period-doubling procedure very slightly changes
the Dirac frequency (see Supplementary Information Section 2).

Next, we introduce a strain pattern in our structure by deforming the lattice
as shown in Fig. 1(c). Here, the term strain  refers not to a strain induced
by a physically applied stress, but to the deformation of the dielectric
pattern that is directly etched into the silicon. The specific strain pattern
is achieved by mapping every point $(x,y,z)$ to $(x,y+a(\kappa x)^{2},z)$,
where $\kappa$ is the strength of the strain. This deformation breaks
periodicity in the $x$ direction, but retains periodicity along the $y$
direction. The spatial scale separation ensured by the assumption of small and
slowly varying strain, $\kappa a\ll 1$, allows us to develop a multiple scale
[24] variant of degenerate perturbation theory to expand the eigenstates and
eigenvalues of the strained system. The eigenstates are, to leading order in
$\kappa$, a slow spatial modulation of the degenerate Bloch modes associated
with the Dirac point of the unstrained ($\kappa=0$) structure.

The resulting effective Hamiltonian, which incorporates the strain, is given
by

$\mathcal{H}_{\rm eff}=E_{\rm D}\sigma_{0}+v_{\rm
D}\left[\left(-i\frac{\partial}{\partial
x}\right)\sigma_{1}+\left(-i\frac{\partial}{\partial
y}+\frac{4ab_{*}\kappa^{2}}{v_{\rm D}}x\right)\sigma_{2}\right],$ (1)

where $E_{\rm D}=\left(\omega_{\rm D}/c\right)^{2}$, $\sigma_{0}$,
$\sigma_{1}$ and $\sigma_{2}$ are Pauli matrices, and $v_{\rm D}=0.915a^{-1}$
and $b_{*}=0.606a^{-2}$ are two parameters calculated from the modes of the
unstrained structure at energy $E_{\rm D}$. A detailed derivation can be found
in Supplementary Information section 3, where explicit expressions for $b_{*}$
and $v_{\rm D}$ in terms of the eigenstates of the periodic structure are
displayed. We note that the effective Hamiltonian displayed in Eq. (1) is
derived directly from the continuum theory of photonic crystals; this is
fundamentally different from the previous work [4] based on the tight-binding
approximation. Our approach extends the the methods of Ref. [15] to the three-
dimensional setting of the slab geometry, where vectorial effects play a role.

Equation (1) corresponds to a two-dimensional Dirac Hamiltonian describing
massless spin-$1/2$ relativistic particles under a constant (pseudo)magnetic
field pointing in the out-of-plane direction, where the magnetic field has a
strength of $B_{\rm eff}=4ab_{*}\kappa^{2}/v_{\rm D}$ and is described by a
vector potential in the Landau gauge. The discrete energies that are
eigenvalues of the Hamiltonian in Eq. (1) for an electron are known as Landau
levels. The energy eigenvalue of the $n^{th}$ level is proportional to
$\sqrt{|n|}$, where $n$ is an integer. Analogously, for our photonic crystal
slabs, the frequency eigenvalues of the electromagnetic eigenmodes are, to
first order in $\kappa$, proportional to $\sqrt{|n|}$ and can be expressed as
$\omega_{n}=\omega_{\rm D}\pm(c^{2}v_{\rm D}/\sqrt{2}\omega_{\rm
D})\sqrt{B_{\rm eff}\lvert n\rvert}$, where $n$ is an integer.

Figure 2: Observation of Landau levels in the spectrum of the strained
photonic crystal. (a) Experimentally measured band structure of the unstrained
honeycomb pattern, showing the Dirac point in the vicinity of $\omega_{\rm
D}=0.32~{}[2\pi ca^{-1}]$. (b) Measured band structure for the strained
pattern. The uniform pseudomagnetic field created by the strain causes the
Dirac point to split into sets of discrete Landau levels. Here
$\kappa=0.0632a^{-1}$, as in the numerical results shown in Fig. 1(d). (c)
Landau level energy spacing
($\Delta\omega_{n,0}=\omega_{n}-\omega_{0}^{\prime}$) is a linear function of
the strain strength $\kappa$. (d) The Landau levels lie at energies
proportional to $\sqrt{\lvert n\rvert}$; this behavior is captured by plotting
the dimensionless quantity
$\frac{\lvert\omega_{n}-\omega_{0}^{\prime}\rvert}{\sqrt{\lvert
n\rvert}\kappa}$ for each Landau level, including six different samples of
different strain strengths, $\kappa$. In both (c) and (d), the solid lines are
the theoretical prediction with no free parameters.

To corroborate our analytical results given in Eq. (1), we also perform
numerical simulations of the strained structure using the guided-mode
expansion method. The strain is implemented in a dielectric profile which
spans 199 period-doubled unit cells in the $x$-direction. Due to the
preservation of lattice periodicity along the $y$-direction, $k_{y}$ is
conserved and the frequencies of the bands can be plotted as functions of
$k_{y}$, as shown in Fig. 1(d). Here, we observe the splitting of the spectrum
near the Dirac point into discrete Landau levels due to the strain-induced
pseudomagnetic field, where the spacing of these levels is proportional to
$\sqrt{|n|}$ for a fixed value of $\kappa$.

To demonstrate the formation of Landau levels in such a system, we use
electron-beam lithography to fabricate both the periodic and the strained
patterns in a silicon slab ($\varepsilon=12.11$) on top of a silica substrate
($\varepsilon=2.25$). A detailed description of the fabrication methods can be
found in Supplementary Material Section 1. Figures 1(a) and (c) show scanning
electron microscope (SEM) images of the fabricated structures. The structure
in Fig. 1(a) has a periodicity along the $x$ direction of $2a=980~{}nm$.

To experimentally characterize the photonic bands of these structures, we
perform angle- and frequency-resolved reflection measurements. The samples are
illuminated by a tunable continuous wave laser (Keysight 81606A) with a
wavelength range of $\lambda=1.45$ \- $1.65~{}\mu m$ ($\pm 1.5~{}pm$ absolute
wavelength resolution accuracy), and a laser linewidth coherence control of 10
kHz. We measure the iso-frequency contours of the fabricated photonic crystal
slabs using back focal plane (BFP) imaging. We then extract the Landau-level
band structures by observing the photonic crystal resonances at a fixed
$k_{x}$ corresponding to the location of the Dirac point of the unstrained
structure. Details of the experimental setup can be found in Supplementary
Material Section 1.

Figure 3: Landau levels are excited at increasing $k_{y}$ in the band
structure as the input beam is moved from left to right through the sample.
Experimental band structures showing the excitation of the Landau levels when
the beam is injected at (a) $x=-20~{}\mu$m, (b) $x=0~{}\mu$m, and (c)
$x=20~{}\mu$m. (d) Center of Landau level excitation in $k_{y}$ plotted
against the position of the input beam, $x$. The linear relationship is
characteristic of Landau level eigenstates in the Landau gauge. The solid line
represents the theoretical prediction, with no free parameters.

Fig. 2(a) shows the bands of the unstrained structure, obtained by BFP
imaging, where we clearly observe linearly dispersing bands near the Dirac
point. We note that a small gap is observed at the Dirac point - this is due
to inevitable fabrication disorder that breaks inversion symmetry. We show in
Supplementary Information Section 4 that the breaking of inversion symmetry
affects the zeroth Landau level significantly more than the others.

Next, we measure the bands of the strained photonic crystal slabs described
above and find the emergence of discrete Landau levels, as shown in Fig. 2(b).
While the effective theory predicts that the Landau levels should be flat, we
see that they are dispersive in both simulation (Fig. 1(d)) and experiment
(Fig. 2(b)), i.e., the bands are concave-up. This arises due to the fact that,
by adding strain, the unit cell is distorted locally as a function of $x$.
This distortion effectively adds a parabolic potential to the Hamiltonian
(i.e., $H\sim x^{2}\sigma_{0}$), which in turn causes the dispersion of the
Landau level bands. A detailed explanation can be found in Supplementary
Information Section 5.

According to the effective theory (Eq. 1), the $n=0$ level should be at the
center of all Landau levels. However, due to the aforementioned inversion-
symmetry breaking, this level is slightly shifted away from the center (see
Supplementary Information Section 4). As a result, we use a new reference
frequency of $\omega_{0}^{\prime}=\frac{1}{2}(\omega_{-1}+\omega_{1})$ as the
Dirac frequency to calculate the Landau level spacings, defined as
$\omega_{n}-\omega_{0}^{\prime}$. In Fig. 2(c), we compare the theoretically
and experimentally obtained level spacings at $k_{y}=0[2\pi a^{-1}]$ under
different strain strengths (characterized by $\kappa$) and observe good
agreement between the two. From the experimental data, we also calculate the
normalized quantity $|\omega_{n}-\omega_{0}^{\prime}|/(\kappa\sqrt{|n|})$,
which should be a constant for all Landau levels. We again observe good
agreement between experiment and the theoretically-predicted value of
$0.0823[2\pi c]$, as shown in Fig. 2(d). In both Figs. 2(c) and (d), the
theoretical plots (solid lines) are obtained directly from analytical
predictions, and have no free parameters.

Figure 4: The introduction of a pseudoelectric potential acts to flatten the
Landau level bands. (a) Schematic figure depicting the strain profile that
corresponds to the pseudoelectric potential used to mitigate the Landau level
dispersion. (b) Simulated band structure with $\kappa=0.0632a^{-1}$ and
$\beta=0.0364$, where the Landau level $n=0$ is predicted to be flattened. (c)
Experimental data for $\kappa=0.0632a^{-1}$ and $\beta=0.0364$, showing the
flattening of the Landau levels.

It is clear from Fig. 2(b) that, as $n$ decreases, the range of $k_{y}$ values
over which the $n^{th}$ Landau level is observed becomes smaller. This can be
intuitively understood as arising from the interaction of the Landau level
states with other states that reside toward the far left and right sides of
the sample. These states rise in energy as one moves away from the sample
center along the $x$ direction. We know from Eq. (1) that the Landau level
states are harmonic oscillator eigenstates centered at $x=k_{y}/B_{\rm eff}$,
with spatial widths of $\Delta x_{n}=\sqrt{(2|n|+\delta_{0,n})/2B_{\rm eff}}$.
As $k_{y}$ is increased, the Landau level center translates and the tail of
the Landau level eventually interacts with the states mentioned above, leading
to an increased linewidth. More details are given in Supplementary Information
Section 3.

The fact that the $x$-position of the Landau level state varies linearly with
$k_{y}$ leads to another clear observable: when the input beam is moved from
left to right in real space along the $x$-direction, the Landau level states
at increasing $k_{y}$ are selectively excited, and therefore appear more
clearly in the band structure. We observe this effect directly, as shown in
Fig. 3(a) through (c): when the input beam is on the left side of the sample
(i.e., $x<0$), we see that the modes on the left side of the band structure
($k_{y}<0$) are more strongly excited, but as the input beam is moved
rightward, we observe that the modes on the right side of the band structure
are increasingly excited.

To further study the relationship between $x$ and $k_{y}$, we extract the
boundary in $k_{y}$-space between the modes that are excited and those that
are not excited. For input beams positioned left of center, we extract the
right boundary, and for input beams positioned right of center, we extract the
left boundary. The boundary values differ from the excitation centers by an
overall offset, which we remove by fitting the data to a line and subtracting
the intercept (one for the left-boundary data and one for the right-boundary
data). Using this procedure, we obtain the relationship between the Landau
level horizontal position and the average vertical momentum, $k_{y}$, of the
excited modes. The linear relationship between these, as shown in Fig. 3(d),
evidences the direct proportionality between the Landau level positions and
$k_{y}$.

We next turn our attention to the mitigation of the Landau level dispersion.
As explained earlier, Eq. 1 predicts flat Landau levels. However, in the
simulations and experiments, the Landau levels exhibit quadratic dispersion as
$k_{y}$ is varied. As shown in Ref. [15], it is possible to mitigate this
dispersion by introducing an additional strain profile, which induces a
pseudo-electric potential. Specifically, we add a cubic term to the
deformation such that the point $(x,y,z)$ is mapped to $(x+a\beta(\kappa
x)^{3},y+a(\kappa x)^{2},z)$. The parameter $\beta$ controls the strength of
this additional strain in the $x$-direction. A schematic of the strained
structure, which induces both pseudomagnetic and pseudoelectric fields, is
shown in Fig. 4(a) (further details are given in Supplementary Information
Section 5). The reason why the pseudoelectric field counters the Landau level
dispersion, to leading order, can be explained as follows. To leading order,
the form of the pseudoelectric field gives rise to a potential $V_{\rm
eff}=3a\beta m\kappa^{2}x^{2}\sigma_{0}$ (to be added to (1)) which is similar
to that which creates the dispersion in the first place (here $m=-3.28a^{-2}$
is a parameter calculated entirely from the states of the periodic structure).
Since the spatial positions of the Landau level eigenstates grow linearly with
$k_{y}$, a quadratic potential in $x$ is equivalent to a parabolic dispersion
in $k_{y}$. An appropriate choice of the field strength (and sign) will then
counteract the original dispersion induced by the strain associated with the
pseudomagnetic field.

By choosing $\beta$ appropriately, the quadratic dispersion of the Landau
levels can be mitigated, leading to nearly flat bands. We note that each
Landau level requires a different value of $\beta$ to counteract its
dispersion. More details are given in Supplementary Information Section 5.
Fig. 4(b) shows numerical simulations of the flattened Landau levels for a
structure with pseudomagnetic and pseudoelectric fields induced by a strain
with $\kappa=0.0632a^{-1}$ and $\beta=0.0364$. Here, the $n=0$ level is
targeted, but other levels are also evidently flatter. Fig. 4 (c) shows the
experimental data for a strained structure with the same values of $\kappa$
and $\beta$ given above, where a good agreement is observed between theory and
experiment.

In conclusion, we have directly observed Landau levels in the spectra of two-
dimensional silicon photonic crystal slabs. As in graphene, the Landau level
energies are proportional to $\sqrt{|n|}$, where $n$ is an integer. The Landau
level bands are found to be dispersive, which can be explained by a distortion
of the unit cell as a result of the strain. We further showed that this
dispersion can be mitigated by adding an additional strain that induces a
position-dependent pseudoelectric field (i.e., a potential). Landau levels
constitute a new methodology for enhancing light-matter interaction which is
distinct from standard slow light or cavity enhancement, because a flat band
acts essentially as a ‘cavity everywhere in space’. The realization of optical
pseudomagnetism prompts several new questions and directions of inquiry,
including: whether Landau-level flat bands can be used to enhance light-matter
coupling more efficiently than conventional photonic crystal flat bands or
other points of high degeneracy (such as van Hove singularities); the question
of the nature of wave mixing processes such as four-wave mixing among Landau
levels; and whether the square-root structure of the eigenvalue spacing can
lead to different properties associated with entangled pair or frequency comb
generation. More broadly, the framework of pseudomagnetism gives an analytical
handle on aperiodic photonic structures, allowing for a new approach to
designing devices and better understanding their behavior.

## Acknowledgements

We gratefully acknowledge funding support from the Office of Naval Research
MURI program under agreement number N00014-20-1-2325, the Air Force Office of
Scientific Research MURI program under agreement number FA9550-22-1-0339, as
well as the Kaufman and Packard foundations under grant numbers KA2020-114794
and 2017-66821, respectively. This research was also supported in part by
National Science Foundation grants DMS-1620422 (MCR), DMS-1620418 (MIW),
DMS-1908657 (MIW) and DMS-1937254 (MIW), as well as Simons Foundation Math + X
Investigator Award #376319 (MIW). The authors acknowledge the Nanofabrication
Lab within the Materials Research Institute at Penn State and the help of
Michael Labella, as well as seed funding from the Center for Nanofabricated
Optics at Penn State University. F.G. thanks GenISys and, in particular, Roger
McCay for his help in optimizing the fracturing of the electron-beam patterns.
M.B. thanks Sebabrata Mukherjee and Alexander Cerjan for fruitful discussions
in the early stages of the project and help with numerical optimization.

We would like to note that the group of Ewold Verhagen has concurrently posted
a similar work on the observation of Landau levels in photonic crystals.

## References

  * Tsui  _et al._ [1982] D. C. Tsui, H. L. Stormer, and A. C. Gossard, Physical Review Letters 48, 1559 (1982).
  * Klitzing  _et al._ [1980] K. v. Klitzing, G. Dorda, and M. Pepper, Physical review letters 45, 494 (1980).
  * Rechtsman  _et al._ [2013] M. C. Rechtsman, J. M. Zeuner, A. Tünnermann, S. Nolte, M. Segev, and A. Szameit, Nature Photonics 7, 153 (2013).
  * Guinea  _et al._ [2010] F. Guinea, M. I. Katsnelson, and A. Geim, Nature Physics 6, 30 (2010).
  * Levy  _et al._ [2010] N. Levy, S. Burke, K. Meaker, M. Panlasigui, A. Zettl, F. Guinea, A. C. Neto, and M. F. Crommie, Science 329, 544 (2010).
  * Jamadi  _et al._ [2020] O. Jamadi, E. Rozas, G. Salerno, M. Milićević, T. Ozawa, I. Sagnes, A. Lemaître, L. Le Gratiet, A. Harouri, I. Carusotto, _et al._ , Light: Science & Applications 9, 144 (2020).
  * Abbaszadeh  _et al._ [2017] H. Abbaszadeh, A. Souslov, J. Paulose, H. Schomerus, and V. Vitelli, Physical review letters 119, 195502 (2017).
  * Brendel  _et al._ [2017] C. Brendel, V. Peano, O. J. Painter, and F. Marquardt, Proceedings of the National Academy of Sciences 114, E3390 (2017).
  * Peri  _et al._ [2019] V. Peri, M. Serra-Garcia, R. Ilan, and S. D. Huber, Nature Physics 15, 357 (2019).
  * Schomerus and Halpern [2013] H. Schomerus and N. Y. Halpern, Physical review letters 110, 013903 (2013).
  * Lledó  _et al._ [2022] C. Lledó, I. Carusotto, and M. Szymanska, SciPost Physics 12, 068 (2022).
  * Li  _et al._ [2008] J. Li, T. P. White, L. O’Faolain, A. Gomez-Iglesias, and T. F. Krauss, Optics express 16, 6227 (2008).
  * Jukić  _et al._ [2012] D. Jukić, H. Buljan, D.-H. Lee, J. D. Joannopoulos, and M. Soljačić, Optics Letters 37, 5262 (2012).
  * Yang  _et al._ [2023] Y. Yang, C. Roques-Carmes, S. E. Kooi, H. Tang, J. Beroz, E. Mazur, I. Kaminer, J. D. Joannopoulos, and M. Soljačić, Nature 613, 42 (2023).
  * Guglielmon  _et al._ [2021] J. Guglielmon, M. C. Rechtsman, and M. I. Weinstein, Physical Review A 103, 013505 (2021).
  * Men  _et al._ [2014] H. Men, K. Y. Lee, R. M. Freund, J. Peraire, and S. G. Johnson, Optics express 22, 22632 (2014).
  * Vuckovic  _et al._ [2002] J. Vuckovic, M. Loncar, H. Mabuchi, and A. Scherer, IEEE Journal of Quantum Electronics 38, 850 (2002).
  * Bendsoe and Sigmund [2003] M. P. Bendsoe and O. Sigmund,  _Topology optimization: theory, methods, and applications_ (Springer Science & Business Media, 2003).
  * Barik  _et al._ [2016] S. Barik, H. Miyake, W. DeGottardi, E. Waks, and M. Hafezi, New Journal of Physics 18, 113013 (2016).
  * Neto  _et al._ [2009] A. C. Neto, F. Guinea, N. M. Peres, K. S. Novoselov, and A. K. Geim, Reviews of modern physics 81, 109 (2009).
  * Fefferman and Weinstein [2012] C. L. Fefferman and M. I. Weinstein, Journal of the American Mathematical Society 25, 1169 (2012).
  * Lee-Thorp  _et al._ [2019] J. P. Lee-Thorp, M. I. Weinstein, and Y. Zhu, Archive for Rational Mechanics and Analysis 232, 1 (2019).
  * Minkov  _et al._ [2020] M. Minkov, I. A. Williamson, L. C. Andreani, D. Gerace, B. Lou, A. Y. Song, T. W. Hughes, and S. Fan, Acs Photonics 7, 1729 (2020).
  * Bender and Orszag [1999] C. M. Bender and S. A. Orszag,  _Advanced mathematical methods for scientists and engineers I: Asymptotic methods and perturbation theory_ , Vol. 1 (Springer Science & Business Media, 1999).