# Chiral TeraHertz surface plasmonics

Ian Aupiais1, Romain Grasset1, Dmitri Daineka2, Javier Briatico3,
Luca Perfetti1, Jean-Paul Hugonin4, Jean-Jacques Greffet4, Yannis Laplace1

(1 LSI, CEA/DRF/IRAMIS, CNRS, Ecole Polytechnique, Institut Polytechnique de
Paris, Palaiseau, France
2 LPICM, CNRS, Ecole Polytechnique, Institut Polytechnique de Paris,
Palaiseau, France
3 Unité Mixte de Physique, CNRS, Thales, Université Paris Saclay, 91767
Palaiseau, France
4 Université Paris-Saclay, Institut d’Optique Graduate School, CNRS,
Laboratoire Charles Fabry, 91127 Palaiseau, France)

###### Abstract

Chiral engineering of TeraHertz (THz) light fields and the use of the
handedness of light in THz light-matter interactions promise many novel
opportunities for advanced sensing and control of matter in this frequency
range. Unlike previously explored methods, this is achieved here by leveraging
the chiral properties of highly confined THz surface plasmon modes. More
specifically, we design ultrasmall surface plasmonic-based THz cavities and
THz metasurfaces that display significant and adjustable chiral behavior under
modest magnetic fields. For such a prototypical example of non-hermitian and
dispersive photonic system, we demonstrate the capacity to magnetic field-tune
both the poles and zeros of cavity resonances, the two fundamental parameters
governing their resonance properties. Alongside the observed handedness-
dependent cavity frequencies, this highlights the remarkable ability to
engineer chiral and tunable radiative couplings for THz resonators and
metasurfaces. The extensive tunability offered by the surface plasmonic
approach paves the way for the development of agile and multifunctional THz
metasurfaces as well as the realization of ultrastrong chiral light-matter
interactions at low energy in matter with potential far-reaching applications
for the design of material properties.

## Introduction

Despite its long history, the role of chirality in light-matter interactions
remains intriguing, drawing intense focus both from a fundamental perspective
[1, 2] and for its technological applications in sensing and controlling
matter across diverse length scales, energy scales and scientific disciplines
[3, 4, 5, 6, 7, 8]. In this context, the design of chiral optical fields and
optically-active devices at THz frequencies holds particular significance for
sensing applications across a large class of systems. This stems from the wide
range of excitations exhibiting chirality within this frequency band, spanning
from chemically and biologically relevant excitations in molecules [9, 10] to
quantum collective modes such as magnons [4, 11, 12] and phonons [12, 13, 14]
in condensed matter systems. Alongside sensing capabilities, the ultrastrong
interaction between matter and chiral optical fields at low energies has
recently emerged as a new paradigm for manipulating material properties
through Floquet engineering [15, 16, 17] and has been proposed to realize
novel phases of matter through hybridization with the vacuum field of chiral
cavities [18, 19, 20, 21, 22]. Notably, the potential to induce magnetism and
non-trivial topologies in otherwise non-magnetic and topologically trivial
materials [18, 21, 23], or photon condensation in the ground state [24],
represent some of the most far-reaching prospects enabled by chiral vacuum
fields. For these reasons, the ability to engineer ultrastrong and chiral
light-matter interactions at THz frequencies, achievable through the design of
cavities with ultrasmall mode volumes, appears as a new frontier in the THz
range [25] that would open up many avenues for chiral sensing and matter
manipulation with applications in chemistry, biology and physics.
In the past years, large efforts have been devoted to the design of chiral and
optically-active devices for polarization control of THz light [26, 27, 28,
29, 30, 31, 32, 33, 34]. These have been obtained following mostly two
distinct routes. The first one is a synthetic or metamaterial-like approach,
in which the geometric chirality of metallic resonators is tailored and
transferred onto the resonating electromagnetic (EM) fields [26, 27, 28, 29,
30, 31]. The second approach relies on the natural magneto-optical activity of
the bulk of semiconductors [32, 33, 34] and is enabled by the propagation of
THz magneto-plasma waves [35, 36]. The latter originate from the microscopic
Lorentz force acting on carriers with low effective mass $m^{*}$, giving rise
to a cyclotron frequency $\nu_{c}=eB/2\pi m^{*}$ that falls within the THz
frequency range for moderate magnetic fields below 1T [37]. Yet, achieving
high EM-field confinement sufficient enough to reach the ultrastrong coupling
regime of light-matter interaction represents to date a challenge for both of
these approaches. On the one hand, shrinking the size of metallic resonators
while preserving operating frequencies in the THz range is a difficult task
that often requires the design of complex and technologically challenging
resonator architectures [38, 39, 40]. On the other hand, utilization of
magneto-plasma waves necessitates light propagation through the bulk of the
semiconductor and hence only offer sizeable chiral effects at the expense of
large interaction volumes [34, 36].
Here, we overcome these challenges by exploiting THz surface plasmons-based
cavities in semiconductors. As demonstrated recently [41], this approach
allows to downscale cavity dimensions at fixed operating frequencies and
provide an efficient way to confine THz EM-fields down to the ultimate limit
allowed by plasmonics. In this work, we show that surface plasmons inherit the
chiral activity of the bulk and use this as a means to create chiral THz
cavities with ultrasmall mode volumes. When assembled together, arrays of such
cavities lead to a novel type of chiral THz metasurfaces, fully tunable by
geometry, magnetic-field and even temperature. This work opens the door to
exciting perspectives for the exploration of ultrastrong chiral light-matter
interactions at low energies in matter.

## Working principle of the chiral THz surface plasmonic cavities

In the Faraday configuration, where a homogeneous magnetic field
$\boldsymbol{B}=B\boldsymbol{z}$ is perpendicular to the plasmonic medium’s
interface and a plane wave is incident normally, the magnetic field acts
exclusively onto the spin-sector of the EM-field (i.e. its polarization state)
and normal modes correspond to left- and right-circularly polarized (LCP/RCP)
photons
${\widetilde{\boldsymbol{E}}}_{{\sigma}_{\pm}}\boldsymbol{=}e^{ik_{z}z}e^{-i\omega
t}\left.\left|{\sigma}_{\pm}\right.\right\rangle$ [35]. Interaction between
light and the plasma involves two distinct magneto-plasmonic permittivities
$\epsilon_{{\sigma}\pm}(\nu,B)$ and results in reflection coefficients denoted
as $\tilde{r}_{{\sigma}_{\pm}}(\nu,B)$ (see supplementary information).
Normal incidence experimental power reflectivity spectra
$R_{{\sigma}_{\pm}}(\nu)=\arrowvert\tilde{r}_{{\sigma}_{\pm}}(\nu,B)\arrowvert^{2}$
as well as the phase
$\phi_{{\sigma}_{\pm}}(\nu,B)=Arg(\tilde{r}_{{\sigma}_{\pm}}(\nu,B))$
accumulated by the electric field of the light upon reflection onto the THz
plasma of bulk semiconductor InSb are shown in Fig. 1a for $\sigma_{\pm}$
polarizations. As the cyclotron frequency ${\nu}_{c}$ varies, the THz plasma
edge initially located at $\nu\approx{\nu}_{p}$ splits into two different
plasma edges. Their frequencies
${\nu}_{p,{\sigma}_{\pm}}\approx\frac{1}{2}\left(\sqrt{{4\nu}^{2}_{p}+{\nu}^{2}_{c}}\pm{\nu}_{c}\right)$
depend on the polarization handedness (see Fig. 1a-b) and highligths the bulk
chirality of the plasmonic medium [35, 36].
In contrast to plane waves, where orbital and spin components of the EM-field
are independent, the normal modes of the proposed ultrasmall chiral THz cavity
illustrated in Fig. 1c are sustained by surface plasmons and their orbital and
spin degrees of freedom are intertwined within the highly subwavelength volume
of the resonator, leading to a more complex chiral behavior. This can be seen
in Fig. 1d where we present the characteristics of the normal modes of
circular patch cavities, which, owing to their cylindrical symmetry, represent
a natural choice of plasmonic resonators in the present context. In the
absence of magnetic field, their normal modes have been studied previously
[42] and are characterized by an orbital angular moment $L\ (L=\
\dots,-2,-1,0,1,2,\dots)\ $ and a radial number $n\ (n=1,2,\ \dots)$. For
lossless cavities, modes of opposite orbital momentum are time-reversal
symmetric of one another
(${\widetilde{\boldsymbol{E}}}_{n,-L}\boldsymbol{=}{\widetilde{\boldsymbol{E}}}^{\boldsymbol{*}}_{n,L}$)
and hence have degenerate frequencies ${\nu}_{n,L}={\nu}_{n,-L}$. In Fig 1d,
we display the EM-fields sustained by the surface plasmon inside the resonator
for a set of representative $(n,L)$ modes: 1) the normal component of the
electric field ($Re(\tilde{E}_{z})$) and 2) the relative density
$\mathrm{\Delta}{\sigma}_{n,L}\left(r\right)={\sigma}_{+,n,L}\left(r\right)-{\sigma}_{-,n,L}\left(r\right)$
of the ${\sigma}_{\pm}$ polarized photons (or spin density of the mode), where
${\sigma}_{\pm,n,L}\left(r\right)\propto\mid\tilde{E}_{x}\pm
i\tilde{E}_{y}\mid^{2}$ depends only on the radial coordinate $r$.

Notably, normal modes are seen to consist in a spatially dependent
superposition of ${\sigma}_{\pm}$ polarized photons which entails the mixing
of orbital and spin degrees of freedom of the resonator. It is possible to
estimate the shifts $\delta{\nu}_{n,L}$ of the resonance frequencies
${\nu}_{n,L}$ by applying a modal perturbation theory, which, to first order
in ${\nu}_{c}$, gives (see supplementary material):

$\frac{\delta{\nu}_{n,L}}{{\nu}_{n,L}}\propto{\nu}_{c}\left(\int^{s/2}_{0}{r{\mathrm{\
}\mathrm{\Delta}\sigma}_{n,L}\left(r\right)dr}\right)$

indicating that the frequency shift is proportional to the spatially
integrated spin density ${\mathrm{\
}\mathrm{\Delta}\sigma}_{n,L}\left(r\right)$ of the ${\sigma}_{\pm}$ polarized
photons ($s$ being the diameter of the patch). As depicted in Fig. 1e, the
$L=0$ mode doesn’t experience any shift as it is linearly polarized, while
modes with finite orbital angular momentum $(L\neq 0)$ split symmetrically
according to $\delta{\nu}_{n,L}=-\delta{\nu}_{n,-L}$ as they have opposite
spin densities of equal magnitude. Because they carry the largest spin
polarization among all, the $(n=1,L=\pm 1)$ modes are expected to provide the
most pronounced mode splitting and hence the largest chiral effect.

## Sample architecture and plasmonic resonances without magnetic field

