# Stable High Power deep-UV Enhancement Cavity in Ultra High Vacuum with
Fluoride Coatings

Zakary Burkley  Institute for Particle Physics and Astrophysics, ETH Zürich,
CH-8093 Zürich, Switzerland  Lucas de Sousa Borges  Institute for Particle
Physics and Astrophysics, ETH Zürich, CH-8093 Zürich, Switzerland  Ben Ohayon
Institute for Particle Physics and Astrophysics, ETH Zürich, CH-8093 Zürich,
Switzerland  Artem Golovozin  Lebedev Physical Institute, Moscow, Russia
Jesse Zhang  Institute for Particle Physics and Astrophysics, ETH Zürich,
CH-8093 Zürich, Switzerland  Paolo Crivelli  Institute for Particle Physics
and Astrophysics, ETH Zürich, CH-8093 Zürich, Switzerland

###### Abstract

We demonstrate the superior performance of fluoride coated versus oxide coated
mirrors in long term vacuum operation of a high power deep-ultraviolet
enhancement cavity. In high vacuum (${10}^{-8}\text{\,}\mathrm{mbar}$), the
fluoride optics can maintain up to a record-high $10\text{\,}\mathrm{W}$ of
stable intracavity power on one hour time scales, whereas for the oxide
optics, we observe rapid degradation at lower intracavity powers with a rate
that increases with power. After observing degradation in high vacuum, we can
recover the fluoride and oxide optics with oxygen; however, this recovery
process becomes ineffective after several applications. For fluoride coatings,
we see that initial UV conditioning in an oxygen environment helps to improve
the performances of the optics. In oxygen-rich environments from
$\sim$${10}^{-4}\text{\,}\mathrm{mbar}1\text{\,}\mathrm{mbar}$, the fluoride
optics can stably maintain up to $20\text{\,}\mathrm{W}$ of intracavity power
on several-hour time scales whereas for the oxide optics there is immediate
degradation with a rate that increases with decreasing oxygen pressure.

##  1 Introduction

Continual innovation of deep-ultraviolet (deep-UV) laser technology is
beneficial for precision scientific measurement, materials processing,
biomedicine, and many other industrial and medicinal applications [1, 2, 3].
Understanding the performance of optics at high average deep-UV powers is of
growing importance. Both picosecond pulsed and continuous wave laser systems
below $250\text{\,}\mathrm{nm}$ operating at Watt-level average powers have
been demonstrated recently [4, 5, 6]. Numerous applications would immediately
gain from increased stability at these powers such as laser cooling of AlF at
$227.5\text{\,}\mathrm{nm}$ [7] and mercury at $253.7\text{\,}\mathrm{nm}$
[8], as well as two-photon cooling of hydrogen at $243\text{\,}\mathrm{nm}$
[9]. A laser-based parity measurement in hydrogen [10] would also benefit from
increased deep-UV Rabi frequencies, as would Rydberg excitation for fast
entanglement of strontium ions [11] or studies of strong long-range
interactions in a lithium optical lattice [12] at $243\text{\,}\mathrm{nm}$ or
$231\text{\,}\mathrm{nm}$, respectively. Anti-matter wave interferometry with
positronium could also be envisaged [13].

One particular area of research that profits notably from advances in deep-UV
radiation sources is precision two-photon spectroscopy of simple and exotic
atoms. For example, comparison of the 1S-2S transition in hydrogen [14] and
anti-hydrogen [15] has tested CPT invariance at a level of
$2\text{\times}{10}^{-12}$, and requires $243\text{\,}\mathrm{nm}$ radiation.
At $205\text{\,}\mathrm{nm}$, differing measurements on the 1S-3S transition
of hydrogen [16, 17] has contributed notably to the proton radius puzzle. A
new experiment underway on the 1S-2S transition in Muonium at
$244\text{\,}\mathrm{nm}$ promises a substantial improvement in the
determination of the muon mass and tests of lepton universality [18]; at
$252.5\text{\,}\mathrm{nm}$, precision two-photon spectroscopy in Xenon [19]
has recently been demonstrated for future neutron EDM measurements at TRIUMF
[20].

While the high precision advantage of two-photon spectroscopy comes from its
removal of the first-order Doppler effect and access to long-lived states,
this second order process suffers from low cross-sections. This can be
overcome with the high peak power offered through pulsed laser systems, but
introduces many systematics—these can be removed through use of a continuous
wave (CW) laser enhanced in an optical cavity. While such enhancement cavities
have been built, e.g., [14, 15, 19], they operate at modest intracavity powers
($<1\text{\,}\mathrm{W}$ CW), and are limited to low buildup by the available
reflectivity of deep-UV optics. Impressively, up to $30\text{\,}\mathrm{W}$ of
CW intracavity power has been demonstrated, but required
$0.67\text{\,}\mathrm{mbar}$ of oxygen on the mirrors to prevent UV-assisted
degradation [9]. This degradation is attributed to surface oxygen depletion
and/or hydrocarbon contamination [21, 22, 23]. In principle, both of these
effects are reversible by keeping O2 present on the mirrors, as was required
in [9, 17, 19]. However, such a solution is not ideal. In [19], the precision
of the measurement was limited by pressure shifts introduced by the required
O2 level. In [9], reaching a vacuum of
$5\text{\times}{10}^{-8}\text{\,}\mathrm{mbar}$ for spectroscopy required a
technically demanding differential pumping setup with geometric constraints
that significantly limited the intracavity power [24]. Furthermore,
implementing such pumping requirements would be challenging in the setups
required to study anti-hydrogen [15] and muonium [18] where complicated
beamline and cryogenic environment constraints exist. For these reasons, a
stable deep-UV cavity operating with Watt-level circulating powers in high
vacuum or reduced oxygen pressures would be of great interest to the laser
spectroscopy community.

The previous enhancement cavities that observed degradation utilized oxide
coatings for their high reflective mirrors. To better understand the
degradation process, an interesting alternative are fluoride coatings. With no
oxides present, the latter should be insensitive to surface oxygen depletion.
Furthermore, low fluence pulsed studies of fluoride coatings show an extremely
high durability with initial UV conditioning, an effect not observed with
oxide coatings [25, 26]. Therefore, in this work, we test the performance of
both fluoride-based and oxide-based coatings at intra-cavity powers from
$1\text{\,}\mathrm{W}10\text{\,}\mathrm{W}$ at $244\text{\,}\mathrm{nm}$ in
high vacuum, and $15\text{\,}\mathrm{W}20\text{\,}\mathrm{W}$ in various
oxygen pressures. At high powers, we observe superior performance in the
stability of the fluoride coated optics relative to the oxide coated optics.

Figure 1: Experimental setup: $244\text{\,}\mathrm{nm}$ radiation is generated
through a frequency quadrupled ytterbium-doped fiber amplified laser system.
This radiation is then coupled into a vacuum enhancement cavity. Extended
cavity diode laser (ECDL), tapered amplifier (TA), second harmonic generation
(SHG), lithium triborate (LBO), polarizing beam splitter (PBS), electro-optic
modulator (EOM), L1-L5 (mode matching lenses), photodiode (PD), power monitor
(PM), caesium lithium borate (CLBO), input coupler (IC), output coupler (OC),
piezoelectric transducer stack (PZT).

##  2 Experimental setup

Inspired by recent demonstrations of high-power Yb-fiber amplifiers [27, 28,
6] near $976\text{\,}\mathrm{nm}$, the system utilizes a frequency doubled
fiber-amplifier developed in collaboration with Toptica for generation of the
fundamental IR radiation at $976\text{\,}\mathrm{nm}$ and second harmonic
generation to $488\text{\,}\mathrm{nm}$. The system can deliver up to
$7\text{\,}\mathrm{W}$ of IR power and $5\text{\,}\mathrm{W}$ in the blue;
with this, up to $1.8\text{\,}\mathrm{W}$ of $244\text{\,}\mathrm{nm}$ power
can be generated in a ”home-built” resonant doubling stage.

The second harmonic enhancement cavity that generates 244 nm radiation uses a
standard bow tie geometry and is maintained on resonance using the Pound-
Drever-Hall (PDH) locking technique [29]. A slow piezoelectric (PZT) stack
(CTS NAC2123-H20-C01) and fast PZT chip (CTS NAC2122-C04) are used for active
length stabilization, with the latter mounted using the technique described in
[30]. The cavity mirrors (LaserOptik Gmbh) consist of a $2.5\%$ transmission
input coupler and three HR mirrors at $488\text{\,}\mathrm{nm}$ ($<99.8\%$);
two of these mirrors have a radius of curvature of 200 mm to focus the beam
within the doubling crystal with a waist of
$47\text{\,}\mathrm{\SIUnitSymbolMicro m}$ by
$44\text{\,}\mathrm{\SIUnitSymbolMicro m}$ . The doubling occurs in a
$10\text{\,}\mathrm{mm}$ long Brewster-cut Caesium Lithium Borate (CLBO)
crystal (Oxide) in a Type-1 critical phase matching configuration
($\theta=76.4^{\circ},\phi=45^{\circ}$). The hygroscopic crystal is kept at
$150\text{\,}\mathrm{\SIUnitSymbolCelsius}$ in an Al oven and flushed
continuously with oxygen. A dichroic Brewster plate (Spectral Optics) couples
the UV radiation out of the cavity.

As was shown in [6], the CLBO can output $>1\text{\,}\mathrm{W}$ of radiation
on the one hour time scale with no evidence of degradation. Enabled by the
increased stability in the IR and blue stages of our laser system, we have
observed similar powers on significantly longer time scales (See Fig. 2). We
measure a slight degradation of $\sim$1-2$\%$ per hour in total power, which
we attribute to damage on the $244\text{\,}\mathrm{nm}$ coating of the output
coupler. This damage is visible by eye upon inspection, and the power can be
recovered by moving to a new spot on the output coupler.

Figure 2: Performance of the $244\text{\,}\mathrm{nm}$ laser system over long
time scales. The laser can output $1\text{\,}\mathrm{W}$ of
$244\text{\,}\mathrm{nm}$ power over several hour time scales without
degradation of the CLBO crystal. We attribute the slow degradation observed to
damage of the output coupler.

The $244\text{\,}\mathrm{nm}$ light is enhanced in a linear, vacuum cavity
composed of two calcium fluoride (CaF2) substrates coated with fluoride
through electron beam evaporation (LaserOptik Gmbh). The input coupler (IC) is
flat with a measured transmission of $0.024$, and the curved (radius of
curvature = $4\text{\,}\mathrm{m}$) high reflector (HR) has a manufacturer
specified reflectivity $>98.5\%$. The mirrors are $0.50\text{\,}\mathrm{m}$
apart, and the Gaussian beam inside the cavity is approximately collimated
with a beam waist of $\sim 0.3\text{\,}\mathrm{mm}$. One spherical lens and
two cylindrical lenses (to compensate for walk-off in the CLBO) are used for
mode matching into cavity. When optimally aligned, measurements of the
rejection signal ratio indicate input coupling efficiency of $80\%$. We
isolate the back-reflected power using a polarizing beam-splitter and a
Faraday rotation stage. A portion of this back-reflected signal is used for
PDH locking, where the HR is mounted on an annular PZT stack (CTS
NAC2123-H20-C01) in the method described in [31] for active stabilization.

The transmitted power through the HR is used to monitor the intracavity power
with a calibrated UV thermal power meter (Thorlabs S401C). This, combined with
the measured HR transmission ($T_{\text{HR}}=0.0012$), yields the intracavity
power, $P_{cav}$. As seen in Fig. 1, the input power, $P_{o}$ is monitored
with a sensitive, calibrated UV photodiode (Thorlabs S120VC) placed directly
after the last mirror before the vacuum enhancement cavity. The mirror is
backside polished, and the leakage light is calibrated to the UV power
entering the cavity. In this way, both the input power and intracavity power
are monitored non-destructively, enabling in situ measurement of the cavity
power enhancement $\beta=P_{\text{\it{cav}}}/P_{\text{\it{o}}}$ in a way which
is independent of any power drifts from the UV laser system leading to the
cavity.