Rather than an isolated cavity, our samples consist in two-dimensional and
subwavelength periodic arrays (period $s+a$) of circular patch cavities of
diameter $s\ $(Fig. 2a). Over a bulk wafer of semiconducting InSb, we deposit
an insulating layer of highly subwavelength thickness $d$ (${Si}_{3}{N}_{4}$
with $d\ll\lambda_{0}$, $\lambda_{0}$ being the wavelength in vacuum) followed
by circular patches of metallic gold on top. In the absence of magnetic field
and due to the four-fold symmetry of the structure, reflectivity spectra are
independent of the incoming THz light polarization. In Fig. 2b, we present the
experimental and simulated power reflectivity spectra for three samples having
circular patches of varying diameters. The major absorption observed at low
frequency corresponds to the resonance frequency of the two degenerate
$(n,L)=(1,\pm 1)$ cavity modes, as confirmed by inspection of the mode profile
obtained from numerical simulations (Fig. 2c). A second pronounced absorption
at high frequency lies slightly below the THz plasma edge. As seen in Fig. 2c,
it corresponds to the resonance of a pair of degenerate plasmonic modes that
are fully expelled from the volume encompassed by the cavities and that are
located between them. These modes are the counterpart in these structures of
the inter-cavity mode observed in striped lattices of plasmonic cavities [41].
For completeness, we also report the presence of: 1) a small absorption for
the $s=43\mu m$ sample (triangle in Fig. 2b) which corresponds to the
excitation of the higher order $(n,L)=(3,\pm 1)$ cavity modes and whose mode
profile is shown in Fig.2c and 2) a resonance corresponding to the excitation
of a continuum of propagating bulk plasma modes within the transparency region
of the plasmonic medium arising the first-order diffraction of the array
($s=43\mu m$ sample, star in Fig. 2b). In the remaining of this paper, we will
not further discuss the latter resonances and will focus our analysis on the
main cavity resonance.

## Surface plasmonic chirality and its magnetic field dependence

Upon application of the magnetic field, we observe the lifting of the
degeneracy of the $\left(n,L\right)=(1,+1)$ and $(1,-1)$ cavity resonance.
This is featured in the polarization resolved THz spectra presented in Fig.
3a-b-c as two frequency-displaced resonant absorptions for opposite
${\sigma}_{\pm}$ polarizations, hence demonstrating the chiral effect that is
looked for.
Besides the polarization handedeness-dependent frequency-shifts, we notice
that the contrast of the cavity resonance varies strongly and evolves
oppositely for opposite ${\sigma}_{\pm}$ polarizations. The large variation of
its contrast can be inferred from the value of the reflectivity minimum at
resonance $R_{min}$ which changes as much as $\sim$50% within the range of
parameters investigated (see top and bottom panels in Fig. 3a-b-c).
This observation suggests that the magnetic field is acting to tune another
crucial cavity parameter: the overall coupling efficiency of external
radiation to the cavity. It can be understood qualitatively from temporal
coupled mode theory [43] which predicts that
$R_{min}={({\mathrm{\Gamma}}_{r}-{\mathrm{\Gamma}}_{\ohm})}^{2}/{({\mathrm{\Gamma}}_{r}+{\mathrm{\Gamma}}_{\ohm})}^{2}$
for single-port cavities [44, 45], where ${\mathrm{\Gamma}}_{r}$ and
${\mathrm{\Gamma}}_{\ohm}$ are the radiative and non-radiative decay rates of
the cavity describing respectively the coupling of the resonator to external
radiation and the decay of the resonance due to ohmic losses. In other words,
$R_{min}$ measures how close the cavity is from critical coupling at which
${\mathrm{\Gamma}}_{i}={\mathrm{\Gamma}}_{r}$ ($R_{min}=0$) and all incident
radiation is coupled and dissipated into the cavity. In order to further
assess quantitatively and reliably all cavity parameters for this highly non-
Hermitian and dispersive photonic system [46, 47, 48], we fitted the two
quadratures of the reflection coefficient $\tilde{r}_{{\sigma}_{\pm}}(\nu,B)$
based on an analytical continuation of the frequency in the complex plane [49,
50, 51] ($\nu\rightarrow\tilde{\nu}$ and
$\tilde{r}_{{\sigma}_{\pm}}(\nu,B)\rightarrow\tilde{r}_{{\sigma}_{\pm}}(\tilde{\nu},B)$,
see supplementary material). This leads to a pole-zero representation of the
cavity resonance shown in Fig. 4a. We define the complex zero
$\tilde{\nu}_{z0}$ by the condition
$\tilde{r}_{{\sigma}_{\pm}}(\tilde{\nu}_{z0})=0$. Following this zero in the
complex plane enables to characterize the undercoupled
($Im(\tilde{\nu}_{z0})<0)$), critically coupled ($Im(\tilde{\nu}_{z0})=0)$ and
overcoupled $(Im(\tilde{\nu}_{z0})>0)$ regimes. Note that the critically
coupled regime is closely linked to the coherent perfect absorption [52, 53,
54]. We observe that as the magnetic field is raised, the zero either moves
towards ($\sigma_{-}$) the real axis (critical coupling) or departs from it
($\sigma_{+}$), signaling not only the tunability of the coupling but also its
chiral behavior. The major results of the analysis of the cavity parameters
are further presented in Fig. 4b.
The relative frequency-shifts of the cavity resonance are defined as
$\frac{\delta{\nu}_{0}}{{\nu}_{0}}=\frac{{\nu}_{0,\sigma_{\pm}}(B)-{\nu}_{0}(B=0)}{{\nu}_{0}(B=0)}$
where ${\nu}_{0,\sigma_{\pm}}(B)=Re(\tilde{\nu}_{0,\sigma_{\pm}}(B))$
corresponds to the real part of the cavity pole
$\tilde{\nu}_{0,\sigma_{\pm}}(B)$, i.e. the cavity resonance frequency. As
anticipated, we observe that they demonstrate opposite variations for
$\sigma_{\pm}$ polarizations as a function of magnetic field $B$. While linear
in $B$ at small fields, they acquire a slight quadratic dependence at larger
fields which is reminiscent of the dependences of bulk magneto-plasma modes
[36]. At $B=500mT$, the relative shifts amount to 14$\%$, 20$\%$ and 24$\%$
for the $21,32$ and $43\mu m$ cavities respectively in between the two
$\sigma_{\pm}$ polarizations. The two decay rates, determined from the
knowledge of the pole and the zero as
$\Gamma_{\ohm}=-Im(\tilde{\nu}_{0}+\tilde{\nu}_{z0})/2$ and
$\Gamma_{r}=-Im(\tilde{\nu}_{0}-\tilde{\nu}_{z0})/2$ (see supplementary
material), also evolve as a function of magnetic-field and exhibit chiral
behaviour:
${\mathrm{\Gamma}}_{r/\ohm,\sigma_{+}}(B)\neq{\mathrm{\Gamma}}_{r/\ohm,\sigma_{-}}(B)$.
While this is somewhat expected for the ohmic decay rate as the cyclotron
resonance modulates material dissipation in the vicinity of $\nu=\pm\nu_{c}$,
the most striking feature is the chiral and large variations observed for the
radiative decay rate by as much as 300$\%$ in the range of parameters
investigated. Such phenomenon is largely unexpected as this quantity is
generally believed to depend solely on geometrical parameters of the arrays
[44, 45, 55]. We attribute it to the magnetic field-tuning of the phase
$\phi_{\sigma_{\pm}}(\nu,B)$ accumulated by the THz fields upon reflection
onto the underlying plasmonic medium. More specifically, denoting
$\tilde{r}_{\sigma_{\pm},0}(\nu,B)$ the background reflectivity of the
structure in the absence of metallic patches (i.e. involving only the
$vacuum/Si_{3}N_{4}/InSb$ planar multilayer), we show that the variation of
the radiative decay rate can be well reproduced and explained by a functionnal
dependence of the form
${\mathrm{\Gamma}}_{r,\sigma_{\pm}}(B)=\alpha\left|1+\tilde{r}_{\sigma_{\pm},0}(\nu,B)\right|^{2}$
where $\tilde{r}_{\sigma_{\pm},0}(\nu,B)$ is evaluated at the
$vacuum/Si_{3}N_{4}$ interface, at the cavity resonance frequency
$\nu=\nu_{0,\sigma_{\pm}}(B)$, and $\alpha$ is a geometry-dependent but
magnetic field-independent scaling factor (black solid lines in Fig. 4b). Such
a relation between background reflectivity and radiative decay rates has been
derived analytically for metal/insulator/metal cavities in a context where the
phase did not play any significant role [56, 57]. In the present case however,
the phase $\phi_{\sigma_{\pm}}(\nu,B)$ of the reflected field is the central
parameter controlling such behavior, as within a good approximation
$\left|\tilde{r}_{\sigma_{\pm},0}(\nu,B)\right|\approx 1$, so that
${\mathrm{\Gamma}}_{r,\sigma_{\pm}}(B)\propto\cos(\phi_{\sigma_{\pm}}(\nu,B)/2)^{2}$.
The basic idea behind the tunability and chirality of the radiative decay
rates that we observe is sketched in Fig. 4c. In the absence of metallic
patches, interference of the incident $E_{i,\sigma_{\pm}}$ and reflected
fields
$E_{r,\sigma_{\pm}}=\tilde{r}_{\sigma_{\pm},0}(\nu,B)E_{i,\sigma_{\pm}}$
result in a total background field of amplitude
$\left|E_{tot,\sigma_{\pm}}\right|=\left|1+\tilde{r}_{\sigma_{\pm},0}(\nu,B)\right|\left|E_{i,\sigma_{\pm}}\right|$
that takes the form of a standing wave whose locations of nodes and anti-nodes
are crucially depending the phase $\phi_{\sigma_{\pm}}(\nu,B)$ accumulated
upon reflection from the plasmonic medium. Immersed into this background
field, the metallic patches in- and out- couple it to the cavity with coupling
rates $\sqrt{2\mathrm{\Gamma}_{r,\sigma_{\pm}}(B)}$ [43] being determined by
the field amplitude at their location and hence scaling as
$\left|1+\tilde{r}_{\sigma_{\pm},0}(\nu,B)\right|$. This leads to an overall
radiation decay rate given by:
${\mathrm{\Gamma}}_{r,\sigma_{\pm}}(B)\propto\left|1+\tilde{r}_{\sigma_{\pm},0}(\nu,B)\right|^{2}$.
This phenomenon is analogous to a resonant nanoparticle located at distance
$L$ above a metallic mirror and whose coupling to external radiation
${\mathrm{\Gamma}}_{r}\propto\sin(2\pi L/\lambda_{0})^{2}$ is adjustable by
varying the location of the nanoparticle with respect to the standing wave
(see Fig. 4d). Here, despite the fixed geometry, the sizeable variations of
the phase provided by the underlying plasmonic medium close to the plasma
frequency as well as their tunability with magnetic field play the same role
as the height for the nanoparticle.
As such, for the proposed plasmonic cavities, we demonstrate that chirality
proceeds not only by frequency splitting of degenerate modes but also by a
chiral and tunable coupling of the resonators to external radiation. This
constitutes an important result of this study.

## Chiral THz surface plasmonic metasurfaces