The cavity is operated in a variable oxygen purge mode or ultra high vacuum
mode. In the O2 mode, two leak valves are adjusted to maintain a fixed flow of
fresh O2 (6.0 grade) on the cavity mirrors. For high O2 pressures of
${10}^{-1}\text{\,}\mathrm{mbar}$ or more, only a roughing pump is used. For
O2 pressures at or below ${10}^{-3}\text{\,}\mathrm{mbar}$, a turbo pump is
used to maintain the pressure equilibrium. The lowest equilibrium pressure of
fresh O2 we can confidently maintain is ${10}^{-4}\text{\,}\mathrm{mbar}$,
limited by the precision of the leak valves. In the UHV mode, there is no
oxygen purge, and the pressures reach below ${10}^{-8}\text{\,}\mathrm{mbar}$.
To limit possible hydrocarbon contaminants, we minimize use of any plastics
and only use UHV compatible components. The vacuum windows are sealed with
indium instead of O-rings, the PZT is UHV compatible, the HR is mounted to the
PZT with Torr seal, and the PZT is mounted following the design in [31] using
PEEK screws, as PEEK is known to have low outgassing rates. All other elements
in the system are either stainless steel, copper, or aluminum.

##  3 Results and discussion

We measured the power enhancement factor $\beta$ versus time for both a SiO2
and CaF2 input coupler in high vacuum and different oxygen pressures, with
various intracavity powers. The high reflector is CaF2 for all measurements.

In Figure 3, we give the enhancement factor in high vacuum as a function of
time for both a SiO2 and CaF2 input couplers for different intracavity powers.
At low power of $\sim 1.5\text{\,}\mathrm{W}$ (Fig. 3(a)), we observe little
or no degradation for both coatings. However, there appears to be slow
oscillations of small amplitude, with the SiO2 input coupler. As we increase
the intracavity powers by inserting more input power (Figs. 3(b)-3(d)), we
observe an increased rate of degradation with the SiO2 IC, limiting the
asymptotic power to a few watts. Small amplitude, slow oscillations are also
observed in these measurements, and appear to relax on a shorter time scale
for increased starting powers. At higher powers, the enhancement factor
decreases by approximately a factor of two in one hour, with a rate that
decreases as the overall enhancement factor decreases. Similar to what was
observed in [9], degradation is easily visible on the coating after
inspection.

(a) $1.5\text{\,}\mathrm{W}$ intra-cavity power for both input couplers.

(b) $3\text{\,}\mathrm{W}$ intra-cavity power for both input couplers.

(c) $5\text{\,}\mathrm{W}$ intra-cavity power for CaF2 IC,
$5.5\text{\,}\mathrm{W}$ start power for SiO2 IC.

(d) $10\text{\,}\mathrm{W}$ intra-cavity power for CaF2 IC,
$9\text{\,}\mathrm{W}$ start power for SiO2 IC.

Figure 3: Enhancement versus time of CaF2 IC versus SiO2 IC for different
intracavity powers at ${10}^{-8}\text{\,}\mathrm{mbar}$.

For CaF2, the data given in Figs. 3(a)-3(d) was measured continuously while
increasing the input power every hour. Furthermore, prior to the high vacuum
measurement, the CaF2 mirror was conditioned in
$1\text{\times}{10}^{-3}\text{\,}\mathrm{mbar}$ of O2 with
$16\text{\,}\mathrm{W}$ of intracavity power for several hours. As seen in
Figs. 3(a)-3(d), performance with the CaF2 IC is stable on one hour timescales
with up to $10\text{\,}\mathrm{W}$ of intracavity power. A slight decrease in
the enhancement from 40 to 34 is visible. We attribute the drop in enhancement
to a slow degradation of the CaF2 input coupler. Rotating the IC is sufficient
to regain enhancement, and upon careful inspection, an inhomogeneity in the
coating is visible by eye after running at these high powers in high vacuum
for an extended period. Nonetheless, the drop in enhancement of 4$\%$ per hour
seen with this IC is an order of magnitude smaller than the 50$\%$ drop seen
with the oxide-coated input coupler. As shown in Fig. 4, we observed that the
rate of enhancement decrease for CaF2 was highly reliant on the initial UV
conditioning of the mirrors in an oxygen environment.

Figure 4: Enhancement versus time for unconditioned and conditioned CaF2 IC.
Intracavity power of $8\text{\,}\mathrm{W}$ for conditioned IC,
$10\text{\,}\mathrm{W}$ start power for unconditioned IC.

We next tested whether higher average intra-cavity powers could be
accommodated, relying on the higher start enhancement of the SiO2 coated IC by
introducing intermittent oxygen revivals. As shown in Fig. 5, in accordance
with [23, 9], the SiO2 optics can be rapidly revived after degradation by
flowing oxygen while maintaining a buildup of ultraviolet light in the
enhancement cavity. However, even though reviving is possible a few times,
eventually the enhancement can no longer be fully recovered. A similar
behavior was observed with the CaF2 input coupler.

Figure 5: Oxygen recovery of enhancement at ${10}^{-8}\text{\,}\mathrm{mbar}$
with $3\text{\,}\mathrm{W}$ start power for SiO2 IC, by momentarily flushing
O2 on the mirrors surface.

As observed, an intermittent revival mode with oxygen in high vacuum is not a
sustainable option. However, as oxygen revives the input couplers, it is
interesting to study their performance in oxygen-rich environments. The
performance of the mirrors at high intracavity powers with a continuous oxygen
purge of ${10}^{-4}\text{\,}\mathrm{mbar}1\text{\,}\mathrm{mbar}$ is shown in
Fig. 6. With SiO2, we observe a rapid initial drop in enhancement, that
increases with decreasing oxygen pressure. However, with CaF2, we see stable
powers over the entire range of oxygen pressures.

(a) $1\text{\,}\mathrm{mbar}$ of O2.
$15\text{\,}\mathrm{W}$ of intra-cavity power for CaF2 IC,
$17\text{\,}\mathrm{W}$ start power for SiO2 IC.

(b) ${10}^{-1}\text{\,}\mathrm{mbar}$ of O2.
$15\text{\,}\mathrm{W}$ intra-cavity power for CaF2 IC,
$12\text{\,}\mathrm{W}$ start power SiO2 IC.

(c) ${10}^{-3}\text{\,}\mathrm{mbar}$ of O2.
$18\text{\,}\mathrm{W}$ intra-cavity power for CaF2 IC,
$13\text{\,}\mathrm{W}$ intra-cavity start power for SiO2 IC.

(d) ${10}^{-4}\text{\,}\mathrm{mbar}$ of O2.
$18\text{\,}\mathrm{W}$ intra-cavity power for CaF2 IC,
$12\text{\,}\mathrm{W}$ start power for SiO2 IC.

Figure 6: Enhancement versus time of CaF2 IC versus SiO2 IC at different O2
pressures.

As was the case with Figs. 3(a)-3(d), Figs. 6(a)-6(d) were measured
continuously for CaF2. Unlike the high vacuum case, we do not observe a
noticeable drop in enhancement over this four hour time scale, indicating
improved performance in oxygen. To test this, we moved to a new spot on the
CaF2 IC, and measured for four hours at ${10}^{-3}\text{\,}\mathrm{mbar}$ as
seen in Fig. 7. After an initial slight decrease in enhancement, we observe
degradation-free performance over this time scale. We attribute the initial
drop in the first 30 minutes to thermal heating of the optics before the
enhancement cavity, which can affect the input coupling efficiency.

Figure 7: ${10}^{-3}\text{\,}\mathrm{mbar}$ of O2, $16\text{\,}\mathrm{W}$
intra-cavity power for CaF2 IC.

To our knowledge, CW intracavity powers between
$5\text{\,}\mathrm{W}30\text{\,}\mathrm{W}$ at these deep-UV wavelengths have
been only demonstrated once before [9]. They used SiO2 mirrors and observed
rapid ($<1\text{\,}\mathrm{min}$) degradation of the mirrors at oxygen
pressures below $0.67\text{\,}\mathrm{mbar}$, attributing the result to
hydrocarbon contamination versus surface oxygen depletion. Using the CaF2
input coupler, our system cannot be affected by surface oxygen depletion.
Therefore, it is reasonable to attribute the degradation of the CaF2 optics to
hydrocarbon contamination. As seen in Fig. 6(a), we observe degradation at
$1\text{\,}\mathrm{mbar}$ for the same SiO2 input coupler used in [9].

UV-ozone cleaning is a commonly used method to remove surface hydrocarbon
contaminants [32]. With an O2 purge, the 244 nm intracavity radiation
continuously generates ozone and atomic oxygen. Simultaneously, the 244 nm
radiation also decomposes hydrocarbon contaminants into a variety of
substances, including excited molecules and free radicals. The latter react
with the atomic oxygen to form simpler volatile molecules, such as H2O, CO2,
and NOx, that are easily desorbed from the surface. The robustness of UV-ozone
cleaning in our system is demonstrated clearly by the revival of the SiO2
input coupler in Fig. 5.

Interestingly, as seen in Figs. 6(b)-6(d) and 7, with the CaF2 input coupler,
we observe degradation-free enhancement on one hour time scales with up to
four orders less oxygen maintained on the mirrors than was required for SiO2
in [9], and on four hour time scales with up to three orders of magnitude less
oxygen. This suggests that the UV-ozone cleaning method is more effective for
fluoride optics or there is an additional mechanism involved. UV conditioning
of fluoride coatings is a well documented method to decrease their absorption
[25, 33, 26], and was clearly observed in our results (see Fig. 4). While the
exact conditioning mechanism is not fully understood, Heber et. al. [33, 26]
attribute it partially to radiation desorption of hydrocarbons on the surface,
which they could not observe with oxide coatings. This desorption was observed
in both air and argon atmospheres, suggesting hydrocarbon contaminants on
fluoride coatings can be removed in oxygen-free or oxygen reduced
environments—removing reliance on UV-ozone cleaning. This would explain why
the CaF2 input coupler in our system can operate degradation free at
significantly reduced O2 pressures. However, as observed in high vacuum, the
radiation desorption cleaning method by itself is unable to fully prevent
hydrocarbon contamination. This suggests the degradation free performance at
low oxygen pressures relies on both cleaning methods. What remains unclear is
the inability of UV-ozone cleaning to completely remove the effects
hydrocarbon contamination in high vacuum over time, as seen in Fig. 5. Future
work should try to characterize the surface contaminants present after this
irreversible degradation, or whether the coating is damaged.

##  4 Summary

In this work we have shown that relative to oxide-coated optics, a high-power,
deep-UV CW enhancement cavity operating at $244\text{\,}\mathrm{nm}$ in high
vacuum with UV conditioned fluoride-coated optics, can have a degradation rate
an order of magnitude smaller. This enables stable operation on hour time
scales with up to $10\text{\,}\mathrm{W}$ of intracavity power. To our
knowledge, these are the highest CW deep-UV intracavity powers maintained on
these time scales in high vacuum, and almost an order of magnitude higher than
has been demonstrated with oxide-coated optics at these wavelengths. While
combining oxygen with the deep-UV radiation can revive the optics after
degradation, this method of revival is not sustainable.

However, with a continuous oxygen purge, we have shown that fluoride-coated
optics can operate degradation-free with intracavity powers of up to
$20\text{\,}\mathrm{W}$ on multi-hour time scales with three to four orders of
magnitude less oxygen than is required to prevent severe degradation with
oxide-coated optics. Demonstration of such high powers on these time scales
was enabled by the stability of our CW $244\text{\,}\mathrm{nm}$ laser system,
which can maintain Watt-level output powers on half-day time scales.

Degradation of optics is a common hindrance to long term stability and power
scaling of systems requiring deep-UV wavelengths. We have shown that fluoride
optics are a promising alternative to oxide optics for development of a robust
deep-UV CW enhancement cavity at increased intracavity powers and reduced
reliance on UV-ozone cleaning.

Funding ERC consolidator grant (818053-Mu-MASS) and Swiss National Foundation
grant 197346.

Acknowledgements The authors would like to thank Dylan Yost (CSU) for
providing some of the optics used here as well as thoughtful discussions about
the research. We also appreciate Samuel Cooper’s (CSU) comments on this
article as well as Michele Masseroni and Rebekka Garreis at ETH for their time
helping us clean and diagnose the mirrors. We are extremely grateful to
Konstantinos Simeonidis for his contributions to the laser system.

Disclosures The authors declare no conflicts of interest.