We discuss here an application where the proposed structures can be exploited
as resonant and tunable reflective THz chiral metasurfaces for polarization
control of the THz light. In order to demonstrate the advantages of
structuring a metasurface over that of using a bulk semiconductor for this
purpose, we compare in Fig. 5 the optical activity of the unstructured bulk
material and that of a structured metasurface. Optical activity is commonly
characterized in terms of: 1) circular dichroism (CD) via the ellipticity
$\eta=\frac{|r_{\sigma_{+}}(\nu)|-|r_{\sigma_{-}}(\nu)|}{|r_{\sigma_{+}}(\nu)|+|r_{\sigma_{-}}(\nu)|}$
and 2) optical rotary dispersion (OR) via the rotation angle
$\theta=\frac{\phi_{\sigma_{+}}(\nu)-\phi_{\sigma_{-}}(\nu)}{2}$[58]. For
linearly polarized incident light, CD and OR result in the present context in
reflected light that is elliptically polarized (CD phenomenon) and rotated by
an angle $\theta$ compared to the incident polarization (OR phenomenon), see
Fig. 5a. In Fig. 5b, we compare the ellipticity and optical rotation angle of
the bulk and structured metasurface at $B=500mT$.
We notice profound differences between the two responses. Firstly, we observe
that CD and OR of the bulk material display rather complex shapes in an
otherwise large frequency window centered around the plasma edge. This results
in overall poor spectral control of the optical activity of the bulk material.
In contrast, the ellipticity of the metasurface presents two sizeable peaks of
opposite polarity at around 1THz and 2THz, which correspond to the frequencies
of the cavity and intercavity resonances respectively for this $s=32\mu m$
cavity metasurface. Ellipticity around the cavity resonance demonstrates in
particular a well-defined unipolar behaviour.
Strikingly, at the cavity resonance, the ellipticity reaches a sizeable value
$\eta\approx 0.5$ that is at a frequency where the bulk ellipticity is almost
null within noise level: this corresponds to a giant enhancement of the
ellipticity provided by the metasurface. The ability to induce such a large
dichroism in a frequency band well below the plasma frequency where the bulk
material only provides a minute effect demonstrates the superiority of
metasurface over the bulk material for the chiral engineering of THz light
fields. The origin of such a metasurface-enabled ellipticity can be traced
back to the difference in the contrast between $\sigma_{+}$ and $\sigma_{-}$
resonances. In this respect, chirality of the radiative decay rates of the
cavities discussed in the previous section is crucial in providing this
effect.
In the supplementary material, we further show the wide tunability of the
optical activity of the metasurfaces both spectrally via the geometrical
tuning of the cavity dimensions and in its magnitude via the tuning of the
magnetic field strenght. Such metasurfaces achieve ellipticities as high as
$\eta\approx 0.7$ at a moderate field $B=500mT$ relevant for applications, and
we expect that full ellipticity ($\eta=1$) is reachable around $600-700mT$.

## Discussion

Our approach offers a major and unique asset with respect to other
alternatives for the chiral engineering of THz fields based on either
synthetic [26, 28, 29, 31] or natural [36, 32, 34] optical activity. In
comparison, the large magnitude of optical activity obtained here is achieved
over much smaller, i.e. deep-subwavelength, scales thanks to the remarkable
confinement properties of surface plasmons. All the reported chiral effects
originate from beneath the surface of the plasmonic medium within the
penetration depth of the surface plasmons that is $\sim 0.01\lambda_{0}$
(where $\lambda_{0}$ is the free space wavelength) and within cavities with
mode volumes of the order of $10^{-6}\lambda_{0}^{3}$. This represents at
least two orders of magnitude enhancement over previous approaches for typical
figure of merits that compare the magnitude of CD/OR and the lenght or volume
over which it is realized. Straightforward improvements can be readily
achieved by reducing cavity mode volumes to $10^{-8}\lambda_{0}^{3}$ [41] or
by increasing the ratio $\nu_{c}/\nu_{p}$ through the magnetic field or the
temperature [41], as this ratio governs the overall magnitude of the chiral
effects. In summary, along with the numerous opportunities it allows for
harnessing ultrastrong chiral light-matter interactions for the sensing and
control of matter at THz frequencies, this system stands out as a singular
platform within the electromagnetic spectrum for unravelling novel phenomena
in chiral plasmonics thanks to the extended parameter space in
$\nu_{c}/\nu_{p}$ which remains beyond reach at mid-infrared and optical
frequencies.

## Acknowledgments

Y.L. acknowledges support from Agence Nationale de la Recherche (grant n∘ ANR-
17-CE30-0011-11) and Labex PALM (grant n∘ ANR-10-LABX-0039-PALM).

## Methods

### Sample manufacturing

Samples were manufactured starting from a $500\ \mu m$ thick, $<100>$
oriented, bulk wafer of InSb commercially available (MTI Corporation). The
InSb wafer was nominally N-type undoped. Each sample consisted of a $\sim 5\
mm*5\ mm*0.5\ mm$ InSb substrate dice-cut from the wafer. An insulating layer
of Si3N4 was deposited via plasma enhanced chemical vapor deposition and the
top metallic patches were realized via photolithography followed by metal
deposition of Ti(15 nm)/Au(200 nm). The geometric and optical parameters of
the samples are provided in the supplementary information.

### THz spectroscopy

Spectroscopy of the samples was performed with a time-domain THz spectrometer
driven by an ultrafast Ti:Sapphire oscillator. THz generation and detection
were achieved via a photoconductive emitter (Tera-SED) and a 1 mm thick
$<110>$ ZnTe crystal, respectively, allowing for spectroscopic coverage from
0.2 THz to 2.5 THz. Delayed THz pulses originating from the optical components
of the setup were eliminated by windowing the main THz pulse reflected off the
samples. The time domain signals were Fourier transformed and analyzed in both
their amplitude and phase. All reflectivity measurements were conducted at
normal incidence and room temperature. For absolute measurements of the
reflectivity, we used the THz pulse reflected off a gold-coated sample as a
reference. In all measurements, the THz spot size ($\sim 2\ mm$ in diameter)
is well below the size of the manufactured samples, so that the samples can be
considered infinite in the plane containing the interface ($x$-$y$ plane in
Fig. 2.a). Polarization-resolved THz measurements were obtained by initially
decomposing the linearly-polarized incident THz pulse into the $\sigma_{\pm}$
circular basis. The reflection coefficient of each $\sigma_{\pm}$ component
was obtained thanks to a polarization-resolved detection. A description of the
experimental setup and of the polarization-resolved THz measurements is
reported in the supplementary information. The magnetic field was applied by a
large static magnet located behind the sample and was calibrated with a Hall
probe. Homogeneity of the magnetic-field and the uncertainty on its absolute
magnitude were better than $3\%$ across the THz spot size.

### Numerical simulations

Numerical simulations of the cavities presented in Fig. 2b were performed with
the Rigorous Coupled Wave Analysis method [59, 60]. The optical and
geometrical parameters used in the simulations are determined from
experimental measurements performed on the samples (see supplementary
information). The permittivity of Gold at THz frequencies was taken from [61]
and the permittivity of $Si_{3}N_{4}$ from [41].

Figure 1: a) ${\sigma}_{\pm}\ $polarization-resolved power spectra (top panel)
and phase (bottom panel) of the reflected THz field from bulk InSb at various
magnetic fields between $0mT$ and $500mT$. Measurements are performed at
normal incidence and the magnetic field is applied perpendicular to the
surface of the plasmonic medium. The dashed line indicates the location of
plasma frequency $\nu_{p}$ at $B=0mT$. b) Cyclotron frequency (black dots) as
a function of magnetic field, determined from fitting the bulk reflectivity
spectra with the magneto-plasmonic model of permittivities (see also
supplementary material). Red and blue dots indicate handedness-dependent
plasma frequencies ${\nu}_{p,{\sigma}_{\pm}}$ associated with the cyclotron
resonance active and inactive bulk magneto-plasma waves ([35]) c) Sketch of a
single subwavelength chiral THz plasmonic cavity with overall dimensions
$s,d\ll\lambda_{0}$. d) Characteristics of the plasmonic modes for some
selected cavity resonances of the circular patch cavity without a magnetic
field: $(n,L)=\\{\left(1,0\right),(1,+1),\ (1,-1),\ (3,+1),\ (3,-1)\\}$. Shown
are color plots of: 1) the normal component of the EM-field
($Re(\tilde{E}_{z})$) (top) and 2) relative density of the ${\sigma}_{\pm}$
polarized photons (or ’spin’ density, bottom) of the corresponding surface
plasmon e) Resonance splitting expected for the different orbital modes of the
subwavelength cavity as a function of the cyclotron frequency.

Figure 2: Figure 2: a) Sample architecture and geometrical parameters of the
two-dimensional array of circular patch cavities (right). Top-view optical
pictures of the actual samples (left). b) Experimental (left) and simulated
(right) normal incidence THz power reflectivity spectra of the 3 samples
investigated at zero magnetic field ($B=0mT$): $s=43\mu m$ (orange solid
curve), $s=32\mu m$ (purple solid curve), $s=21\mu m$ (green solid curve). The
normal incidence THz power reflectivity spectrum of bulk InSb is reported for
comparison (black solid curve). The different plasmonic resonances are
reported: main cavity resonance (dots), second cavity resonance (triangle),
intercavity resonance (squares), first-order diffraction of the array (star)
c) Simulated EM mode profiles of the $s=43\mu m$ sample performed at
frequencies corresponding to the resonances indicated in panel b. All mode
profiles are shown for a $z-$cut taken in the middle of the dielectric
insulator at $z=d/2$. Top: Mode profile of the main cavity resonance (plot of
$Re(\tilde{E}_{z})$). Center: Mode profile of the second cavity resonance
(plot of $Re(\tilde{E}_{z})$). Bottom: Mode profile of the inter-cavity
resonance (plots of $|\tilde{E}_{x}|$ and $|\tilde{E}_{y}|$)

Figure 3: Experimental measurements of the ${\sigma}_{\pm}$ polarization-
resolved THz power reflectivity spectra as a function of magnetic field
between $0mT$ and $500mT$ for the three samples investigated: a) $s=21\mu m$
b) $s=32\mu m$ and c) $s=43\mu m$. The central panels display the THz
reflectivity spectra in color scale. For each magnetic field, the
corresponding cyclotron frequency is reported on the right axis. The top and
bottom panels show the THz spectra at $B=500mT$ obtained for the
${\sigma}_{+}$ (bottom panel, red solid curves) and ${\sigma}_{-}$ (top panel,
blue solid curves) polarizations, together with the THz spectrum obtained at
$B=0mT$ for comparison (grey solid curves).