## References

  * [1] Neil Savage.  Ultraviolet lasers.  Nature Photonics, 1(2):83–85, 2007.
  * [2] Haoyuan Lu, Hao Xu, Jianye Zhao, and Dong Hou.  A deep ultraviolet mode-locked laser based on a neural network.  Scientific reports, 10(1):1–6, 2020.
  * [3] Callum R Smith, Asbjørn Moltke, Abubakar I Adamu, Mattia Michieletto, Patrick Bowen, Peter M Moselund, Christos Markos, and Ole Bang.  Low-noise tunable deep-ultraviolet supercontinuum laser.  Scientific Reports, 10(1):1–11, 2020.
  * [4] Benjamin Willenberg, Fabian Brunner, Christopher R Phillips, and Ursula Keller.  High-power picosecond deep-uv source via group velocity matched frequency conversion.  Optica, 7(5):485–491, 2020.
  * [5] Hana Turcicova, Ondrej Novak, Lukas Roskot, Martin Smrz, Jiri Muzik, Michal Chyla, Akira Endo, and Tomas Mocek.  New observations on duv radiation at 257 nm and 206 nm produced by a picosecond diode pumped thin-disk laser.  Optics express, 27(17):24286–24299, 2019.
  * [6] Z. Burkley, A. D. Brandt, C. Rasor, S. F. Cooper, and D. C. Yost.  Highly coherent, watt-level deep-uv radiation via a frequency-quadrupled yb-fiber laser system.  Appl. Opt., 58(7):1657–1661, Mar 2019.
  * [7] Stefan Truppe, Silvio Marx, Sebastian Kray, Maximilian Doppelbauer, Simon Hofsäss, Hanns Christian Schewe, Nicole Walter, Jesús Pérez-Ríos, Boris G Sartakov, and Gerard Meijer.  Spectroscopic characterization of aluminum monofluoride with relevance to laser cooling and trapping.  Physical Review A, 100(5):052513, 2019.
  * [8] Ruchen Zhao, Xiaohu Fu, Lei Zhang, Su Fang, Jianfang Sun, Yan Feng, Zhen Xu, and Yuzhu Wang.  High-power continuous-wave narrow-linewidth 253.7 nm deep-ultraviolet laser.  Applied optics, 56(32):8973–8977, 2017.
  * [9] SF Cooper, Z Burkley, AD Brandt, C Rasor, and DC Yost.  Cavity-enhanced deep ultraviolet laser for two-photon cooling of atomic hydrogen.  Optics letters, 43(6):1375–1378, 2018.
  * [10] C Rasor and DC Yost.  Laser-based measurement of parity violation in hydrogen.  Physical Review A, 102(3):032801, 2020.
  * [11] Chi Zhang, Fabian Pokorny, Weibin Li, Gerard Higgins, Andreas Pöschl, Igor Lesanovsky, and Markus Hennrich.  Submicrosecond entangling gate between trapped ions via rydberg interaction.  Nature, 580(7803):345–349, 2020.
  * [12] Elmer Guardado-Sanchez, Benjamin Spar, Peter Schauss, Ron Belyansky, Jeremy T Young, Przemyslaw Bienias, Alexey V Gorshkov, Thomas Iadecola, and Waseem S Bakr.  Quench dynamics of a fermi gas with strong long-range interactions.  arXiv preprint arXiv:2010.05871, 2020.
  * [13] M. K. Oberthaler.  Antimatter wave interferometry with positronium.  Nucl. Instrum. Meth. B, 192:129–134, 2002.
  * [14] Christian G Parthey, Arthur Matveev, Janis Alnis, Birgitta Bernhardt, Axel Beyer, Ronald Holzwarth, Aliaksei Maistrou, Randolf Pohl, Katharina Predehl, Thomas Udem, et al.  Improved measurement of the hydrogen 1 s–2 s transition frequency.  Physical review letters, 107(20):203001, 2011.
  * [15] M Ahmadi, BXR Alves, CJ Baker, W Bertsche, A Capra, C Carruth, CL Cesar, M Charlton, S Cohen, R Collister, et al.  Characterization of the 1s–2s transition in antihydrogen.  Nature, 557(7703):71–75, 2018.
  * [16] Hélene Fleurbaey, Sandrine Galtier, Simon Thomas, Marie Bonnaud, Lucile Julien, François Biraben, François Nez, Michel Abgrall, and Jocelyne Guéna.  New measurement of the 1 s- 3 s transition frequency of hydrogen: contribution to the proton charge radius puzzle.  Physical review letters, 120(18):183001, 2018.
  * [17] Alexey Grinin, Arthur Matveev, Dylan C Yost, Lothar Maisenbacher, Vitaly Wirthl, Randolf Pohl, Theodor W Hänsch, and Thomas Udem.  Two-photon frequency comb spectroscopy of atomic hydrogen.  Science, 370(6520):1061–1066, 2020.
  * [18] Paolo Crivelli.  The mu-mass (muonium laser spectroscopy) experiment.  Hyperfine Interactions, 239(1):1–9, 2018.
  * [19] Emily Altiere, Eric R Miller, Tomohiro Hayamizu, David J Jones, Kirk W Madison, and Takamasa Momose.  High-resolution two-photon spectroscopy of a transition of xenon.  Physical Review A, 97(1):012507, 2018.
  * [20] Ruediger Picker.  How the minuscule can contribute to the big picture: the neutron electric dipole moment project at triumf.  In Proceedings of the 14th International Conference on Meson-Nucleon Physics and the Structure of the Nucleon (MENU2016), page 010005, 2017.
  * [21] Roderick R Kunz, Vladimir Liberman, and Deanna K Downs.  Experimentation and modeling of organic photocontamination on lithographic optics.  Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures Processing, Measurement, and Phenomena, 18(3):1306–1313, 2000.
  * [22] Jeromy Hollenshead and Leonard Klebanoff.  Modeling radiation-induced carbon contamination of extreme ultraviolet optics.  Journal of Vacuum Science & Technology B: Microelectronics and Nanometer Structures Processing, Measurement, and Phenomena, 24(1):64–82, 2006\.
  * [23] Dorian Gangloff, Molu Shi, Tailin Wu, Alexei Bylinskii, Boris Braverman, Michael Gutierrez, Rosanna Nichols, Junru Li, Kai Aichholz, Marko Cetina, et al.  Preventing and reversing vacuum-induced optical losses in high-finesse tantalum (v) oxide mirror coatings.  Optics express, 23(14):18014–18028, 2015.
  * [24] Samuel F Cooper, Adam D Brandt, Cory Rasor, Zakary Burkley, and Dylan C Yost.  Cryogenic atomic hydrogen beam apparatus with velocity characterization.  Review of Scientific Instruments, 91(1):013201, 2020.
  * [25] E Eva, K Mann, N Kaiser, B Anton, R Henking, D Ristau, P Weissbrodt, D Mademann, L Raupach, and E Hacker.  Laser conditioning of laf 3/mgf 2 dielectric coatings at 248 nm.  Applied optics, 35(28):5613–5619, 1996.
  * [26] Joerg Heber, Roland Thielsch, Holger Blaschke, Norbert Kaiser, Uwe Leinhos, and A Görtler.  Changes in optical interference coatings exposed to 193-nm excimer laser radiation.  In Laser-Induced Damage in Optical Materials: 1998, volume 3578, pages 83–96. International Society for Optics and Photonics, 1999.
  * [27] Z Burkley, C Rasor, SF Cooper, AD Brandt, and DC Yost.  Yb fiber amplifier at 972.5 nm with frequency quadrupling to 243.1 nm.  Applied Physics B, 123(1):1–6, 2017.
  * [28] Jingwei Wu, Xiushan Zhu, Hua Wei, Kort Wiersma, Michael Li, Jie Zong, Arturo Chavez-Pirson, Valery Temyanko, L. J. LaComb, R. A. Norwood, and N. Peyghambarian.  Power scalable 10 w 976 nm single-frequency linearly polarized laser source.  Opt. Lett., 43(4):951–954, Feb 2018.
  * [29] RWP Drever, John L Hall, FV Kowalski, J_ Hough, GM Ford, AJ Munley, and H Ward.  Laser phase and frequency stabilization using an optical resonator.  Applied Physics B, 31(2):97–105, 1983.
  * [30] Travis C Briles, Dylan C Yost, Arman Cingöz, Jun Ye, and Thomas R Schibli.  Simple piezoelectric-actuated mirror with 180 khz servo bandwidth.  Optics express, 18(10):9739–9746, 2010.
  * [31] A. Chadi, G. Méjean, R. Grilli, and D. Romanini.  Note: Simple and compact piezoelectric mirror actuator with 100 khz bandwidth, using standard components.  Review of Scientific Instruments, 84(5):056112, 2013.
  * [32] Rajiv Kohli.  Chapter 9 - applications of uv-ozone cleaning technique for removal of surface contaminants.  In Rajiv Kohli and K.L. Mittal, editors, Developments in Surface Contamination and Cleaning: Applications of Cleaning Techniques, pages 355–390. Elsevier, 2019.
  * [33] Joerg Heber, Roland Thielsch, Holger Blaschke, Norbert Kaiser, Klaus R Mann, Eric Eva, Uwe Leinhos, and Andreas Goertler.  Stability of optical interference coatings exposed to low-fluence 193-nm arf radiation.  In Optical Microlithography XI, volume 3334, pages 1041–1047. International Society for Optics and Photonics, 1998.


# Global asymptotic stability of an active disassembly model of flagellar
length control

Thomas G. Fai  Youngmin Park

###### Abstract

Organelle size control is a fundamental question in biology that demonstrates
the fascinating ability of cells to maintain homeostasis within their highly
variable environments. Theoretical models describing cellular dynamics have
the potential to help elucidate the principles underlying size control. Here,
we perform a detailed mathematical study of the active disassembly model
proposed in [Fai et al, Length regulation of multiple flagella that self-
assemble from a shared pool of components, _eLife_ , 8, (2019): e42599]. We
construct a Lyapunov function to show that the dynamical equations for the
flagellar lengths are well-behaved and rule out the possibility of
oscillations in the model. We prove global asymptotic stability in the case of
two flagella and show that this result generalizes to arbitrary flagellar
number. Our analysis reveals a defect in the original active disassembly
model, which is manifested as corners in the phase space that yield unphysical
solutions. To address this defect, we construct a hybrid system. Although this
hybrid system leads to behaviors such as time delays and the accumulation of
trajectories in the corners of phase space, we show that these scenarios do
not typically arise in applications.

Key words: length control; flagella; mathematical modeling; dynamical systems

Robust size control is essential in biology at the levels of organs, cells,
and organelles. The absence of size regulation is a marker of pathologies such
as cancer [1], and the varied sizes of biological structures highlight
differences in the fundamental physical constraints on biological systems [2,
3, 4].

At the level of organelles, a particular useful model organism for flagellar
length control is the unicellular green algae _Chlamydomonas reinhardtii_
(Figure 1(a))[5, 6, 7, 8, 9]. _Chlamydomonas_ uses its two flagella to propel
itself through its aqueous environment and the lengths of its flagella are
important for maintaining swimming speed and direction [10, 11].

Experiments on _Chlamydomonas_ show that flagellar maintenance and length
control is a dynamic process. During severing experiments in which a single
flagellum is amputated, there is an initial phase in which the two flagellar
lengths equalize to a shorter intermediate length and a subsequent phase in
which the two flagella grow together to the initial steady-state length [12]
(shown schematically in Figure 1(b)). This reveals the coupling between the
length dynamics of the two flagella.

Figure 1: _Chlamydomonas reinhardtii_ and the active disassembly model of
flagellar length control. (a) Experimental image of _Chlamydomonas_ from [7]
provided by the author. The two flagella have equal lengths of approximately
12 µm. (b) Upon severing one flagellum, the severed flagellum begins to
regenerate while the unsevered flagellum shortens. The flagella rapidly
equalize at an intermediate length before returning to their initial lengths
on a slower timescale. (c) Schematic of the active disassembly model [13], in
which the concentration gradient generated by a diffusing depolymerase leads
to a length-dependent rate of disassembly at the flagellar tip.