Figure 4: a) Location of the poles ($\tilde{\nu}_{0}$) and zeros
($\tilde{\nu}_{z0}$) of the cavity resonance in the complex plane and their
trajectory under application of the magnetic field. b) Chirality of the cavity
parameters as a function of magnetic field for $\sigma_{+}$ (red dots) and
$\sigma_{-}$ (blue dots) polarizations. Top panels: relative frequency-shifts
$\delta{\nu}_{0,\sigma_{\pm}}/{\nu}_{0}$. Bottom panels: radiative decay rates
${\mathrm{\Gamma}}_{r,\sigma_{\pm}}$. The experimental data (dots) is compared
to the model describing the variation of the radiative decay rates in a
magnetic field (black solid lines):
${\mathrm{\Gamma}}_{r,\sigma_{\pm}}(B)=\alpha\left|1+\tilde{r}_{\sigma_{\pm},0}(\nu,B)\right|^{2}$(see
main text). From left to right: $s=43,32$ and $21\mu m$ cavity samples. c
Sketch of the origin of the chirality and of the tuning of the radiative decay
rates with the magnetic field. In the absence of the gold patch, incoming and
reflected light from the plasmonic medium interfere to form standing waves
with total amplitude $\left|E_{tot,\sigma_{\pm}}\right|$ (the background
field) whose location of nodes and anti-nodes depend on the phase
$\phi_{\sigma_{\pm}}(\nu,B)$ imposed by the medium upon reflection, as shown
by the dashed lines which indicate the location of the anti-nodes of the
background field. When the gold patch is then introduced, it acts both as a
support to the cavity resonance and as a coupler to it. Its location with
respect to the background field determines how much light is coupled in and
out of the cavity. The chirality and tunability of the phase
$\phi_{\sigma_{\pm}}(\nu,B)$ of the bulk plasmonic medium, as shown in Fig.
1a, is hence rooting that of the radiative decay rates. d Analogy of the
tuning of the radiative decay rates with a resonant nanoparticle located above
a metallic mirror. In the absence of the nanoparticle, incoming and reflected
light from the metallic mirror form standing waves. When the nanoparticle is
then immersed into this background field, its location determines the amount
of light that is scattered and used to drive its resonance, i.e. the radiative
decay rates. If it is located at the anti-node of the background field, the
absence of light inhibits the resonance, while if it is located at a node, the
resonance is maximally driven. The height of the nanoparticle controls in this
case the radiative decay rates and tunability is achieved geometrically.

Figure 5: Comparison between the optical activity arising from the bulk
semiconductor and a THz metasurface. a Sketch of the change in the
polarization of the light as due to optical activity from either the bulk
semiconductor or a structured metasurface (here shown is the metasurface).
Optical activity manifests in circular dichroism (CD) and optical rotary
dispersion (OD). For linearly polarized incident light, CD results in
elliptical polarization while OR results in optical rotation of the
polarization by an angle $\theta$ of the reflected light. b Comparision of CD
and OR between the bulk semiconductor (top panel) and $s=32\mu m$ metasurface
(bottom panel) at $B=500mT$.

## References

  * [1] Yiqiao Tang and Adam E. Cohen.  Optical chirality and its interaction with matter.  Physical Review Letters, 104(16):163901, April 2010.
  * [2] Jungho Mun, Minkyung Kim, Younghwan Yang, Trevon Badloe, Jincheng Ni, Yang Chen, Cheng-Wei Qiu, and Junsuk Rho.  Electromagnetic chirality: from fundamentals to nontraditional chiroptical phenomena.  Light: Science and Applications, 9(1), September 2020.
  * [3] Yiqiao Tang and Adam E. Cohen.  Enhanced enantioselectivity in excitation of chiral molecules by superchiral light.  Science, 332(6027):333–336, April 2011.
  * [4] S. Bordács, I. Kézsmárki, D. Szaller, L. Demkó, N. Kida, H. Murakawa, Y. Onose, R. Shimano, T. Rõõm, U. Nagel, S. Miyahara, N. Furukawa, and Y. Tokura.  Chirality of matter shows up via spin excitations.  Nature Physics, 8(10):734–738, August 2012.
  * [5] Peter Lodahl, Sahand Mahmoodian, Søren Stobbe, Arno Rauschenbeutel, Philipp Schneeweiss, Jürgen Volz, Hannes Pichler, and Peter Zoller.  Chiral quantum optics.  Nature, 541(7638):473–480, jan 2017.
  * [6] Jochen R. Brandt, Francesco Salerno, and Matthew J. Fuchter.  The added value of small-molecule chirality in technological applications.  Nature Reviews Chemistry, 1(6), June 2017.
  * [7] Andrew Lininger, Giovanna Palermo, Alexa Guglielmelli, Giuseppe Nicoletta, Madhav Goel, Michael Hinczewski, and Giuseppe Strangi.  Chirality in light–matter interaction.  Advanced Materials, 35(34), May 2022.
  * [8] Cyriaque Genet.  Chiral light–chiral matter interactions: an optical force perspective.  ACS Photonics, 9(2):319–332, January 2022.
  * [9] Won Jin Choi, Keiichi Yano, Minjeong Cha, Felippe M. Colombari, Ji-Young Kim, Yichun Wang, Sang Hyun Lee, Kai Sun, John M. Kruger, André F. de Moura, and Nicholas A. Kotov.  Chiral phonons in microcrystals and nanofibrils of biomolecules.  Nature Photonics, 16(5):366–373, mar 2022.
  * [10] Won Jin Choi, Sang Hyun Lee, Bum Chul Park, and Nicholas A. Kotov.  Terahertz circular dichroism spectroscopy of molecular assemblies and nanostructures.  Journal of the American Chemical Society, 144(50):22789–22804, dec 2022.
  * [11] Tobias Kampfrath, Alexander Sell, Gregor Klatt, Alexej Pashkin, Sebastian Mährlein, Thomas Dekorsy, Martin Wolf, Manfred Fiebig, Alfred Leitenstorfer, and Rupert Huber.  Coherent terahertz control of antiferromagnetic spin waves.  Nature Photonics, 5(1):31–34, nov 2010.
  * [12] Jun Cui, Emil Viñas Boström, Mykhaylo Ozerov, Fangliang Wu, Qianni Jiang, Jiun-Haw Chu, Changcun Li, Fucai Liu, Xiaodong Xu, Angel Rubio, and Qi Zhang.  Chirality selective magnon-phonon hybridization and magnon-induced chiral phonons in a layered zigzag antiferromagnet.  Nature Communications, 14(1), jun 2023.
  * [13] Kyosuke Ishito, Huiling Mao, Yusuke Kousaka, Yoshihiko Togawa, Satoshi Iwasaki, Tiantian Zhang, Shuichi Murakami, Jun ichiro Kishine, and Takuya Satoh.  Truly chiral phonons in $\alpha$-HgS.  Nature Physics, 19(1):35–39, oct 2022.
  * [14] Andrey Baydin, Felix G. G. Hernandez, Martin Rodriguez-Vega, Anderson K. Okazaki, Fuyang Tay, G. Timothy Noe, Ikufumi Katayama, Jun Takeda, Hiroyuki Nojiri, Paulo H. O. Rappl, Eduardo Abramof, Gregory A. Fiete, and Junichiro Kono.  Magnetic control of soft chiral phonons in PbTe.  Physical Review Letters, 128(7):075901, feb 2022.
  * [15] Y. H. Wang, H. Steinberg, P. Jarillo-Herrero, and N. Gedik.  Observation of floquet-bloch states on the surface of a topological insulator.  Science, 342(6157):453–457, October 2013.
  * [16] J. W. McIver, B. Schulte, F.-U. Stein, T. Matsuyama, G. Jotzu, G. Meier, and A. Cavalleri.  Light-induced anomalous hall effect in graphene.  Nature Physics, 16(1):38–41, November 2019.
  * [17] Gabriel E. Topp, Gregor Jotzu, James W. McIver, Lede Xian, Angel Rubio, and Michael A. Sentef.  Topological floquet engineering of twisted bilayer graphene.  Physical Review Research, 1(2):023031, September 2019.
  * [18] Xiao Wang, Enrico Ronca, and Michael A. Sentef.  Cavity quantum electrodynamical chern insulator: Towards light-induced quantized anomalous hall effect in graphene.  Physical Review B, 99(23):235156, June 2019.
  * [19] Hannes Hübener, Umberto De Giovannini, Christian Schäfer, Johan Andberger, Michael Ruggenthaler, Jerome Faist, and Angel Rubio.  Engineering quantum materials with chiral optical cavities.  Nature Materials, 20(4):438–442, nov 2020.
  * [20] I. V. Tokatly, D. R. Gulevich, and I. Iorsh.  Vacuum anomalous hall effect in gyrotropic cavity.  Physical Review B, 104(8):l081408, August 2021.
  * [21] D. D. Sedov, V. Shirobokov, I. V. Iorsh, and I. V. Tokatly.  Cavity-induced chiral edge currents and spontaneous magnetization in two-dimensional electron systems.  Physical Review B, 106(20):205114, nov 2022.
  * [22] John Clai Owens, Margaret G. Panetta, Brendan Saxberg, Gabrielle Roberts, Srivatsan Chakram, Ruichao Ma, Andrei Vrajitoarea, Jonathan Simon, and David I. Schuster.  Chiral cavity quantum electrodynamics.  Nature Physics, 18(9):1048–1052, jul 2022.
  * [23] Kanta Masuki and Yuto Ashida.  Berry phase and topology in ultrastrongly coupled quantum light-matter systems.  Physical Review B, 107(19):195104, May 2023.
  * [24] Alberto Mercurio, Gian Marcello Andolina, Francesco M. D. Pellegrino, Omar Di Stefano, Pablo Jarillo-Herrero, Claudia Felser, Frank H. L. Koppens, Salvatore Savasta, and Marco Polini.  Photon condensation, van vleck paramagnetism, and chiral cavities.  Physical Review Research, 6(1):013303, March 2024.
  * [25] Johan Andberger, Lorenzo Graziotto, Luca Sacchi, Mattias Beck, Giacomo Scalari, and Jerome Faist.  Terahertz chiral metamaterial cavities breaking time-reversal symmetry.  arXiv.2308.03195, 2023.
  * [26] Shuang Zhang, Jiangfeng Zhou, Yong-Shik Park, Junsuk Rho, Ranjan Singh, Sunghyun Nam, Abul K. Azad, Hou-Tong Chen, Xiaobo Yin, Antoinette J. Taylor, and Xiang Zhang.  Photoinduced handedness switching in terahertz chiral metamolecules.  Nature Communications, 3(1), jul 2012.
  * [27] Jianfeng Wu, Binghao Ng, Haidong Liang, Mark B.H. Breese, Minghui Hong, Stefan A. Maier, Herbert O. Moser, and Ortwin Hess.  Chiral metafoils for terahertz broadband high-contrast flexible circular polarizers.  Physical Review Applied, 2(1):014005, July 2014.
  * [28] Tetsuo Kan, Akihiro Isozaki, Natsuki Kanda, Natsuki Nemoto, Kuniaki Konishi, Hidetoshi Takahashi, Makoto Kuwata-Gonokami, Kiyoshi Matsumoto, and Isao Shimoyama.  Enantiomeric switching of chiral metamaterial for terahertz polarization modulation employing vertically deformable MEMS spirals.  Nature Communications, 6(1), oct 2015.
  * [29] Teun-Teun Kim, Sang Soon Oh, Hyeon-Don Kim, Hyun Sung Park, Ortwin Hess, Bumki Min, and Shuang Zhang.  Electrical access to critical coupling of circularly polarized waves in graphene chiral metamaterials.  Science Advances, 3(9), September 2017.
  * [30] Won Jin Choi, Gong Cheng, Zhengyu Huang, Shuai Zhang, Theodore B. Norris, and Nicholas A. Kotov.  Terahertz circular dichroism spectroscopy of biomaterials enabled by kirigami polarization modulators.  Nature Materials, 18(8):820–826, jul 2019.
  * [31] Longqing Cong, Prakash Pitchappa, Nan Wang, and Ranjan Singh.  Electrically programmable terahertz diatomic metamolecules for chiral optical control.  Research, 2019:1–11, February 2019.
  * [32] Qianyi Mu, Fei Fan, Sai Chen, Shitong Xu, Chuanzhong Xiong, Xin Zhang, Xianghui Wang, and Shengjiang Chang.  Tunable magneto-optical polarization device for terahertz waves based on insb and its plasmonic structure.  Photonics Research, 7(3):325, February 2019.
  * [33] Zhiyu Tan, Fei Fan, Dan Zhao, Shanshan Li, Xianghui Wang, and Shengjiang Chang.  Linear-polarized terahertz isolator by breaking the gyro-mirror symmetry in cascaded magneto-optical metagrating.  Nanophotonics, 10(16):4141–4148, October 2021.
  * [34] Xuewei Ju, Zhiqiang Hu, Guofeng Zhu, Feng Huang, Yanqing Chen, Cuixia Guo, Alexey Belyanin, Junichiro Kono, and Xiangfeng Wang.  Creating a near-perfect circularly polarized terahertz beam through the nonreciprocity of a magnetoplasma.  Optics Express, 31(23):38540, October 2023.
  * [35] E D Palik and J K Furdyna.  Infrared and microwave magnetoplasma effects in semiconductors.  Reports on Progress in Physics, 33(3):1193–1322, September 1970\.
  * [36] X. Wang, A. A. Belyanin, S. A. Crooker, D. M. Mittleman, and J. Kono.  Interference-induced terahertz transparency in a semiconductor magneto-plasma.  Nature Physics, 6(2):126–130, December 2009.
  * [37] R. Shimano, Y. Ino, Yu. P. Svirko, and M. Kuwata-Gonokami.  Terahertz frequency hall measurement by magneto-optical kerr spectroscopy in inas.  Applied Physics Letters, 81(2):199–201, July 2002.
  * [38] B. Paulillo, J. M. Manceau, A. Degiron, N. Zerounian, G. Beaudoin, I. Sagnes, and R. Colombelli.  Circuit-tunable sub-wavelength thz resonators: hybridizing optical cavities and loop antennas.  Optics Express, 22(18):21302, August 2014.
  * [39] Janine Keller, Giacomo Scalari, Sara Cibella, Curdin Maissen, Felice Appugliese, Ennio Giovine, Roberto Leoni, Mattias Beck, and Jérôme Faist.  Few-electron ultrastrong light-matter coupling at 300 ghz with nanogap hybrid lc microcavities.  Nano Letters, 17(12):7410–7415, November 2017.
  * [40] Alireza Mottaghizadeh, Yanko Todorov, Mathis Cameau, Djamal Gacemi, Angela Vasanelli, and Carlo Sirtori.  Nanoscale electromagnetic confinement in thz circuit resonators.  Optics Express, 25(23):28718, November 2017.
  * [41] Ian Aupiais, Romain Grasset, Tingwen Guo, Dmitri Daineka, Javier Briatico, Sarah Houver, Luca Perfetti, Jean-Paul Hugonin, Jean-Jacques Greffet, and Yannis Laplace.  Ultrasmall and tunable terahertz surface plasmon cavities at the ultimate plasmonic limit.  Nature Communications, 14(1), November 2023.
  * [42] Fred Minkowski, Feng Wang, Ayan Chakrabarty, and Qi-Huo Wei.  Resonant cavity modes of circular plasmonic patch nanoantennas.  Applied Physics Letters, 104(2):021111, jan 2014.
  * [43] Shanhui Fan, Wonjoo Suh, and J. D. Joannopoulos.  Temporal coupled-mode theory for the fano resonance in optical resonators.  Journal of the Optical Society of America A, 20(3):569, March 2003\.
  * [44] J.-M. Manceau, S. Zanotto, I. Sagnes, G. Beaudoin, and R. Colombelli.  Optical critical coupling into highly confining metal-insulator-metal resonators.  Applied Physics Letters, 103(9), August 2013.
  * [45] Che Qu, Shaojie Ma, Jiaming Hao, Meng Qiu, Xin Li, Shiyi Xiao, Ziqi Miao, Ning Dai, Qiong He, Shulin Sun, and Lei Zhou.  Tailor the functionalities of metasurfaces based on a complete phase diagram.  Physical Review Letters, 115(23):235503, December 2015.
  * [46] Filippo Alpeggiani, Nikhil Parappurath, Ewold Verhagen, and L. Kuipers.  Quasinormal-mode expansion of the scattering matrix.  Physical Review X, 7(2):021035, June 2017.
  * [47] Hanwen Zhang and Owen D. Miller.  Quasinormal coupled mode theory.  arXiv:2010.08650, 2020.
  * [48] Mohammed Benzaouia, John D. Joannopoulos, Steven G. Johnson, and Aristeidis Karalis.  Quasi-normal mode theory of the scattering matrix, enforcing fundamental constraints for truncated expansions.  Physical Review Research, 3(3):033228, September 2021.
  * [49] E. Popov, L. Mashev, and D. Maystre.  Theoretical study of the anomalies of coated dielectric gratings.  Optica Acta: International Journal of Optics, 33(5):607–619, May 1986.
  * [50] M. Nevière, R. Reinisch, and E. Popov.  Electromagnetic resonances in linear and nonlinear optics: phenomenological study of grating behavior through the poles and zeros of the scattering operator.  Journal of the Optical Society of America A, 12(3):513, March 1995\.
  * [51] Felix Binkowski, Fridtjof Betz, Rémi Colom, Patrice Genevet, and Sven Burger.  Poles and zeros in non-hermitian systems: Application to photonics.  Physical Review B, 109(4):045414, January 2024.
  * [52] Y. D. Chong, Li Ge, Hui Cao, and A. D. Stone.  Coherent perfect absorbers: Time-reversed lasers.  Physical Review Letters, 105(5):053901, July 2010.
  * [53] Wenjie Wan, Yidong Chong, Li Ge, Heeso Noh, A. Douglas Stone, and Hui Cao.  Time-reversed lasing and interferometric control of absorption.  Science, 331(6019):889–892, February 2011.
  * [54] William R. Sweeney, Chia Wei Hsu, and A. Douglas Stone.  Theory of reflectionless scattering modes.  Physical Review A, 102(6):063511, December 2020.
  * [55] Cheryl Feuillet-Palma, Yanko Todorov, Angela Vasanelli, and Carlo Sirtori.  Strong near field enhancement in thz nano-antenna arrays.  Scientific Reports, 3(1), March 2013.
  * [56] Patrick T. Bowen, Alexandre Baron, and David R. Smith.  Theory of patch-antenna metamaterial perfect absorbers.  Physical Review A, 93(6):063849, June 2016.
  * [57] P. T. Bowen, A. Baron, and D. R. Smith.  Effective-medium description of a metasurface composed of a periodic array of nanoantennas coupled to a metallic film.  Physical Review A, 95(3):033822, March 2017.
  * [58] Jiangfeng Zhou, Dibakar Roy Chowdhury, Rongkuo Zhao, Abul K. Azad, Hou-Tong Chen, Costas M. Soukoulis, Antoinette J. Taylor, and John F. O’Hara.  Terahertz chiral metamaterials with giant and dynamically tunable optical activity.  Physical Review B, 86(3):035448, jul 2012.
  * [59] M. G. Moharam, T. K. Gaylord, Eric B. Grann, and Drew A. Pommet.  Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings.  Journal of the Optical Society of America A, 12(5):1068, May 1995\.
  * [60] Jean Paul Hugonin and Philippe Lalanne.  Reticolo software for grating analysis.  arXiv:2101.00901, 2021.
  * [61] Y. Todorov, L. Tosetto, J. Teissier, A. M. Andrews, P. Klang, R. Colombelli, I. Sagnes, G. Strasser, and C. Sirtori.  Optical properties of metal-dielectric-metal microcavities in the thz frequency range.  Optics Express, 18(13):13886, June 2010.


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