The balance-point model, in which the steady-state length is determined by the
competition between assembly and disassembly processes, has been proposed to
account for the length regulation observed in severing experiments [14, 15].
Here, we focus on a particular balance point model suitable for organisms such
as _Chlamydomonas reinhardtii_ , in which multiple flagella are assembled
simultaneously from a shared pool of building blocks. In the original balance-
point model proposed for _Chlamydomonas_ , the assembly rate is taken to be a
monotonically decreasing function of length whereas the disassembly rate is
taken to be length-independent [14, 15]. This yields a unique point of
intersection that determines the steady-state length.

We have shown previously [13] that it is important to consider not only the
predictions of the balance-point model for a single flagellum, but also its
predictions for the coupled behavior of two (or more) flagella in light of
severing experiments. Whereas constant disassembly models are unable to
account for the length-equalization observed in severing experiments, we have
proposed an active disassembly model based on physical first principles of
flux balance (Figure 1(c)) in which length-dependent disassembly leads to
simultaneous length control.

In this work, we consider the mathematical properties of the active
disassembly model as a system of nonlinear ordinary differential equations. We
show that these equations have a global asymptotically stable steady-state
solution. To prove this, we find a Lyapunov function for the dynamics. The
Lyapunov function also implies there are no limit cycles or periodic
oscillations. However, there is still the possibility of nontrivial dynamics.
We construct a hybrid system to ensure that the flagellar lengths remain non-
negative, and this yields regions of state space from which all trajectories
contract to a single point on the boundary before flowing toward steady-state.
We interpret these trajectories in light of existing experimental data. In
addition, we provide alternate proofs of global asymptotic stability and the
non-existence of periodic oscillations that do not rely on the construction of
a Lyapunov function. We discuss how these results generalize to other
situations.

By characterizing the mathematical properties of the active disassembly model,
we aim to provide a strong foundation for subsequent applications and analysis
of the model. In the Discussion, we note several open questions regarding the
mathematical model and its application to flagellar length control in
_Chlamydomonas_.

##  1 Active disassembly model of flagellar length control

Flagella are dynamic structures. During the process known as intraflagellar
transport (IFT), protein assemblies known as IFT trains are transported by
motor molecules, some that carry trains toward the tip (anterograde direction)
and others that carry IFT trains toward the base of the flagella (retrograde
direction) [16]. In the balance-point model, this turnover is captured by the
continual assembly and disassembly processes which reach a dynamic steady-
state.

Let $L_{1}(t)$ and $L_{2}(t)$ denote the lengths of two flagella on a single
_Chlamydomonas_ organism at a particular time $t$. The flagellum is composed
of tubulin filaments, and we assume there is a total amount $T$ of tubulin
available.

Similarly, we assume there are a total number $M$ of anterograde molecular
motors, of which a number $M_{f}$ are present in the pool at the flagellar
base while the rest are undergoing IFT. The motors move ballistically in the
anterograde direction with velocity and diffusively in the retrograde
direction with diffusion coefficient $D$ [17]. The injection flux $J$ of IFT
particles into the flagella satisfies mass action kinetics, so that

$J=k_{\text{on}}M_{f}.$ (1)

There is a separation of timescales between the motion of IFT particles, which
traverse the flagellum in several seconds, and the flagellar length dynamics,
which varies over a timescale of several minutes [18]. As shown in [13], upon
combining (1) with conservation of protein number, this separation of
timescales leads to a quasi-steady state approximation for the flux:

$J=\frac{k_{\text{on}}M/2}{1+k_{\text{on}}(L_{1}+L_{2})/2v+k_{\text{on}}(L_{1}^{2}+L_{2}^{2})/4D}.$
(2)

The flux $J$ represents the number of IFT particles injected per second into
the flagella, and these trains are assumed to carry an amount of tubulin
proportional to the free tubulin in the basal pool. Therefore, the overall
speed of assembly is equal to:

$\text{assembly}=\gamma J\left(T-L_{1}-L_{2}\right),$ (3)

where $\gamma$ is a proportionality constant that represents the fraction of
tubulin arriving at the flagellar tip that is incorporated into the flagellum.
Disassembly is taken to be a linear function of the depolymerase concentration
at the flagellar tip, i.e. each depolymerizing enzyme located near the
flagellar tip removes a constant number of tubulin subunits per second. This
implies that for the $i$th flagellum, the speed of disassembly satisfies

$\text{disassembly}_{i}=d_{0}+d_{1}\frac{JL_{i}}{D},$ (4)

for $i=1,\,2$, assuming that the depolymerase has the same pattern of motion
as the anterograde motors [13]. Setting the growth rate to be the difference
between assembly and disassembly speeds leads to a system of coupled nonlinear
differential equations:

$\displaystyle L_{1}^{\prime}$ $\displaystyle=\gamma
J(T-L_{1}-L_{2})-d_{0}-d_{1}\frac{JL_{1}}{D}$ (5) $\displaystyle
L_{2}^{\prime}$ $\displaystyle=\gamma
J(T-L_{1}-L_{2})-d_{0}-d_{1}\frac{JL_{2}}{D}.$ (6)

where $J\equiv J(L_{1},L_{2})$ satisfies

$J(L_{1},L_{2})=\frac{k_{\text{on}}M/2}{1+k_{\text{on}}(L_{1}+L_{2})/2v+k_{\text{on}}(L_{1}^{2}+L_{2}^{2})/4D}.$
(7)

###  1.1 Nondimensionalization

We next rewrite eqs. 5 and 6 in dimensionless form. In terms of the
nondimensional length $\widetilde{L}=L/T$ and nondimensional time
$\widetilde{t}=t\gamma k_{\text{on}}M/2$ as well as the dimensionless
parameters $\pi_{0}=2d_{0}/\gamma k_{\text{on}}MT$, $\pi_{1}=d_{1}/\gamma D$,
$\pi_{2}=k_{\text{on}}T/2v$, and $\pi_{3}=k_{\text{on}}T^{2}/4D$, the system
of equations eqs. 5 and 6 becomes

$\displaystyle\frac{\mathrm{d}\widetilde{L}_{1}}{\mathrm{d}\widetilde{t}}=\widetilde{J}\left(1-\widetilde{L}_{1}-\widetilde{L}_{2}\right)-\pi_{0}-\pi_{1}\widetilde{J}\widetilde{L}_{1}$
(8)
$\displaystyle\frac{\mathrm{d}\widetilde{L}_{2}}{\mathrm{d}\widetilde{t}}=\widetilde{J}\left(1-\widetilde{L}_{1}-\widetilde{L}_{2}\right)-\pi_{0}-\pi_{1}\widetilde{J}\widetilde{L}_{2},$
(9)

where the dimensionless flux $\widetilde{J}=2J/k_{\text{on}}M$ satisfies

$\widetilde{J}=\frac{1}{1+\pi_{2}(\widetilde{L}_{1}+\widetilde{L}_{2})+\pi_{3}(\widetilde{L}_{1}^{2}+\widetilde{L}_{2}^{2})}.$
(10)

We interpret the dimensionless parameters as follows: $\pi_{0}$ is the ratio
of disassembly and assembly rates, $\pi_{1}$ is the depolymerase activity
level, $\pi_{2}=\tau_{b}/\tau_{i}$ is the ratio of the ballistic timescale
$\tau_{b}:=T/v$ of IFT transport to the injection timescale
$\tau_{i}:=k_{\text{on}}^{-1}$, and $\pi_{3}=\tau_{d}/\tau_{i}$ is an
analogous ratio of the diffusive timescale $\tau_{d}=T^{2}/2D$ to the
injection timescale. (We could equivalently think of $\pi_{2}$ and $\pi_{3}$
as ratios of lengthscales related to the same physical processes.)

It is a straightforward calculation to verify that both flagella have the same
steady-state length $\widetilde{L}_{\text{ss}}$ given by

$\widetilde{L}_{\text{ss}}=\frac{1+\pi_{0}\pi_{2}+\pi_{1}/2}{2\pi_{0}\pi_{3}}\left(\sqrt{1+\frac{2\left(1-\pi_{0}\right)\pi_{0}\pi_{3}}{\left(1+\pi_{0}\pi_{2}+\pi_{1}/2\right)^{2}}}-1\right).$
(11)

##  2 Nonlinear behavior

To simplify notation, we drop tildes and use the dimensionless formulation for
the remainder of the article. The active disassembly model eqs. 8 and 9 is not
a gradient system, i.e. it cannot be written as the gradient of an energy
function $E(L_{1},L_{2})$. This is proven by contradiction. Suppose there is
an energy function such that $L_{1}^{\prime}=\partial E/\partial L_{1}$ and
$L_{2}^{\prime}=\partial E/\partial L_{2}$. Then, by equality of mixed
partials

$\frac{\partial^{2}E}{\partial L_{1}\partial L_{2}}=\frac{\partial
L_{1}^{\prime}}{\partial L_{2}}=\frac{\partial L_{2}^{\prime}}{\partial
L_{1}}=\frac{\partial^{2}E}{\partial L_{2}\partial L_{1}},$ (12)

i.e.

$\frac{\partial\left(J\left(1-L_{1}-L_{2}\right)-\pi_{0}-\pi_{1}JL_{1}\right)}{\partial
L_{2}}=\frac{\partial\left(J\left(1-L_{1}-L_{2}\right)-\pi_{0}-\pi_{1}JL_{2}\right)}{\partial
L_{1}}.$ (13)

However, the left-hand side of the the inequality is equal to

$-J+\left(1-L_{1}-L_{2}\right)\frac{\partial J}{\partial
L_{2}}-\pi_{1}L_{1}\frac{\partial J}{\partial L_{2}},$ (14)

whereas the right-hand side is equal to

$-J+\left(1-L_{1}-L_{2}\right)\frac{\partial J}{\partial
L_{1}}-\pi_{1}L_{2}\frac{\partial J}{\partial L_{1}}.$ (15)

Because $J(L_{1},L_{2})=J(L_{2},L_{1})$, we have $\partial J/\partial
L_{1}=\partial J/\partial L_{2}$. Assuming $\pi_{1}\neq 0$, follows that there
is equality if and only if $L_{1}=L_{2}$. We have thereby reached the
contradiction since $L_{1}\neq L_{2}$ in general. It follows that no energy
function exists.

However, using a technique applicable to many nonlinear dynamical systems, we
may construct a Lyapunov function to show that the steady state is globally
asymptotically stable. First, we specify the subspace $\Omega$ with
$(L_{1},L_{2})\in\Omega\subset\mathbb{R}^{2}$ over which the model is defined.
Because we must have $L_{1},\,L_{2}\geq 0$ for the variables to have physical
meaning as lengths, $\Omega$ is contained within the first quadrant. Moreover,
the sum of lengths cannot exceed the total pool size so that $L_{1}+L_{2}\leq
1$. To summarize, $\Omega$ is defined to be the triangular region in the first
quadrant bound by the lines $L_{1}=0$, $L_{2}=0$, and $L_{1}+L_{2}\leq 1$.

We define the deviations from steady-state $\Delta L_{1}:=L_{1}-L_{\text{ss}}$
and $\Delta L_{2}:=L_{2}-L_{\text{ss}}$ on $\Omega$ and express the dynamics
in terms of $\Delta L_{1}$ and $\Delta L_{2}$:

$\displaystyle\frac{\mathrm{d}L_{1}}{\mathrm{d}t}$
$\displaystyle=J\left(1-2L_{\text{ss}}-\Delta L_{1}-\Delta
L_{2}\right)-\pi_{0}-\pi_{1}J\left(L_{\text{ss}}+\Delta L_{1}\right)$ (16)
$\displaystyle\frac{\mathrm{d}L_{2}}{\mathrm{d}t}$
$\displaystyle=J\left(1-2L_{\text{ss}}-\Delta L_{1}-\Delta
L_{2}\right)-\pi_{0}-\pi_{1}J\left(L_{\text{ss}}+\Delta L_{2}\right),$ (17)

or equivalently

$\displaystyle\frac{\mathrm{d}L_{1}}{\mathrm{d}t}$
$\displaystyle=J\left(1-2L_{\text{ss}}\right)-J\left(\Delta L_{1}+\Delta
L_{2}\right)-\pi_{0}-\pi_{1}JL_{\text{ss}}-\pi_{1}J\Delta L_{1}$ (18)
$\displaystyle\frac{\mathrm{d}L_{2}}{\mathrm{d}t}$
$\displaystyle=J\left(1-2L_{\text{ss}}\right)-J\left(\Delta L_{1}+\Delta
L_{2}\right)-\pi_{0}-\pi_{1}JL_{\text{ss}}-\pi_{1}J\Delta L_{2}.$ (19)