# Limited Resources Induce Bistability in Microtubule Length Regulation

Matthias Rank  Arnold-Sommerfeld-Center for Theoretical Physics and Center for
NanoScience, Ludwig-Maximilians-Universität München, Theresienstraße 37, 80333
München, Germany  Aniruddha Mitra  B CUBE – Center for Molecular
Bioengineering and Center for Advancing Electronics Dresden (cfaed),
Technische Universität Dresden, Arnoldstraße 18, 01307 Dresden, Germany  Max
Planck Institute of Molecular Cell Biology and Genetics, Pfotenhauerstraße
108, 01307 Dresden, Germany  Louis Reese  Department of Bionanoscience, Kavli
Institute of Nanoscience, Faculty of Applied Sciences, Delft University of
Technology, Van der Maasweg 9, 2629 HZ Delft, The Netherlands  Stefan Diez
<EMAIL_ADDRESS>B CUBE – Center for Molecular Bioengineering and
Center for Advancing Electronics Dresden (cfaed), Technische Universität
Dresden, Arnoldstraße 18, 01307 Dresden, Germany  Max Planck Institute of
Molecular Cell Biology and Genetics, Pfotenhauerstraße 108, 01307 Dresden,
Germany  Erwin Frey<EMAIL_ADDRESS>Arnold-Sommerfeld-Center for Theoretical
Physics and Center for NanoScience, Ludwig-Maximilians-Universität München,
Theresienstraße 37, 80333 München, Germany

###### Abstract

The availability of protein is an important factor for the determination of
the size of the mitotic spindle. Involved in spindle-size regulation is
kinesin-8, a molecular motor and microtubule (MT) depolymerase, which is known
to tightly control MT length. Here, we propose and analyze a theoretical model
in which kinesin-induced MT depolymerization competes with spontaneous
polymerization while supplies of both tubulin and kinesin are limited. In
contrast to previous studies where resources were unconstrained, we find that,
for a wide range of concentrations, MT length regulation is bistable. We test
our predictions by conducting _in vitro_ experiments, and find that the
bistable behavior manifests in a bimodal MT length distribution.

The absolute and relative abundance of particular sets of proteins is
important for a wide range of processes in cells. For example, during _Xenopus
laevis_ embryogenesis, importin $\alpha$ becomes progressively localized to
the cell membrane Wilbur and Heald (2013). As a consequence of importin’s
depletion from the cytoplasm, the protein kif2a escapes inactivation, and
decreases the size of the mitotic spindle. Similarly, formation of the mitotic
spindle reduces the concentration of free tubulin dimers, the building blocks
of microtubules (MTs). Thus, up to $60\%$ of all tubulin heterodimers Borisy
and Taylor (1967); Borisy (1967) may be incorporated into the spindle Good
_et al._ (2013). In addition, it has been shown _in vivo_ and _in vitro_ that
both spindle size Good  _et al._ (2013); Hazel  _et al._ (2013) and the
lengths of its constituent MTs Winey  _et al._ (1995) scale with cytoplasmic
volume.