By symmetry, the steady-state lengths $L_{1,\text{ss}}$ and $L_{2,\text{ss}}$
are equal. We denote the common steady-state length by $L_{\text{ss}}$. It is
determined by setting the derivative equal to zero in (18) or (19), which
yields the identity

$0=J_{\text{ss}}\left(1-2L_{\text{ss}}\right)-\pi_{0}-\pi_{1}J_{\text{ss}}L_{\text{ss}},$
(20)

where

$J_{\text{ss}}=\frac{1}{1+2\pi_{2}L_{\text{ss}}+2\pi_{3}L_{\text{ss}}^{2}}.$
(21)

(Eq. (20) yields a quadratic equation for $L_{\text{ss}}$, which has a unique
positive root given by (11).)

We next use the steady-state relationship (20) to rewrite eqs. 18 and 19 in a
different form. Before we proceed, we first note the identities given by

$\displaystyle J\left(1-2L_{\text{ss}}\right)$
$\displaystyle=(J-J_{\text{ss}})\left(1-2L_{\text{ss}}\right)+J_{\text{ss}}\left(1-2L_{\text{ss}}\right)$
(22) $\displaystyle\pi_{1}JL_{\text{ss}}$
$\displaystyle=\pi_{1}(J-J_{\text{ss}})L_{\text{ss}}+\pi_{1}J_{\text{ss}}L_{\text{ss}},$
(23)

where the first line follows from adding and subtracting
$J_{\text{ss}}\left(1-2L_{\text{ss}}\right)$ and the second line follows from
adding and subtracting $\pi_{1}J_{\text{ss}}L_{\text{ss}}$. Plugging the above
expressions into eqs. 18 and 19 and using the steady-state relation (20)
yields

$\displaystyle\frac{\mathrm{d}L_{1}}{\mathrm{d}t}$
$\displaystyle=\left(J-J_{\text{ss}}\right)\left(1-2L_{\text{ss}}\right)-J\left(\Delta
L_{1}+\Delta
L_{2}\right)-\pi_{1}\left(J-J_{\text{ss}}\right)L_{\text{ss}}-\pi_{1}J\Delta
L_{1}$ (24) $\displaystyle\frac{\mathrm{d}L_{2}}{\mathrm{d}t}$
$\displaystyle=\left(J-J_{\text{ss}}\right)\left(1-2L_{\text{ss}}\right)-J\left(\Delta
L_{1}+\Delta
L_{2}\right)-\pi_{1}\left(J-J_{\text{ss}}\right)L_{\text{ss}}-\pi_{1}J\Delta
L_{2}.$ (25)

Adding and subtracting the above equations yields

$\displaystyle\frac{\mathrm{d}L_{1}}{\mathrm{d}t}+\frac{\mathrm{d}L_{2}}{\mathrm{d}t}$
$\displaystyle=2\left(J-J_{\text{ss}}\right)\left(1-2L_{\text{ss}}-\pi_{1}L_{\text{ss}}\right)-\left(2+\pi_{1}\right)J\left(\Delta
L_{1}+\Delta L_{2}\right)$ (26)
$\displaystyle\frac{\mathrm{d}L_{1}}{\mathrm{d}t}-\frac{\mathrm{d}L_{2}}{\mathrm{d}t}$
$\displaystyle=-\pi_{1}J\left(\Delta L_{1}-\Delta L_{2}\right)$ (27)

It follows that

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\Delta
L_{1}+\Delta L_{2}\right)^{2}$
$\displaystyle=2\left(J-J_{\text{ss}}\right)\left(1-2L_{\text{ss}}-\pi_{1}L_{\text{ss}}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)-\left(2+\pi_{1}\right)J\left(\Delta L_{1}+\Delta
L_{2}\right)^{2}$ (28)
$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\Delta
L_{1}-\Delta L_{2}\right)^{2}$ $\displaystyle=-\pi_{1}J\left(\Delta
L_{1}-\Delta L_{2}\right)^{2}.$ (29)

where we have used
$\mathrm{d}L_{i}/\mathrm{d}t=\mathrm{d}(L_{\text{ss}}+\Delta
L_{i})/\mathrm{d}t=\mathrm{d}\Delta L_{i}/\Delta t$ for $i=1,\,2$ and

$\displaystyle\frac{1}{2}\frac{\mathrm{d}(\Delta L_{1}+\Delta
L_{2})^{2}}{\mathrm{d}t}$ $\displaystyle=(\Delta L_{1}+\Delta
L_{2})\frac{\mathrm{d}(\Delta L_{1}+\Delta L_{2})}{\mathrm{d}t}$ (30)
$\displaystyle\frac{1}{2}\frac{\mathrm{d}(\Delta L_{1}-\Delta
L_{2})^{2}}{\mathrm{d}t}$ $\displaystyle=(\Delta L_{1}-\Delta
L_{2})\frac{\mathrm{d}(\Delta L_{1}-\Delta L_{2})}{\mathrm{d}t}.$ (31)

We use eqs. 28 and 29 to form a Lyapunov function for the dynamics. Let

$\phi(L_{1},L_{2})=\frac{1}{2}\left[\left(\Delta L_{1}+\Delta
L_{2}\right)^{2}+\frac{\pi_{3}}{\pi_{1}}\left(\Delta L_{1}-\Delta
L_{2}\right)^{2}\right].$ (32)

######  Claim.

$\phi$ is a Lyapunov function on $\Omega$.

We first comment on some implications of the above claim before proceeding
with its proof. It follows from the claim that there is a unique steady-state
$L_{\text{ss}}$ at which
$(L_{1},L_{2})=(L_{\text{ss}},L_{\text{ss}})\in\Omega$ and that no limit
cycles exist in $\Omega$.

###### Proof.

Since it is the sum of squared terms with positive factors, $\phi\geq 0$. It
remains to show that $\phi^{\prime}\leq 0$ with equality achieved only at the
unique steady-state $L_{1}=L_{2}=L_{\text{ss}}$. Because the second term of
$\phi$ has a non-positive derivative (by (29)), a sufficient condition for
$\phi^{\prime}\leq 0$ is that the right-hand side of (28) is negative away
from the steady-state, i.e.

$2\left(J-J_{\text{ss}}\right)\left(1-2L_{\text{ss}}-\pi_{1}L_{\text{ss}}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)-\left(2+\pi_{1}\right)J\left(\Delta L_{1}+\Delta
L_{2}\right)^{2}<0,$ (33)

except when $\Delta L_{1}=\Delta L_{2}=0$ in which case it is zero. As we show
next, in some cases (33) holds and the claim immediately follows. Later on, we
show how, if (33) does not hold, we may use the second term of $\phi$ to
obtain the desired result $\phi^{\prime}\leq 0$.

It follows from the steady-state condition (20) that

$1-2L_{\text{ss}}-\pi_{1}L_{\text{ss}}=\frac{\pi_{0}}{J_{\text{ss}}}.$ (34)

Therefore, the first term on the left-hand side of (33) may be rewritten as

$\displaystyle
2\left(J-J_{\text{ss}}\right)\left(1-2L_{\text{ss}}-\pi_{1}L_{\text{ss}}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)$
$\displaystyle=\frac{2\pi_{0}}{J_{\text{ss}}}\left(J-J_{\text{ss}}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)$
$\displaystyle=2\pi_{0}J_{\text{ss}}^{-1}\left(\frac{1}{J^{-1}}-\frac{1}{J_{\text{ss}}^{-1}}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)$
$\displaystyle=2\pi_{0}J_{\text{ss}}^{-1}\left(\frac{J_{\text{ss}}^{-1}-J^{-1}}{J_{\text{ss}}^{-1}J^{-1}}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)$
$\displaystyle=2\pi_{0}J\left(J_{\text{ss}}^{-1}-J^{-1}\right)\left(\Delta
L_{1}+\Delta L_{2}\right).$ (35)

To simplify the right-hand side of the above expression, we use the definition
of the flux (10) to get

$\displaystyle J_{\text{ss}}^{-1}-J^{-1}$
$\displaystyle=1+2\pi_{2}L_{\text{ss}}+2\pi_{3}L_{\text{ss}}^{2}-\left(1+\pi_{2}\left(L_{1}+L_{2}\right)+\pi_{3}\left(L_{1}^{2}+L_{2}^{2}\right)\right)$
$\displaystyle=-\pi_{2}\left(\Delta L_{1}+\Delta
L_{2}\right)-\pi_{3}\left(L_{\text{ss}}^{2}-L_{1}^{2}+L_{\text{ss}}^{2}-L_{2}^{2}\right).$
(36)

Note that

$\displaystyle L_{1}^{2}-L_{\text{ss}}^{2}+L_{2}^{2}-L_{\text{ss}}^{2}$
$\displaystyle=(L_{1}+L_{\text{ss}})\Delta L_{1}+(L_{2}+L_{\text{ss}})\Delta
L_{2}$ $\displaystyle=L_{1}\Delta L_{1}+L_{2}\Delta L_{2}+L_{\text{ss}}(\Delta
L_{1}+\Delta L_{2}),$ (37)

and moreover

$\displaystyle L_{1}\Delta L_{1}+L_{2}\Delta L_{2}$
$\displaystyle=\frac{L_{1}+L_{2}}{2}\left(\Delta L_{1}+\Delta
L_{2}\right)+\left(L_{1}-\frac{L_{1}+L_{2}}{2}\right)\Delta
L_{1}+\left(L_{2}-\frac{L_{1}+L_{2}}{2}\right)\Delta L_{2}$
$\displaystyle=\frac{L_{1}+L_{2}}{2}\left(\Delta L_{1}+\Delta
L_{2}\right)+\frac{L_{1}-L_{2}}{2}\Delta L_{1}-\frac{L_{1}-L_{2}}{2}\Delta
L_{2}$ $\displaystyle=\frac{L_{1}+L_{2}}{2}\left(\Delta L_{1}+\Delta
L_{2}\right)+\frac{1}{2}\left(\Delta L_{1}-\Delta L_{2}\right)^{2}.$ (38)

Combining eqs. 35, 36, 37 and 38, for the left-hand side of (33) we have

$\displaystyle 2\pi_{0}J\left(J_{\text{ss}}^{-1}-J^{-1}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)$ $\displaystyle=-2\pi_{0}J\left(\pi_{2}\left(\Delta
L_{1}+\Delta L_{2}\right)^{2}+\pi_{3}L_{\text{ss}}\left(\Delta L_{1}+\Delta
L_{2}\right)^{2}\right.$
$\displaystyle\left.+\frac{\pi_{3}}{2}\left(L_{1}+L_{2}\right)\left(\Delta
L_{1}+\Delta L_{2}\right)^{2}+\frac{\pi_{3}}{2}\left(\Delta L_{1}-\Delta
L_{2}\right)^{2}\left(\Delta L_{1}+\Delta L_{2}\right)\right).$ (39)

To show the above expression is negative, we consider two cases. In the case
$\Delta L_{1}+\Delta L_{2}>0$, the term in parentheses on the right-hand side
of (39) is a sum of squared terms multiplied by positive numbers. Therefore
the right-hand side of (39) is negative.

In the case $\Delta L_{1}+\Delta L_{2}<0$, we may bound the positive term by
the negative terms in the right-hand side of (39) to show that the overall
expression is negative. In this case, $(L_{1},L_{2})$ must reside in a
triangular subset $\Omega^{t}\subset\Omega$ (the striped region in Fig. 2).

Figure 2: Schematic of the triangular region $\Omega^{t}$ used in the proof.
$\Omega^{t}$ is illustrated by the diagonal pattern; it is defined as the
region in which $\Delta L_{1}+\Delta L_{2}<0$. The subset of $\Omega^{t}$ in
which both $\Delta L_{1}<0$ and $\Delta L_{2}<0$ is illustrated by the
crosshatched pattern.

On this subset, if both $\Delta L_{1}$ and $\Delta L_{2}$ are negative, then

$\left\lvert\Delta L_{1}-\Delta L_{2}\right\rvert<\left\lvert\Delta
L_{1}+\Delta L_{2}\right\rvert,$ (40)

and moreover $\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert<\left\lvert\Delta L_{1}\right\rvert+\left\lvert\Delta
L_{2}\right\rvert<2L_{\text{ss}}$, so that

$\frac{\pi_{3}}{2}\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert^{2}\left\lvert\Delta L_{1}+\Delta
L_{2}\right\rvert\leq\pi_{3}L_{\text{ss}}\left\lvert\Delta L_{1}+\Delta
L_{2}\right\rvert^{2},$ (41)

and since this term already appears in (39) with a minus sign it follows that
the right-hand side of (39) is negative.

On the other hand, if $\Delta L_{1}$ and $\Delta L_{2}$ have opposite signs,
then we must use a different bound. From the geometry of the relevant
triangular subset $\Omega^{t}\subset\Omega$, we have as before
$\left\lvert\Delta L_{1}+\Delta L_{2}\right\rvert\leq 2L_{\text{ss}}$ so that

$\frac{\pi_{3}}{2}\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert^{2}\left\lvert\Delta L_{1}+\Delta
L_{2}\right\rvert\leq\pi_{3}L_{\text{ss}}\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert^{2}.$ (42)

From the condition that the flagella initially grow if they start at zero
length, we have $\pi_{0}<1$, and from the steady-state equation (20) it
follows that $L_{\text{ss}}\leq 1/2$. Therefore

$\displaystyle 2\pi_{0}J\pi_{3}L_{\text{ss}}\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert^{2}<\pi_{3}J\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert^{2}.$ (43)

It follows that $(\pi_{3}/2)\left\lvert\Delta L_{1}-\Delta
L_{2}\right\rvert^{2}\left\lvert\Delta L_{1}+\Delta L_{2}\right\rvert$ is
bounded above by $\pi_{3}/\pi_{1}\left(\pi_{1}J\left(\Delta L_{1}-\Delta
L_{2}\right)^{2}\right)$, where $\pi_{1}J\left(\Delta L_{1}-\Delta
L_{2}\right)^{2}$ is the right-hand side of (29). Therefore,

$\frac{\mathrm{d}\phi}{\mathrm{d}t}=\frac{1}{2}\left(\frac{\mathrm{d}(\Delta
L_{1}+\Delta
L_{2})^{2}}{\mathrm{d}t}+\frac{\pi_{3}}{\pi_{1}}\frac{\mathrm{d}(\Delta
L_{1}-\Delta L_{2})^{2}}{\mathrm{d}t}\right)<0.$ (44)

∎

Next, we consider whether $\Omega$ is a trapping region closed to the
dynamics. That is, we must consider whether any trajectories originating in
$\Omega$ may eventually leave $\Omega$. First, we show that no trajectories
cross the $L_{1}+L_{2}=1$ boundary. To show that all trajectories move inward
from this boundary, we calculate

$\displaystyle\begin{pmatrix}L_{1}\\\
L_{2}\end{pmatrix}^{\prime}\Big{|}_{L_{1}+L_{2}=1}\cdot\begin{pmatrix}-1\\\
-1\end{pmatrix}$
$\displaystyle=-(J(1-1)-\pi_{0}-\pi_{1}JL_{1})-(J(1-1)-\pi_{0}-\pi_{1}JL_{2})$
$\displaystyle=2\pi_{0}+\pi_{1}J(L_{1}+L_{2})>0.$ (45)

This proves that trajectories flow inward cannot leave $\Omega$ across the
$L_{1}+L_{2}=1$ boundary. See Figure 3 for an illustration of $\Omega$ as well
as the slope field and contours of $\phi(L_{1},L_{2})$.

Figure 3: Dynamics of the active disassembly model. The slope field is plotted
and contours $\phi(L_{1},L_{2})=C$ of the Lyapunov function
$\phi(L_{1},L_{2})$ are shown for several values of $C$. The black asterisk
denotes the unique steady state, the black curves represent boundaries of
$\Omega$, and the red diamonds represent the points $(L_{\text{crit}},0)$ and
$(0,L_{\text{crit}})$ to which trajectories contract from the corner regions.
The parameters used are $\pi_{0}=0.1$, $\pi_{1}=0.15$, $\pi_{2}=1$, and
$\pi_{3}=2$.

In addition, for $\Omega$ to be a trapping region, no trajectories should
escape the first quadrant. In terms of $\Omega$, this means that no
trajectories cross the $L_{1}$ axis on the interval $0<L_{1}\leq 1$, and no
trajectories cross the $L_{2}$ axis on the interval $0<L_{2}\leq 1$. However,
the model eqs. 8 and 9 does not satisfy this requirement. To the contrary,

$\displaystyle\begin{pmatrix}L_{1}\\\
L_{2}\end{pmatrix}^{\prime}\Big{|}_{0<L_{1}\leq
1,L_{2}=0}\cdot\begin{pmatrix}0\\\ 1\end{pmatrix}$
$\displaystyle=J(1-L_{1})-\pi_{0},$ (46)

and upon solving the quadratic equation resulting from $J(T-L_{1})-\pi_{0}=0$
one finds that there is a unique positive root

$L_{\text{crit}}=\frac{1+\pi_{0}\pi_{2}}{2\pi_{0}\pi_{3}}\left(\sqrt{1+\frac{4\left(1-\pi_{0}\right)\pi_{0}\pi_{3}}{\left(1+\pi_{0}\pi_{2}\right)^{2}}}-1\right).$
(47)

Therefore, to prevent trajectories that cross the $L_{1}$-axis resulting in
unphysical negative lengths, we construct a hybrid system by setting
$L_{2}^{\prime}=0$ when $L_{1}>L_{\text{crit}}$ on $L_{2}=0$. This results in
trajectories that flow down to the $L_{1}$-axis, move horizontally according
to this modified dynamics with $L_{2}^{\prime}=0$ until they reach
$(L_{1},L_{2})=(L_{\text{crit}},0)$, at which point they rejoin the full
dynamics and move toward steady state. The dynamics are modified symmetrically
to prevent trajectories from crossing the $L_{2}$-axis as well.

To characterize the trapping region of the full dynamics, we therefore define
the trajectory $\omega_{1}(t)$ going backward in time from the point
$(L_{\text{crit}},0)$, as defined by

$\displaystyle\omega_{1}(0)=(L_{\text{crit}},0),$ (48)

up until it hits the diagonal $L_{1}+L_{2}=1$. We define another trajectory
$\omega_{2}(t)$ similarly as the one going backward in time from the point
$(0,L_{\text{crit}})$. The points from which the trajectories $\omega_{1}(t)$
and $\omega_{2}(t)$ originate are illustrated by red diamonds and the
trajectories are shown as red curves in Fig. 3.

As a consequence of the modified dynamics, trajectories originating in the
corner regions of phase space bounded by $\omega_{1}(t)$ and $\omega_{2}(t)$
consist of three distinct phases: first, there is flow down to the $L_{1}$ or
$L_{2}$-axis, followed next by horizontal or vertical motion to the point
$(L_{\text{crit}},0)$ or $(0,L_{\text{crit}})$, respectively, and finally
leading to recovery by the full dynamics back to steady state. One interesting
consequence of this behavior is its interpretation in the context of severing
experiments. We find that the typical severing protocol, in which the two
flagella first reach steady-state before severing takes place, does not lead
to states within the corner region (purple trajectory in Fig. 4).

Figure 4: Upon severing, the pool size instantaneously decreases by the
steady-state length $L_{\text{ss}}$. This leads to a rescaling of parameters.
(a) A severing experiment is illustrated by plotting two trajectories in the
$L_{1}$-$L_{2}$ phase plane in green and purple, (b) and (c) Length vs. time
corresponding respectively to the green and purple trajectories in (a). To
allow for direct comparison to Fig. 3, the nondimensional parameters are
obtained using the initial tubulin pool size prior to severing. The green
trajectory is an example in which the flagella are initially equal at a value
greater than the steady-state length. There is an initial delay $\tau$ as the
longer flagellum shortens, during which the severed flagellum remains at zero
length (see inset in (c)).

Fig. 4 shows a representative simulation of a severing experiment. The code
used to generate this simulation and the corresponding figures is included in
our public code repository located at
https://github.com/thomasgfai/StableFlagellarLengths.

If the system is initially out of steady-state, however, severing may result
in a state outside the trapping region. This yields two phases of flow after
severing: first, a horizontal flow to the accumulation point
$L_{\text{crit}}$, and second a flow by the full dynamics to reach the global
steady-state (green trajectory in Fig. 4). As shown in the inset of Fig. 4(c),
this first phase of horizontal flow leads to an apparent delay in the response
of the cut flagellum.

One may reasonably ask the question of how a model derived based on a physical
principles can lead to unphysical negative lengths. In this case, the origin
of unphysical trajectories is the basal term $d_{0}$ in the disassembly model,
which assumes there is a nonzero rate of disassembly even when the flagellum
is very small. In the corner regions of the phase space, this basal term is
greater than the rate of assembly and there is a negative growth rate at zero
length because the long second flagellum leads to a tubulin pool that is
nearly depleted. Note that this issue does not arise in agent-based models in
which each tubulin monomer is represented explicitly, since in that case only
assembly can occur when the flagellar length is zero. From this perspective,
the modified dynamics essentially ensure that no more tubulin can be removed
once there are no more tubulin monomers remaining on the flagellum.

###  2.1 Nonexistence of Oscillations: Alternate Proof

Recall that we have an invariant set on the domain $\Omega$ by construction.
In this section, we use an alternative approach to exclude oscillations using
standard methods in dynamical systems. The proof will primarily focus on the
region of $\Omega$ containing the fixed point $(L_{\text{ss}},L_{\text{ss}})$
that is bounded by the lines $L_{2}=L_{1}\pm\varepsilon$:

$\Omega_{\varepsilon}:=\\{(L_{1},L_{2})\,|\,0\leq L_{1},\,0\leq
L_{2},\,L_{1}+L_{2}\leq T,\,L_{2}\leq
L_{1}+\varepsilon,\,L_{1}-\varepsilon\leq L_{2}\\},$

where $0<\varepsilon\leq L_{\text{crit}}$. An example of this domain is shown
in Figure 5 for $\varepsilon=L_{\text{crit}}$ (checkerboard region). The
boundary of $\Omega_{\varepsilon}$ for $\varepsilon=0.2$ is shown using dotted
lines, labeled $L_{2}=L_{1}\pm\varepsilon$. The smaller the value of
$\varepsilon$, the thinner the domain $\Omega_{\varepsilon}$ will be. We will
first show that $\Omega_{\varepsilon}$ is a trapping region (invariant set)
for each $\varepsilon$ arbitrarily small, by showing that the flows across the
sections defined by the lines $L_{2}=L_{1}\pm\varepsilon$ are always inward
towards the fixed point $L_{ss}$. This condition is sufficient because we have
already shown that flow along the relevant boundaries of $\Omega$ point
inwards, specifically, $(L_{1},0)$ and $(0,L_{2})$ for
$L_{1},L_{2}\in[0,L_{\text{crit}}]$, and $L_{1}+L_{2}=1$. To prove the non-
existence of oscillations, we assume it exists and show that a portion of the
oscillation must escape $\Omega_{\varepsilon}$ for some $\varepsilon$ and can
not return, which is a contradiction.

Figure 5: Domains of the alternative proof. $\Omega_{L_{\text{crit}}}$
(checkerboard region) denotes the region $\Omega_{\varepsilon}$ for
$\varepsilon=L_{\text{crit}}$. Dotted lines show an example boundary of
$\Omega_{\varepsilon}$ for $\varepsilon=0.2$. $\Pi_{1}$ (white region) denotes
the region between $\Omega_{L_{\text{crit}}}$ and the trajectories
$\omega_{1}(t)$ and $\omega_{2}(t)$ (red curves). $\Pi_{2}$ (striped region)
denotes the region outside of $\Omega_{L_{\text{crit}}}$ and $\pi_{1}$ but
still inside the overall triangular domain $\Omega$. The parameters are the
same as in Figure 4.

######  Claim.

$\Omega_{\varepsilon}$ is an invariant set.

In the ensuing proof, we show that the flows across the sections defined by
the lines $L_{2}=L_{1}\pm\varepsilon$ are always inward towards the fixed
point $L_{ss}$. Combined with the inward flow along the boundaries of
$\Omega$, this inward flow across the lines $L_{2}=L_{1}\pm\varepsilon$ is
sufficient to prove the claim.