Assembly and disassembly of MTs are regulated by a set of proteins that
interact with the plus ends of protofilaments Howard and Hyman (2003, 2007).
One of these factors, the molecular motor kinesin-8, acts as a depolymerase
Varga  _et al._ (2006); Howard and Hyman (2007). As a consequence, spindle
size increases in its absence Goshima  _et al._ (2005), and decreases upon
overexpression of the protein Stumpff  _et al._ (2008). Moreover, the
kinesin-8 homolog Kip3 from _Saccharomyces cerevisiae_ has been shown to
depolymerize MTs in a length-dependent fashion Varga  _et al._ (2006, 2009).
This is facilitated by a density gradient on the MT, caused by the interplay
between the processive motion of Kip3 along the MT and its depolymerase
activity at the plus end, which effectively enables the MT to “sense” its own
length Varga  _et al._ (2009); Reese  _et al._ (2011). In combination with
spontaneous MT polymerization, the Kip3 gradient leads to a length regulation
mechanism Melbinger  _et al._ (2012); Reese  _et al._ (2014).

Here we explore the combined effect of limited resources and Kip3-induced
depolymerization on the length regulation of MTs. As seen in theoretical
studies on the collective motion of molecular motors, resource limitation
affects the density profile on the MT: Regions of low and high motor density
separate, as a localized domain wall emerges on the MT Lipowsky  _et al._
(2001); Cook and Zia (2009); Greulich  _et al._ (2012); Ciandrini  _et al._
(2014). This is a direct result of resource limitation, and does not rely on
the existence of a motor density gradient, as necessary for domain wall
localization in the presence of unlimited resources Parmeggiani  _et al._
(2003, 2004); Leduc  _et al._ (2012); Subramanian  _et al._ (2013). So far,
most work on the role of limited resources has focused on single components of
the relevant system Adams  _et al._ (2008); Cook and Zia (2009); Cook  _et
al._ (2009); Brackley  _et al._ (2010); Greulich  _et al._ (2012); Ciandrini
_et al._ (2014); Mohapatra  _et al._ (2016, 2017). Only a few studies have
considered simultaneous limitation of two resources Brackley  _et al._ (2012).
In particular, the role of resource limitation has not been explored when two
processes with antagonistic actions are concurrently affected by the limited
availability of protein.

In this Letter, we study the impact of limitations in the supply of both
tubulin and the depolymerizing molecular motor Kip3 on the regulation of MT
length. We build on a recently validated quantitative model of MT dynamics
Melbinger  _et al._ (2012), and extend it to include the constraint of
resource limitations. We find that Kip3 can tightly control MT length,
irrespective of the specific parameter choice. Over a broad range of tubulin
and kinesin concentrations, length regulation is bistable, i.e., the MT can
assume one of two stationary states. We corroborate these findings by
performing _in vitro_ experiments, which show that the MT length distribution
is indeed bimodal for certain concentrations of the components of interest, in
accordance with the theoretical expectations.

(a)

(b)

Figure 1:  Sketch of the model. LABEL:sub@sf_model_schematic A MT in a closed
volume interacts with molecular motors. LABEL:sub@sf_model_lk Motors attach to
the MT lattice at rate $\omega_{\textnormal{A}}$, and detach at rate
$\omega_{\textnormal{D}}$. Motors proceed stepwise toward the plus end at rate
$\nu$, provided the next site is unoccupied. At the tip, motor-induced lattice
depolymerization (rate $\delta$) competes with spontaneous polymerization
(rate $\gamma$). $\omega_{\textnormal{A}}$ and $\gamma$ and depend on the
concentrations of the proteins available in the closed volume, Eqs. (1).

To investigate the impact of limited resources on MT dynamics, we employ a
driven diffusive lattice gas model Krapivsky  _et al._ (2010); Chou  _et al._
(2011); Appert-Rolland  _et al._ (2015) for spontaneous MT polymerization and
kinesin-catalyzed MT depolymerization Melbinger  _et al._ (2012); Reese  _et
al._ (2014), as illustrated in Fig. 1. As kinesin-8 motors predominantly move
along single protofilaments Bormuth  _et al._ (2012); Bugiel  _et al._ (2015),
it suffices to consider a one-dimensional lattice of dynamic length $L(t)$.
The state of each site, $i$, is described by its occupation number, $n_{i}$,
where $n_{i}=0$ and $n_{i}=1$ signify an empty and occupied site,
respectively. On the MT lattice the dynamics follow the totally asymmetric
simple exclusion process with Langmuir kinetics (TASEP/LK) Lipowsky  _et al._
(2001); Klumpp and Lipowsky (2003); Parmeggiani  _et al._ (2003, 2004): Motors
can attach to any empty site on the MT lattice at rate
$\omega_{\textnormal{A}}$, and detach at rate $\omega_{\textnormal{D}}$. Since
binding of motors to the MT depletes the volume concentration of motors
$c_{m}$, the attachment rate $\omega_{\textnormal{A}}$ decreases as

$\displaystyle\omega_{\textnormal{A}}=\omega_{\textnormal{A}}^{0}(c_{m}-m/V)~{}.$
(1a) Here $m$ is the number of motors attached to a protofilament, and $V$ is
the effective volume available to the motors, see Sec. S.III in the Supp. Mat.
Sup . We are specifically interested in the molecular motor Kip3 Varga  _et
al._ (2006); Gupta  _et al._ (2006), which is the kinesin-8 homolog in _S.
cerevisiae_. Based on published _in vitro_ single-molecule experiments, we
estimate its detachment rate to be $\omega_{\textnormal{D}}=4.9\cdot
10^{-3}~{}\textnormal{s}^{-1}$ and the attachment rate to any vacant site as
$\omega_{\textnormal{A}}^{0}=6.7\cdot
10^{-4}~{}\textnormal{nM}^{-1}~{}\textnormal{s}^{-1}$ Varga  _et al._ (2009);
see Sec. S.III. On a protofilament, motors move toward the plus end at rate
$\nu=6.35~{}\textnormal{s}^{-1}$ provided that the next site is empty Varga
_et al._ (2009). At the plus end, Kip3 catalyzes MT shrinkage Walczak  _et
al._ (2013). This is described as a stochastic process where a motor arriving
at the last site removes it at rate $\delta=2.3~{}\textnormal{s}^{-1}$ Varga
_et al._ (2009). At the same time, MTs polymerize spontaneously through
attachment of single tubulin heterodimers to their plus ends. As tubulin
resources are limited, this decreases the volume concentration of tubulin
$c_{T}$ and the polymerization rate,
$\displaystyle\gamma=\gamma^{0}[c_{T}-L/(aV)]~{},$ (1b)

decreases with increasing MT length; here, $a=8.4~{}\textnormal{nm}$ is the
size of a tubulin dimer Hyman  _et al._ (1995), the (net) polymerization rate
per protofilament is
$\gamma_{0}=0.38~{}\mu\textnormal{M}^{-1}\textnormal{s}^{-1}$ Hyman  _et al._
(1992); Brouhard  _et al._ (2008), and the effective volume is $V\approx
1.66~{}\mu\textnormal{m}^{3}$ (Sec. S.III).

(a)

(b)

(c)

(d)

(e)

Figure 2:  Basic phenomenology. Differently colored traces depict different
simulation runs under the indicated starting conditions.
LABEL:sub@sf_trajectories_a–LABEL:sub@sf_trajectories_b MTs evolve toward a
stationary length, which depends on the concentrations of Kip3 and tubulin.
LABEL:sub@sf_densprof_a–LABEL:sub@sf_densprof_b The corresponding motor
density $\rho$ on the MT is shown; see Fig. S2e: $\rho$ increases with
distance from the minus end, and peaks at the plus end.
LABEL:sub@sf_trajectories_bistability The MT length is bistable for a range of
concentrations. Dotted lines: Results of the full mean-field theory.

We performed extensive stochastic simulations Gillespie (1977), and explored
how the MT dynamics depend on the volume concentration of the motor Kip3,
$c_{m}$, and tubulin, $c_{T}$. Figure 2a and 2b show the dynamics of MT length
for two representative concentrations of both components, and for various
initial MT lengths. In all cases, the MT length reaches a stationary state,
albeit at different values. Moreover, the corresponding motor density $\rho$
also differs (Figs. 2c–2d): While for a short stationary length, the overall
motor density is relatively high, it remains low when the MT length is long.
For both cases, the motor density peaks at the plus end (forming a “spike”).

(a)

(b)

(c)

(d)

(e)

Figure 3:  Theoretical results.
LABEL:sub@sf_parameter_scan_short–LABEL:sub@sf_parameter_scan_long _In silico_
scans of the stationary length of MTs, shown in color, as a function of
$c_{m}$ and $c_{T}$. Simulations start from a fully depolymerized lattice
(short) in panel LABEL:sub@sf_parameter_scan_short, in
LABEL:sub@sf_parameter_scan_long the MT is initially fully polymerized (i.e.,
long). In the region bounded by the red lines (obtained from the full MF
theory, Sec. S.II), the stationary length differs for these two cases: Here,
MT dynamics is bistable. LABEL:sub@sf_stability Rate of change of the MT
length, $\partial_{t}L$, as a function of $L$ at
$c_{T}=1.5~{}\mu\textnormal{M}$ for three different motor concentrations, as
obtained from the approximate MF theory. For low and high motor
concentrations, MT length is monostable, while for intermediate
concentrations, two stable stationary states are separated by an unstable
state (bistability). LABEL:sub@sf_comparison_theory_sims Comparison of the
steady-state length obtained from simulations (blue) and the full MF theory
(orange) at $c_{T}=2~{}\mu\textnormal{M}$. LABEL:sub@sf_pd_theory Stability
diagram as obtained from the full MF theory.

We observe that the stationary MT length depends on its initial value for a
certain range of tubulin and kinesin concentrations. Here, depending on
whether the MT starts from a single tubulin dimer or a fully polymerized
filament with all tubulin resources depleted, the stationary length is short
or long, respectively, i.e. the MT dynamics is bistable (Fig. 2e). This is
fundamentally different from MT length regulation with unlimited resources
Melbinger  _et al._ (2012); Reese  _et al._ (2014), where only one stationary
state of finite length is observed. A parameter scan of the tubulin and Kip3
concentrations (Figs. 3a–3b) shows that bistable length regulation occurs over
a broad parameter range.

What physical processes determine MT length and lead to bistability? To answer
this question one needs to analyze the intricate interplay between the
crowding of molecular motors in the lattice bulk Lipowsky  _et al._ (2001);
Parmeggiani  _et al._ (2004) and the (de)polymerization kinetics at the MT tip
Reese  _et al._ (2011), as well as the exchange of resources between filament
and solution. The rate of change of protofilament length is determined by the
antagonism between spontaneous polymerization and Kip3-driven depolymerization
kinetics,