###### Proof.

To show that the vector field points inwards for a given
$\Omega_{\varepsilon}$, we show that $dL_{1}/dt>dL_{2}/dt$ along
$L_{2}=L_{1}+\varepsilon$ and
$\mathrm{d}L_{1}/\mathrm{d}t<\mathrm{d}L_{2}/\mathrm{d}t$ along
$L_{2}=L_{1}-\varepsilon$. First, we plug in $L_{2}=L_{1}+\varepsilon$ into
the right hand side of eqs. 8 and 9 and take the difference
$\mathrm{d}L_{1}/\mathrm{d}t-\mathrm{d}L_{2}/\mathrm{d}t$:

$\displaystyle\frac{\mathrm{d}L_{1}}{\mathrm{d}t}-\frac{\mathrm{d}L_{2}}{\mathrm{d}t}$
$\displaystyle=J(1-L_{1}-L_{2})-\pi_{0}-\pi_{1}JL_{1}-\left[J(1-L_{1}-L_{2})-\pi_{0}-\pi_{1}JL_{2}\right]$
$\displaystyle=-\pi_{1}JL_{1}+\pi_{1}J(L_{1}+\varepsilon)$
$\displaystyle=\pi_{1}J\varepsilon>0.$

Therefore, $\mathrm{d}L_{1}/\mathrm{d}t>\mathrm{d}L_{2}/\mathrm{d}t$
independent of the choice of $\varepsilon$. Inward flow along the lower
boundary follows using the same argument. Thus, for each $\varepsilon$, the
set ${\Omega}_{\varepsilon}$ is invariant.

We now proceed with a proof by contradiction. Suppose that there exists a
periodic solution. Every oscillating solution must surround at least one fixed
point (see index theory, [19]), and the only fixed point of this system exists
within ${\Omega}_{\varepsilon}$. Therefore, the oscillating solution surrounds
the fixed point $(L_{\text{ss}},L_{\text{ss}})$, and there exists some
$\varepsilon^{*}$ such that parts of the periodic solution lie outside of
${\Omega}_{\varepsilon^{*}}$. (If an oscillating solution is fully contained
within ${\Omega}_{\varepsilon}$ for some $\varepsilon$, to obtain
$\varepsilon^{*}$ we may simply reduce $\varepsilon$ until parts of the
periodic solution lie outside of ${\Omega}_{\varepsilon^{*}}$.) It follows
that solutions must exit ${\Omega}_{\varepsilon^{*}}$, but this contradicts
${\Omega}_{\varepsilon^{*}}$ being an invariant set. ∎

###  2.2 Global Asymptotic Stability: Alternate Proof

######  Claim.

The fixed point $(L_{\text{ss}},L_{\text{ss}})$ is globally asymptotically
stable in $\Omega$.

###### Proof.

We will consider three parts of the domain. First,
$\Omega_{L_{\text{crit}}}:=\Omega_{\varepsilon}$ for
$\varepsilon=L_{\text{crit}}$ (Figure 5, checkerboard region), second, the
region between $\Omega_{L_{\text{crit}}}$ and $\omega_{1}$ which we call
$\Pi_{1}$ (Figure 5, white region within $\Omega$), and third, the region
$\Pi_{2}:=\Omega\setminus(\Omega_{L_{\text{crit}}}\cup\Pi_{1})$ (Figure 5,
striped region). The sets $\Omega_{L_{\text{crit}}}$, $\Pi_{1}$, and $\Pi_{2}$
are taken to be closed.

We begin our proof with $\Pi_{1}$ and $\Pi_{2}$ because they are the most
straightforward. In $\Pi_{1}$, all solutions exit through one of the lines
$L_{2}=L_{1}\pm L_{\text{crit}}$. This observation follows because there are
no fixed points in $\Pi_{1}$ and flows cannot cross the solution $\omega_{1}$.
In $\Pi_{2}$, we force all flows to exit through $(L_{\text{crit}},0)$ or
$(0,L_{\text{crit}})$ by definition. Therefore, all flows originating in
$\Pi_{1}\cup\Pi_{2}$ will enter the trapping region
$\Omega_{L_{\text{crit}}}$.

It remains to show asymptotic stability in $\Omega_{L_{\text{crit}}}$. Note
that our system is continuously differentiable on all of
$\Omega_{L_{\text{crit}}}$ and therefore existence and uniqueness of solutions
holds for all time. Before proceeding with the proof of asymptotic stability,
we recall two basic definitions.

######  Definition 1.

The positive semiorbit is defined as the set

$\gamma(\bm{p})=\\{\bm{x}\in\mathbb{R}^{2}\,:\,\bm{x}=\phi(t,\bm{p}),\,\text{for
some}\,\,\,t\in(0,\infty)\\}.$

######  Definition 2.

The $\omega$-limit set of a point $\bm{x}\in M$, $\omega(\bm{x})$, is the set
of those points $\bm{y}\in M$ for which there exists a sequence
$t_{n}\rightarrow\infty$ with $\phi(t_{n},\bm{x})\rightarrow\bm{y}$, where
$\phi$ is the flow of an autonomous ODE $\dot{\bm{x}}=\bm{F}(\bm{x})$, and $M$
is an open subset of $\mathbb{R}^{2}$ [20]. The omega-limit set
$\omega(\gamma)$ for a trajectory $\gamma$ is defined analogously.

In the case that all solutions are bounded, only certain types of
$\omega$-limit sets are possible. This is a consequence of the following
theorem (Theorem 7.4.2 in [21]):

######  Theorem 1.

Let $\gamma$ be a positive semiorbit contained in a compact subset $K$ of
$\mathbb{R}^{2}$ and suppose that $K$ contains only a finite number of
equilibrium points. Then

  1. 1.

$\omega(\gamma)$ consists of a single point $\bm{p}$ which is an equilibrium
point of $\bm{F}$, and for an initial condition $\bm{x}(0)\in\gamma$,
$\bm{x}(t)\rightarrow\bm{p}$ as $t\rightarrow\infty$, or

  2. 2.

$\omega(\gamma)$ is a periodic orbit, or

  3. 3.

$\omega(\gamma)$ consists of a finite set of equilibrium points together with
their connecting orbits. Each such orbit approaches an equilibrium point on
$t\rightarrow\infty$ and as $t\rightarrow-\infty$.

This theorem, which follows from the Poincaré-Bendixson theorem, tells us that
if a solution is trapped in a compact domain for all time, then a single
equilibrium must be either a unique asymptotically stable point or a saddle-
like point with at least one pair of connected unstable and stable manifolds.
(The saddle-like observation follows by the third part of the theorem; a
single equilibrium can also be stable in backwards time.)

Because we have eliminated periodic orbits in $\Omega_{L_{\text{crit}}}$ and
shown that our stable fixed point $(L_{\text{ss}},L_{\text{ss}})$ is locally
asymptotically stable and therefore cannot be saddle-like, we must have an
asymptotically stable equilibrium point in $\Omega_{L_{\text{crit}}}$.
Finally, because all solutions in $\Omega\setminus\Omega_{L_{\text{crit}}}$
eventually reach the trapping region $\Omega_{L_{\text{crit}}}$,
$(L_{\text{ss}},L_{\text{ss}})$ is the globally asymptotically stable
equilibrium point in $\Omega$. ∎

###  2.3 Generalization to organisms with $N>2$ flagella

Here we generalize the Lyapunov function to arbitrary $N$. In this case, the
equations corresponding to eqs. 5, 6 and 7 are:

$\displaystyle\frac{\mathrm{d}L_{i}}{\mathrm{d}t}$ $\displaystyle=\gamma
J\left(T-\sum_{j=1}^{N}L_{j}\right)-d_{0}-d_{1}\frac{JL_{i}}{D},\quad
i=1,\dots,N,$ (49) $\displaystyle J(L_{1},\dots,L_{N})$
$\displaystyle=\frac{k_{\text{on}}M/N}{1+\left(k_{\text{on}}/Nv\right)\sum_{j=1}^{N}L_{j}+\left(k_{\text{on}}/2ND\right)\sum_{j=1}^{N}L_{j}^{2}}.$
(50)

To put the equations in dimensionless form, we follow the approach taken in
the case $N=2$ above by introducing the dimensionless variables
$\widetilde{L}_{i}=L_{i}/T$ and $\widetilde{t}=t\gamma k_{\text{on}}M/N$ as
well as the dimensionless parameters $\pi_{0}=Nd_{0}/\gamma k_{\text{on}}MT$,
$\pi_{1}=d_{1}/\gamma D$, $\pi_{2}=k_{\text{on}}T/Nv$, and
$\pi_{3}=k_{\text{on}}T^{2}/ND$. This leads to the following system, which
generalizes eqs. 8, 9 and 10 above:

$\displaystyle\frac{\mathrm{d}\widetilde{L}_{I}}{\mathrm{d}\widetilde{t}}$
$\displaystyle=\widetilde{J}\left(1-\sum_{j=1}^{N}\widetilde{L}_{j}\right)-\pi_{0}-\pi_{1}\widetilde{J}\widetilde{L}_{i},\quad
i=1,\dots,N,$ (51) $\displaystyle\widetilde{J}$
$\displaystyle=\frac{1}{1+\pi_{2}\sum_{j=1}^{N}\widetilde{L}_{j}+\pi_{3}\sum_{j=1}^{N}\widetilde{L}_{j}^{2}}.$
(52)

As in the case $N=2$, (51) leads to a unique fixed point with
$\widetilde{L}_{i}=\widetilde{L}_{\text{ss}}$ for $i=1,\dots,N$ with

$\widetilde{L}_{\text{ss}}=\frac{1+\pi_{0}\pi_{2}+\pi_{1}/N}{2\pi_{0}\pi_{3}}\left(\sqrt{1+\frac{4\left(1-\pi_{0}\right)\pi_{0}\pi_{3}}{N\left(1+\pi_{0}\pi_{2}+\pi_{1}/N\right)^{2}}}-1\right),$
(53)

which decays as $1/N$ as $N\to 0$.

Note that, in the limit $N\to\infty$ for fixed $T$ and $M$, the dimensionless
basal disassembly parameter $\pi_{0}$ satisfies $\pi_{0}\to\infty$. This is
because the injection rate into any one flagellum, which goes as $1/N$,
becomes smaller and smaller in comparison to the $N$-independent basal
disassembly $d_{0}$, resulting in steady-state lengths of zero beyond a
critical value of $N$. However, in the special case $d_{0}=0$, one obtains
flagella of progressively smaller but nonzero lengths as $N\to 0$.

Next, we prove global asymptotic stability by constructing a Lyapunov
function, following similar logic to that used in the $N=2$ case. In the
following, we once again drop the tildes for convenience. Define $\Delta
L_{i}=L_{i}-L_{\text{ss}}$ as before for $i=1,\dots,N$ and let

$\phi(L_{1},\dots,L_{N})=\frac{1}{2}\left[\left(\sum_{i=1}^{N}\Delta
L_{i}\right)^{2}+\frac{2\pi_{3}}{\pi_{1}}\sum_{i\neq j}\left(\Delta
L_{i}-\Delta L_{j}\right)^{2}\right].$ (54)

######  Claim.

$\phi$ is a Lyapunov function on $\Omega$.

###### Proof.

We first note that setting the derivative to zero in (51) yields the identify

$\displaystyle
0=J_{\text{ss}}\left(1-NL_{\text{ss}}\right)-\pi_{0}-\pi_{1}J_{\text{ss}}L_{\text{ss}}.$
(55)

where

$\displaystyle
J_{\text{ss}}=\frac{1}{1+N\pi_{2}L_{\text{ss}}+N\pi_{3}L_{\text{ss}}^{2}}.$
(56)

Following the logic of the proof in the case $N=2$ results in the equations

$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\sum_{i=1}^{N}\Delta
L_{i}\right)^{2}$
$\displaystyle=N\left(J-J_{\text{ss}}\right)\left(1-NL_{\text{ss}}-\pi_{1}L_{\text{ss}}\right)\sum_{i=1}^{N}\Delta
L_{i}-\left(N+\pi_{1}\right)J\left(\sum_{i=1}^{N}\Delta L_{i}\right)^{2}$ (57)
$\displaystyle\frac{1}{2}\frac{\mathrm{d}}{\mathrm{d}t}\left(\Delta
L_{i}-\Delta L_{j}\right)^{2}$ $\displaystyle=-\pi_{1}J\left(\Delta
L_{i}-\Delta L_{j}\right)^{2},\quad\text{for }i\neq j.$ (58)