$\displaystyle\partial_{t}L(t)=\bigl{[}\gamma(t)-\rho_{+}(t)\delta\bigr{]}a~{},$
(2a) where $\rho_{+}$ is the probability that the terminal site, i.e. the site
directly at the MT tip, is occupied by a motor. The number of motors on the
protofilament changes when a motor attaches to one of the empty lattice sites,
or any of the motors on it detaches; the number also decreases when a motor
falls off the plus end, taking the last tubulin heterodimer with it. Together,
this yields
$\displaystyle\partial_{t}m(t)\\!=\omega_{\textnormal{A}}(t)\\!\bigl{[}L(t)/a-m(t)\bigr{]}\\!-\omega_{\textnormal{D}}m(t)\\!-\rho_{+}(t)\delta.$
(2b)

In Eqs. (2), the tip density $\rho_{+}$ drives the loss of tubulin dimers and
motors due to depolymerization. This density, in turn, is determined by the
flux of motors along the protofilament toward the MT tip. We assume that these
bulk dynamics are fast in comparison to MT length changes due to
polymerization and depolymerization. Given this time scale separation, the
bulk density can be assumed to be stationary (Sec. S.II A), such that the tip
density is determined by a balance between bulk current and depolymerization
current. Neglecting correlations in the motor density, $\langle
n_{i}n_{j}\rangle\approx\langle n_{i}\rangle\langle n_{j}\rangle$, and
imposing a continuum limit, the mean-field (MF) bulk current is given by
$j(x)=\nu\rho(x)[1-\rho(x)]$, where $\rho(x)$ denotes the average motor
density at position $x$. On length scales of the order of the size of a
tubulin dimer $a$, this current is constant since
$\omega_{\textnormal{A}},\omega_{\textnormal{D}}\ll\nu$, such that the motor
flux in the MT bulk equals the flux off the tip:
$\nu\rho_{L-a}(1-\rho_{L-a})\approx\rho_{+}\delta$. Here, the subscript $L-a$
signifies that the density is evaluated very close to the MT plus end, just
before the density spike begins (cf. Fig. 2c–2d); note that in general
$\rho_{+}\neq\rho_{L-a}$.

In order to determine the bulk density $\rho_{L-a}$, one needs to consider the
combined effects of steric exclusion and motor exchange between filament and
cytosol along the complete MT. In the stationary state, changes in motor
density caused by transport are balanced by attachment-detachment kinetics,
i.e.,

$\nu
a\bigl{(}2\rho-1\bigr{)}\partial_{x}\rho=-\omega_{\textnormal{A}}(1-\rho)+\omega_{\textnormal{D}}\rho~{}.$
(3)

This differential equation has solutions in terms of Lambert $W$-functions
Parmeggiani  _et al._ (2003, 2004) which allow one to compute $\rho_{L-a}$
without any further approximations (Sec. S.II). However, much can already be
learned from an approximate solution, where the density is approximated as a
Taylor series, $\rho(x)\approx Ax+Bx^{2}$; note that $\rho(0)=0$. Upon
inserting this expression into Eq. (3), $A$ and $B$ can be read off by
comparing the coefficients in the ensuing power series, and using
$\rho_{L-a}\approx\rho(L)$. The motor current off the MT, $\rho_{+}\delta$, is
now readily computed, and one obtains to second order in
$\omega_{\textnormal{A,D}}L/a$:

$\rho_{+}\delta\approx\omega_{\textnormal{A}}L/a-(\omega_{\textnormal{A}}+\omega_{\textnormal{D}})(L/a)^{2}\omega_{\textnormal{A}}/(2\nu)~{}.\\\
$ (4)

With Eqs. (2) and (4), we have arrived at a closed set of (nonlinear)
equations for the dynamics of the MT length and the number of motors bound to
a protofilament. It can be viewed as a dynamical system which, as a function
of the control parameters $c_{m}$ and $c_{T}$, may show bifurcations in the
number and nature of its steady states.

The dynamics of nonlinear systems is best visualized by the flow field
$(\partial_{t}m,\partial_{t}L)$ in phase space. Here, the MT state, described
by $L$ and $m$, evolves along the lines drawn in a stream plot (Figs. S17).
This analysis shows that the number of motors bound to the MT equilibrates
almost instantaneously, much more rapidly than the MT length changes.
Therefore, we can assume that the dynamics reduces to the subspace (nullcline)
$\partial_{t}m=0$. This adiabatic elimination of $m$ yields an effective
dynamics of the MT length $L(t)$, as shown in Fig. 3c. Keeping the tubulin
concentration fixed at a typical value of $1.5~{}\mu\textnormal{M}$, we find
that if the motor concentration is either low ($c_{m}=10~{}\textnormal{nM}$)
or high ($c_{m}=30~{}\textnormal{nM}$)), there is only a single state where
the MT length becomes stationary. Hence, regardless of its initial length, a
MT will always reach a uniquely defined stable steady length (monostability).
By contrast, for intermediate motor concentrations
($c_{m}=20~{}\textnormal{nM}$), we observe bistability: Here, three stationary
states exist, two stable states for long and short MT lengths, respectively,
and one unstable state at intermediate MT length, see Fig. 3c. This implies
that, depending on its initial length, a MT may either grow long or remain
short. The same behavior is observed for the full MF analysis, which includes
an exact solution of Eq. (3); see Sec. S.II.

Figure 3d shows that the results obtained from the full MF theory compare very
well with those of the stochastic simulations. In particular, we consistently
observe a bistable regime, with two stable solutions separated by an unstable
solution (separatrix). The stability diagram shown in Fig. 3e summarizes the
different regimes of length regulation as a function of protein
concentrations. In the regimes dominated by depolymerization or
polymerization, the stationary MT length will be short or long, respectively.
At intermediate protein concentrations, the MT length may be short or long
depending on the initial length (bistable regime).

(a)

(b)

(c)

Figure 4:  Experimental results. LABEL:sub@sf_exp_histograms Length
distribution of MTs grown for $3$ hours and subsequently incubated with
various concentrations of Kip3 for one hour. The distribution is unimodal for
$0~{}\textnormal{nM}$, $4~{}\textnormal{nM}$, and $400~{}\textnormal{nM}$
Kip3; it is bimodal for $20~{}\textnormal{nM}$ Kip3, indicating that length
regulation is bistable. LABEL:sub@sf_exp_box Box plots for the MT length show
that the median MT length decreases as the Kip3 concentration is increased
(left axis). The dashed line (right axis) indicates the average number of MTs
per field of view (fov). LABEL:sub@sf_varying_inc_time Length distributions of
MTs initially grown for 3 hours and subsequently incubated with
$20~{}\textnormal{nM}$ Kip3 for various amounts of time. Within an hour, a
bistable distribution is established, and its shape is conserved as the
incubation time is increased.

While these results have been obtained for a single MT, they are not limited
to this case. We find that when many MTs globally access proteins in a well-
mixed pool, length regulation is still accurate, see Sec. S.V. Moreover, the
total length of MTs is bistable in a concentration regime similar to the
single-MT case, Fig. S8. Here, all MTs jointly become short or long, and their
average initial length determines which of these states is reached.

Because diffusion in a real system is fast only on short length scales, and
large systems are not well-mixed, we decided to test the actual behavior
directly by performing a set of _in vitro_ experiments. We grew GMP-CPP
stabilized MTs from a MT polymerization solution based on
$2~{}\mu\textnormal{M}$ tubulin at 27 ∘C (see Sec. S.VI for details). The
resulting MTs had a length distribution similar to a Schulz distribution
Jeune-Smith and Hess (2010), and their median length could be influenced by
varying the incubation time ($1.5$ or $3$ hours). Subsequent to initial MT
polymerization, different Kip3 concentrations were added to samples of the
same polymerization solution and, as a control, no Kip3 was added to the final
aliquot; all parts were incubated for another hour, so that MT polymerization
from the remaining free tubulin and Kip3-induced depolymerization could occur
simultaneously. The resultant MT length distributions in the samples were
imaged as described in Sec. S.VI.

In the first experiment, MTs grown for $3$ hours were incubated with 0, 4, 20,
and 400 nM Kip3 for another hour. In the absence of Kip3 (0 nM), the length
distribution of MTs peaked around $11~{}\mu\textnormal{m}$ (Fig. 4a and Fig.
S20). The presence of Kip3 reduced the median MT length (Fig. 4b, box
plots/left axis), and also decreased the number of MTs per field of view (Fig.
4b, dashed line/right axis). The latter indicates that a number of MTs were
completely depolymerized or shrank below the detection limit of our setup.
Notably, at low and high Kip3 concentrations of $4~{}\textnormal{nM}$ and
$400~{}\textnormal{nM}$, the length distributions were unimodal with peaks
around $11~{}\mu\textnormal{m}$ and $2.5~{}\mu\textnormal{m}$ (Fig. 4a),
respectively. This corresponds to the monostable regimes at low and high motor
concentrations derived above, where polymerization and depolymerization,
respectively, dominate and the final length is independent of the initial
length.

In contrast, at a Kip3 concentration of $20~{}\textnormal{nM}$, the MT length
distribution was qualitatively different: Here, two distinct populations of
MTs (peaks around $2.5~{}\mu\textnormal{m}$ and $11~{}\mu\textnormal{m}$) were
observed, resulting in a bimodal length distribution (Fig. 4a). We could
exclude that the short MTs observed in this experiment were additionally
nucleated after the addition of Kip3, Fig. S19. Furthermore, the two peaks are
not transient. In contrast, a bimodal distribution is fully established within
an hour, and the qualitative distribution remained intact until the end of our
experiment, Fig. 4c. The bimodal length distribution must therefore result
from bistable length regulation: According to Fig. S7, MTs in a well-mixed
many-filament system will jointly become short or long in the bistable regime,
and their local average initial length distinguishes between these cases.
However, diffusion of protein is fast only on short length scales, and slow in
large systems, such as our experimental setup. In addition, the association of
motors with MTs, on which they may remain for minutes or longer, significantly
slows down diffusion in crowded environments, e.g., inside cells. The
resulting separation of length scales of the small well-mixed range and the
large system size may hence allow different regions of a system to develop
independently. Given a broad initial distribution of MT lengths, the local
average length of MTs in some regions is therefore in the domain of attraction
of the steady state with long length, while in other regions MTs are attracted
towards the short length. Hence, MTs in distant spatial regions evolve towards
the different fixed points and domains with long and short filaments are
formed, which coexist at stationarity. This interpretation is supported by the
length distribution of MTs resulting from a solution of Kip3 and tubulin which
is incubated for 1 hour in a shaker at the same conditions otherwise, Fig.
S18. Because constant mixing leads to a global well-mixed reservoir, the
resulting length distribution is unimodal, confirming our expectations.