From the steady-state identify (55), we have

$1-NL_{\text{ss}}-\pi_{1}L_{\text{ss}}=\frac{\pi_{0}}{J_{\text{ss}}},$ (59)

so that the first term on the right-hand side of (57) may be rewritten as

$\displaystyle
N\left(J-J_{\text{ss}}\right)\left(1-NL_{\text{ss}}-\pi_{1}L_{\text{ss}}\right)\sum_{i=1}^{N}\Delta
L_{i}$
$\displaystyle=\frac{N\pi_{0}}{J_{\text{ss}}}\left(J-J_{\text{ss}}\right)\sum_{i=1}^{N}\Delta
L_{i}$
$\displaystyle=N\pi_{0}J\left(J_{\text{ss}}^{-1}-J^{-1}\right)\sum_{i=1}^{N}\Delta
L_{i}.$ (60)

The expression above may be simplified further by using the definition of the
flux to obtain

$\displaystyle J_{\text{ss}}^{-1}-J^{-1}=-\pi_{2}\sum_{i=1}^{N}\Delta
L_{i}-\pi_{3}\sum_{i=1}^{N}\left(L_{i}^{2}-L_{\text{ss}}^{2}\right),$ (61)

and noting that

$\displaystyle\sum_{i=1}^{N}\left(L_{i}^{2}-L_{\text{ss}}^{2}\right)$
$\displaystyle=\sum_{i=1}^{N}\left(L_{i}+L_{\text{ss}}\right)\Delta L_{i}$
$\displaystyle=\sum_{i=1}^{N}L_{i}\Delta
L_{i}+L_{\text{ss}}\sum_{i=1}^{N}\Delta L_{i}.$ (62)

Moreover, we may rewrite the first term on the right-hand side above by adding
and subtracting $(\frac{1}{N}\sum_{i=1}^{N}L_{i})\sum_{i=1}^{N}\Delta L_{i}$
to obtain

$\displaystyle\sum_{i=1}^{N}L_{i}\Delta L_{i}$
$\displaystyle=\left(\frac{1}{N}\sum_{i=1}^{N}L_{i}\right)\sum_{i=1}^{N}\Delta
L_{i}+\sum_{i=1}^{N}\left(L_{i}-\frac{1}{N}\sum_{i=1}^{N}L_{i}\right)\Delta
L_{i}$
$\displaystyle=\left(\frac{1}{N}\sum_{i=1}^{N}L_{i}\right)\sum_{i=1}^{N}\Delta
L_{i}+\frac{1}{N}\sum_{i=1}^{N}\left(\sum_{j\neq i}\Delta L_{i}-\Delta
L_{j}\right)\Delta L_{i}.$ (63)

Combining eqs. 60, 61, 62 and 63 yields

$\displaystyle
N\pi_{0}J\left(J_{\text{ss}}^{-1}-J^{-1}\right)\sum_{i=1}^{N}\Delta L_{i}=$
$\displaystyle-N\pi_{0}J\left(\pi_{2}\left(\sum_{i=1}^{N}\Delta
L_{i}\right)^{2}+\pi_{3}L_{\text{ss}}\left(\sum_{i=1}^{N}\Delta
L_{i}\right)^{2}\right.$
$\displaystyle\left.+\frac{\pi_{3}}{N}\left(\sum_{i=1}^{N}L_{i}\right)\left(\sum_{i=1}^{N}\Delta
L_{i}\right)^{2}+\frac{\pi_{3}}{N}\left(\sum_{i=1}^{N}\Delta
L_{i}\right)\sum_{i\neq j}\left(\Delta L_{i}-\Delta L_{j}\right)^{2}\right).$
(64)

The first three terms in parentheses on the right-hand side of (64) are
clearly positive. The final term may be bounded according to

$\left|\frac{\pi_{3}}{N}\left(\sum_{i=1}^{N}\Delta L_{i}\right)\sum_{i\neq
j}\left(\Delta L_{i}-\Delta
L_{j}\right)^{2}\right|\leq\frac{\pi_{3}}{N}\left|\sum_{i=1}^{N}\Delta
L_{i}\right|\left|\sum_{i\neq j}\left(\Delta L_{i}-\Delta
L_{j}\right)^{2}\right|\leq\frac{2\pi_{3}}{N}\sum_{i\neq j}\left(\Delta
L_{i}-\Delta L_{j}\right)^{2},$ (65)

where the final inequality follows from $\left|\sum_{i=1}^{N}\Delta
L_{i}\right|\leq\left|\sum_{i=1}^{N}L_{i}\right|+\left|\sum_{i=1}^{N}L_{\text{ss}}\right|\leq
2$, by (55) and the definition of $\Omega$, which asserts that the sum of
lengths must not exceed the total pool.

Therefore, substituting eqs. 64 and 65 into eqs. 57 and 58 results in

$\displaystyle\frac{\mathrm{d}\phi}{\mathrm{d}t}$ $\displaystyle\leq
2\pi_{3}\sum_{i\neq j}\left(\Delta L_{i}-\Delta
L_{j}\right)^{2}-\frac{2\pi_{3}}{\pi_{1}}\sum_{i\neq
j}\left(\pi_{1}J\left(\Delta L_{i}-\Delta L_{j}\right)^{2}\right)$
$\displaystyle\leq 2\pi_{3}\sum_{i\neq j}\left(\Delta L_{i}-\Delta
L_{j}\right)^{2}-2\pi_{3}\sum_{i\neq j}\left(\Delta L_{i}-\Delta
L_{j}\right)^{2}\leq 0,$ (66)

with equality obtained at the unique steady-state, i.e. $\phi$ satisfies the
requirement of a Lyapunov function. ∎

Note that the Lyapunov function we construct for general $N$ does not reduce
to the Lyapunov function constructed previously in the case $N=2$. There is an
additional factor of two in the second term of $\phi$. However, there is no
inconsistency; the Lyapunov function is not unique. Still, the Lyapunov
function constructed previously provides a tighter bound in the case $N=2$,
and we conjecture that the factor of two in the second term of the Lyapunov
function in (54) is not essential.

##  3 Discussion

We have studied the properties of the dynamical system proposed as an active
disassembly model of flagellar length control proposed in [13]. We have
described the region of validity for the model and used the method of Lyapunov
function as well as the Poincare-Bendixson to prove the existence of an
asymptotically stable steady-state. Although the Poincare-Bendixson method is
limited to systems of two equations, the Lyapunov approach is shown to
generalize to $N>2$.

Interestingly, we find that the phase space of physically-allowable lengths is
not a trapping region with respect to the dynamics. Therefore, we must modify
the dynamics in the corner regions from which trajectories are carried outside
the first quadrant. This results in two accumulation points; however, we show
that the unique dynamics associated with these accumulation points are not
likely to arise in application.

The generality of the approach taken here remains an open question. The proofs
given here are independent of the model parameters. However, if the model were
perturbed slightly to involve additional terms, it is not clear whether the
same Lyapunov function approach would still apply. Characterizing the family
of models for which this approach succeeds would be a valuable step toward
establishing the mathematical properties of balance-point models more
generally.

A related open question is whether there exists an optimal Lyapunov function
in the sense of providing the tightest possible bound, and if so, whether this
may provide biological intuition into the energetics of length control.

Finally, we discuss aspects of length control in the case of a large number of
organelles. Indeed, this is a biologically-relevant case given that single-
cell organisms such as _Paramecium_ cells swim using on the order of five
thousand cilia. Depending on which parameters are fixed in the limit
$N\to\infty$, our dimensional analysis of the active disassembly model yields
different limiting behaviors. A remarkable consequence of the model is that it
provides a universal length-control mechanism regardless of organelle number,
and an important next step would be to establish the limits of the model by
testing its predictions on different organisms over a range of $N$.

##  4 Acknowledgments

We thank Ariel Amir and Te Cao for providing feedback on an early version of
this manuscript. We acknowledge useful discussions with Prathitha Kar,
Lishibanya Mohapatra, and Jane Kondev. We acknowledge support under National
Science Foundation grant DMS-1913093 (TGF) and National Institute of Health
grant T32 NS007292 (YP).

## References

  * [1] D. Zink, A. H. Fischer, J. A. Nickerson, Nuclear structure in cancer cells, Nature Reviews Cancer 4 (9) (2004) 677–687.
  * [2] J. B. S. Haldane, On being the right size, Harper’s Magazine 152 (1926) 424–427.
  * [3] W. F. Marshall, Cell geometry: How cells count and measure size, Annual Review of Biophysics 45 (2016) 49–64.
  * [4] R. Milo, R. Phillips, Cell biology by the numbers, Garland Science, 2015.
  * [5] K. A. Wemmer, W. F. Marshall, Flagellar length control in chlamydomonas—a paradigm for organelle size regulation, International review of cytology 260 (2007) 175–212.
  * [6] B. L. Gutman, K. K. Niyogi, Chlamydomonas and arabidopsis. a dynamic duo, Plant Physiology 135 (2) (2004) 607–610.
  * [7] R. E. Goldstein, Green algae as model organisms for biological fluid dynamics, Annual review of fluid mechanics 47 (2015) 343–375.
  * [8] P. A. Salomé, S. S. Merchant, A series of fortunate events: Introducing chlamydomonas as a reference organism, The Plant Cell 31 (8) (2019) 1682–1707.
  * [9] K. Y. Wan, G. Jékely, On the unity and diversity of cilia (2020).
  * [10] R. L. Nguyen, L.-W. Tam, P. A. Lefebvre, The lf1 gene of chlamydomonas reinhardtii encodes a novel protein required for flagellar length control, Genetics 169 (3) (2005) 1415–1424.
  * [11] D. K. Khona, V. G. Rao, M. J. Motiwalla, P. S. Varma, A. R. Kashyap, K. Das, S. M. Shirolikar, L. Borde, J. A. Dharmadhikari, A. K. Dharmadhikari, et al., Anomalies in the motion dynamics of long-flagella mutants of chlamydomonas reinhardtii, Journal of biological physics 39 (1) (2013) 1–14.
  * [12] J. L. Rosenbaum, J. E. Moulder, D. L. Ringo, Flagellar elongation and shortening in chlamydomonas: the use of cycloheximide and colchicine to study the synthesis and assembly of flagellar proteins, The Journal of Cell Biology 41 (2) (1969) 600–619.
  * [13] T. G. Fai, L. Mohapatra, P. Kar, J. Kondev, A. Amir, Length regulation of multiple flagella that self-assemble from a shared pool of components, Elife 8 (2019) e42599.
  * [14] W. F. Marshall, J. L. Rosenbaum, Intraflagellar transport balances continuous turnover of outer doublet microtubules implications for flagellar length control, The Journal of cell biology 155 (3) (2001) 405–414.
  * [15] W. F. Marshall, H. Qin, M. R. Brenni, J. L. Rosenbaum, Flagellar length control system: testing a simple model based on intraflagellar transport and turnover, Molecular biology of the cell 16 (1) (2005) 270–278.
  * [16] K. G. Kozminksi, K. A. Johnson, P. Forscher, J. L. Rosenbaum, A motility in the eukaryotic flagellum unrelated to flagellar beating, Proceedings of the National Academy of Sciences 90 (12) (1993) 5519–5523.
  * [17] A. Chien, S. M. Shih, R. Bower, D. Tritschler, M. E. Porter, A. Yildiz, Dynamics of the ift machinery at the ciliary tip, Elife 6 (2017) e28606.
  * [18] N. L. Hendel, M. Thomson, W. F. Marshall, Diffusion as a ruler: modeling kinesin diffusion as a length sensor for intraflagellar transport, Biophysical journal 114 (3) (2018) 663–674.
  * [19] J. Guckenheimer, P. Holmes, Local bifurcations, in: Nonlinear oscillations, dynamical systems, and bifurcations of vector fields, Springer, 1983, pp. 117–165.
  * [20] G. Teschl, Ordinary differential equations and dynamical systems, Vol. 140, American Mathematical Soc., 2012.
  * [21] N. Lebovitz, Ordinary differential equations, Cengage Learning, 2019.