We then sought to obtain further information about the domains of attraction
of the respective stationary states and the corresponding separatrix marking
the boundary between these domains (Fig. S17b). If the MT length distribution
at which the length regulation process starts is short, MTs in all regions
will be in the domain of attraction of the short stationary length. To test
this prediction, we stopped MT growth after $1.5$ hours and subsequently added
the same amounts of Kip3 to the polymerization solution as before. The median
MT length in the absence of Kip3 was significantly shorter (Fig. S15b) than
the corresponding value for MTs grown for $3$ hours. We observed that the
length distribution remained unimodal when Kip3 was added, irrespective of its
concentration, Fig. S15a. This indicates that, after 1.5 hours of initial MT
polymerization, filament lengths still lie below the separatrix in Fig. S15b.
Taken together, our experimental findings qualitatively confirm our
theoretical predictions, including the existence of a regime where MT length
regulation by Kip3 gives rise to two populations of filaments with clearly
distinct lengths.

Taking a broader perspective, we believe that – similar to the case considered
here – effects of resource limitation are of relevance to other aspects of
mitotic spindle formation and disassembly, and other processes in which
protein availability in the cytosol constrains dynamic interactions.

###### Acknowledgements.

E.F. and S.D. acknowledge support by the German Excellence Initiative via the
programs “NanoSystems Initiative Munich” (NIM) and “Center for Advancing
Electronics Dresden” (cfaed), respectively. L.R. acknowledges support by the
Netherlands Organisation for Scientific Research (NWO), via the researc
program “Membrane-cytoplasm protein shuttling”, Project No. 13PMP03.

## References

  * Wilbur and Heald (2013) J. D. Wilbur and R. Heald, eLife 2, e00290 (2013).
  * Borisy and Taylor (1967) G. G. Borisy and E. W. Taylor, Journal of Cell Biology 34, 525 (1967).
  * Borisy (1967) G. G. Borisy, The Journal of Cell Biology 34, 535 (1967).
  * Good  _et al._ (2013) M. C. Good, M. D. Vahey, A. Skandarajah, D. A. Fletcher, and R. Heald, Science 342, 856 (2013).
  * Hazel  _et al._ (2013) J. Hazel, K. Krutkramelis, P. Mooney, M. Tomschik, K. Gerow, J. Oakey, and J. C. Gatlin, Science 342, 853 (2013).
  * Winey  _et al._ (1995) M. Winey, C. L. Mamay, E. T. O’Toole, D. N. Mastronarde, T. H. Giddings, K. L. McDonald, and J. R. McIntosh, The Journal of Cell Biology 129, 1601 (1995).
  * Howard and Hyman (2003) J. Howard and A. A. Hyman, Nature 422, 753 (2003).
  * Howard and Hyman (2007) J. Howard and A. A. Hyman, Current Opinion in Cell Biology 19, 31 (2007).
  * Varga  _et al._ (2006) V. Varga, J. Helenius, K. Tanaka, A. A. Hyman, T. U. Tanaka, and J. Howard, Nature Cell Biology 8, 957 (2006).
  * Goshima  _et al._ (2005) G. Goshima, R. Wollman, N. Stuurman, J. M. Scholey, and R. D. Vale, Current Biology 15, 1979 (2005).
  * Stumpff  _et al._ (2008) J. Stumpff, G. von Dassow, M. Wagenbach, C. Asbury, and L. Wordeman, Developmental Cell 14, 252 (2008).
  * Varga  _et al._ (2009) V. Varga, C. Leduc, V. Bormuth, S. Diez, and J. Howard, Cell 138, 1174 (2009).
  * Reese  _et al._ (2011) L. Reese, A. Melbinger, and E. Frey, Biophysical Journal 101, 2190 (2011).
  * Melbinger  _et al._ (2012) A. Melbinger, L. Reese, and E. Frey, Physical Review Letters 108, 258104 (2012).
  * Reese  _et al._ (2014) L. Reese, A. Melbinger, and E. Frey, Interface Focus 4, 1 (2014).
  * Lipowsky  _et al._ (2001) R. Lipowsky, S. Klumpp, and T. M. Nieuwenhuizen, Physical Review Letters 87, 108101 (2001).
  * Cook and Zia (2009) L. J. Cook and R. K. P. Zia, Journal of Statistical Mechanics: Theory and Experiment 2009, P02012 (2009).
  * Greulich  _et al._ (2012) P. Greulich, L. Ciandrini, R. J. Allen, and M. C. Romano, Physical Review E 85, 011142 (2012).
  * Ciandrini  _et al._ (2014) L. Ciandrini, I. Neri, J. C. Walter, O. Dauloudet, and A. Parmeggiani, Physical Biology 11, 056006 (2014).
  * Parmeggiani  _et al._ (2003) A. Parmeggiani, T. Franosch, and E. Frey, Physical Review Letters 90, 086601 (2003).
  * Parmeggiani  _et al._ (2004) A. Parmeggiani, T. Franosch, and E. Frey, Physical Review E 70, 046101 (2004).
  * Leduc  _et al._ (2012) C. Leduc, K. Padberg-Gehle, V. Varga, D. Helbing, S. Diez, and J. Howard, Proceedings of the National Academy of Sciences 109, 6100 (2012).
  * Subramanian  _et al._ (2013) R. Subramanian, S.-C. Ti, L. Tan, S. A. Darst, and T. M. Kapoor, Cell 154, 377 (2013).
  * Adams  _et al._ (2008) D. A. Adams, B. Schmittmann, and R. K. P. Zia, Journal of Statistical Mechanics: Theory and Experiment 2008, P06009 (2008).
  * Cook  _et al._ (2009) L. J. Cook, R. K. P. Zia, and B. Schmittmann, Physical Review E 80, 031142 (2009).
  * Brackley  _et al._ (2010) C. A. Brackley, M. C. Romano, and M. Thiel, Physical Review E 82, 051920 (2010).
  * Mohapatra  _et al._ (2016) L. Mohapatra, B. L. Goode, P. Jelenkovic, R. Phillips, and J. Kondev, Annu. Rev. Biophys. 45, 85 (2016).
  * Mohapatra  _et al._ (2017) L. Mohapatra, T. J. Lagny, D. Harbage, P. R. Jelenkovic, and J. Kondev, Cell Systems 4, 559 (2017).
  * Brackley  _et al._ (2012) C. A. Brackley, L. Ciandrini, and M. C. Romano, Journal of Statistical Mechanics: Theory and Experiment 2012, P03002 (2012).
  * Krapivsky  _et al._ (2010) P. L. Krapivsky, S. Redner, and E. Ben-Naim,  _A Kinetic View of Statistical Physics_ (Cambridge University Press, 2010).
  * Chou  _et al._ (2011) T. Chou, K. Mallick, and R. K. P. Zia, Reports on Progress in Physics 74, 116601 (2011).
  * Appert-Rolland  _et al._ (2015) C. Appert-Rolland, M. Ebbinghaus, and L. Santen, Physics Reports 593, 1 (2015), arXiv:1507.06166 .
  * Bormuth  _et al._ (2012) V. Bormuth, B. Nitzsche, F. Ruhnow, A. Mitra, M. Storch, B. Rammner, J. Howard, and S. Diez, Biophysical Journal 103, L4 (2012).
  * Bugiel  _et al._ (2015) M. Bugiel, E. Böhl, and E. Schäffer, Biophysical Journal 108, 2019 (2015).
  * Klumpp and Lipowsky (2003) S. Klumpp and R. Lipowsky, Journal of Statistical Physics 113, 233 (2003).
  * (36) See Supplementary Materials which includes Refs. Krug (1991); Kolomeisky  _et al._ (1998); Frey  _et al._ (2004); Strogatz (2014); Gardner  _et al._ (2011); Telley  _et al._ (2009); Schneider  _et al._ (2015); Graf and Frey (2017); Halatek and Frey (2012, 2018); Salmon  _et al._ (1984); Castoldi and Popov (2003); Korten  _et al._ (2011); Ruhnow  _et al._ (2011) .
  * Gupta  _et al._ (2006) M. L. Gupta, P. Carvalho, D. M. Roof, and D. Pellman, Nature Cell Biology 8, 913 (2006).
  * Walczak  _et al._ (2013) C. E. Walczak, S. Gayek, and R. Ohi, Annual Review of Cell and Developmental Biology 29, 417 (2013).
  * Hyman  _et al._ (1995) A. A. Hyman, D. Chrétien, I. Arnal, and R. H. Wade, The Journal of cell biology 128, 117 (1995).
  * Hyman  _et al._ (1992) A. A. Hyman, S. Salser, D. N. Drechsel, N. Unwin, and T. J. Mitchison, Molecular Biology of the Cell 3, 1155 (1992).
  * Brouhard  _et al._ (2008) G. J. Brouhard, J. H. Stear, T. L. Noetzel, J. Al-Bassam, K. Kinoshita, S. C. Harrison, J. Howard, and A. A. Hyman, Cell 132, 79 (2008).
  * Gillespie (1977) D. T. Gillespie, The Journal of Physical Chemistry 81, 2340 (1977).
  * Jeune-Smith and Hess (2010) Y. Jeune-Smith and H. Hess, Soft Matter 6, 1778 (2010).
  * Krug (1991) J. Krug, Physical Review Letters 67, 1882 (1991).
  * Kolomeisky  _et al._ (1998) A. B. Kolomeisky, G. M. Schütz, E. B. Kolomeisky, and J. P. Straley, Journal of Physics A: Mathematical and General 31, 6911 (1998).
  * Frey  _et al._ (2004) E. Frey, A. Parmeggiani, and T. Franosch, Genome informatics. International Conference on Genome Informatics 15, 46 (2004).
  * Strogatz (2014) S. H. Strogatz,  _Nonlinear Dynamics and Chaos_ (Westview Press, 2014).
  * Gardner  _et al._ (2011) M. K. Gardner, B. D. Charlebois, I. M. Jánosi, J. Howard, A. J. Hunt, and D. J. Odde, Cell 146, 582 (2011).
  * Telley  _et al._ (2009) I. A. Telley, P. Bieling, and T. Surrey, Biophysical Journal 96, 3341 (2009).
  * Schneider  _et al._ (2015) R. Schneider, T. Korten, W. J. Walter, and S. Diez, Biophysical Journal 108, 2249 (2015).
  * Graf and Frey (2017) I. R. Graf and E. Frey, Physical Review Letters 118, 128101 (2017).
  * Halatek and Frey (2012) J. Halatek and E. Frey, Cell Reports 1, 741 (2012).
  * Halatek and Frey (2018) J. Halatek and E. Frey, Nat. Phys. (2018), 10.1038/s41567-017-0040-5.
  * Salmon  _et al._ (1984) E. D. Salmon, W. M. Saxton, R. J. Leslie, M. L. Karow, and J. R. McIntosh, Journal of Cell Biology 99, 2157 (1984).
  * Castoldi and Popov (2003) M. Castoldi and A. V. Popov, Protein Expression and Purification 32, 83 (2003).
  * Korten  _et al._ (2011) T. Korten, B. Nitzsche, C. Gell, F. Ruhnow, C. Leduc, and S. Diez, “Fluorescence Imaging of Single Kinesin Motors on Immobilized Microtubules,” in  _Single Molecule Analysis: Methods and Protocols_, edited by E. J. G. Peterman and G. J. L. Wuite (Humana Press, Totowa, NJ, 2011) pp. 121–137.
  * Ruhnow  _et al._ (2011) F. Ruhnow, D. Zwicker, and S. Diez, Biophysical Journal 100, 2820 (2011).