# Characterization of laser-induced ionization dynamics in solid dielectrics

Peter Jürgens1<EMAIL_ADDRESS>Benjamin Liewehr2 Björn Kruse2 Christian
Peltz2 Tobias Witting1 Anton Husakou1 Arnaud Rouzeé1 Mikhail Ivanov1 Thomas
Fennel2 Marc J. J. Vrakking1 Alexandre Mermillod-Blondin1 1Max-Born-Institute
for Nonlinear Optics and Short Pulse Spectroscopy, Max-Born-Strasse 2A,
D-12489 Berlin, Germany
2Institute of Physics, University of Rostock, Albert-Einstein-Str. 23, D-18059
Rostock, Germany

###### Abstract

The formation of an electron-hole plasma during the interaction of intense
femtosecond laser pulses with transparent solids lies at the heart of
femtosecond laser processing. Advanced micro- and nanomachining applications
require improved control over the excitation characteristics. Here, we relate
the emission of low-order harmonics to the strong laser-field-induced plasma
formation. Together with a measurement of the total plasma density we identify
the contribution of two competing ionization mechanisms - strong-field and
electron-impact ionization.

††preprint: APS/123-QED

##  I Introduction

The precise modification of solid targets with ultrashort laser pulses (USP)
has become a key technology in various fields ranging from data storage Glezer
_et al._ (1995) to waveguide writing Davis  _et al._ (1996); Miura  _et al._
(1997); Mittholiya  _et al._ (2017) and a multitude of cutting purposes Nisar
_et al._ (2013); Ahmed  _et al._ (2008). Every permanent modification of the
optical properties of a target irradiated with USP is initiated by the
formation of an electron-hole plasma von der Linde and Schüler (1996);
Schaffer  _et al._ (2001). Motivated by the tremendous potential of ultrashort
pulse laser micro- and nanomachining Gattass and Mazur (2008); Mottay  _et
al._ (2016); Sugioka and Cheng (2014); Phillips  _et al._ (2015); Malinauskas
_et al._ (2016); Furch  _et al._ (2019), petahertz electronics Garg  _et al._
(2016); Yang  _et al._ (2020); Goulielmakis  _et al._ (2007); Schiffrin  _et
al._ (2013); Schultze  _et al._ (2013); Sederberg  _et al._ (2020) and
coherent control Greenland  _et al._ (2010); Horstmann  _et al._ (2020), huge
theoretical and experimental effort has been devoted to the description of
laser-matter interactions at intensities around the threshold for an
irreversible laser-induced modification Lebugle, M. and Sanner, N. and
Varkentina, N. and Sentis, M. and Utéza, O. (2014); Quéré  _et al._ (2001);
Mao, S. S. and Quéré, F. and Guizard, S. and Mao, X. and Russo, R.E., and
Petite, G. and Martin, P. (2004); Winkler  _et al._ (2006). Although the main
excitation mechanisms relevant for USP-driven laser material processing are
unequivocally established as strong-field ionization (SFI) and electron-impact
ionization (IMP), their relative importance generally remains undetermined.
Experimentally distinguishing the influence of SFI and IMP is rather
challenging since both mechanisms lead to an intrapulse increase of the
electron-hole plasma density Kaiser  _et al._ (2000). Indirect measurements
performed to reveal the relative role of SFI and IMP only allow the
qualitative identification of the prevailing excitation mechanism Jupé  _et
al._ (2009); Gertsvolf  _et al._ (2010).
The population transfer provided by SFI not only depends on the intensity
envelope of the laser pulse but also on its sub-cycle structure Yudin and
Ivanov (2001); Zhokhov and Zheltikov (2018); Gertsvolf  _et al._ (2010). While
femtosecond transient absorption and reflectivity measurements have emerged as
reliable tools for the characterization of the cycle-averaged plasma density
buildup Audebert  _et al._ (1994); Martin  _et al._ (1997); Mero  _et al._
(2003); Garcia-Lechuga  _et al._ (2017); Winkler  _et al._ (2017); Møller  _et
al._ (2020), sub-cycle carrier dynamics in solids have been studied by
attosecond transient absorption spectroscopy Schultze  _et al._ (2013, 2014);
Schlaepfer  _et al._ (2018) and attosecond polarization spectroscopy Sommer
_et al._ (2016). Due to the concomitant manifestation of other electronic
processes such as the dynamical Franz-Keldysh effect Lucchini  _et al._
(2016); Otobe  _et al._ (2016), and IMP Kaiser  _et al._ (2000), these time-
domain methods are not capable of isolating the contribution of SFI to the
ionization yield.
Several theoretical as well as experimental studies support the existence of a
unique optical signature of SFI Brunel (1990); Mitrofanov  _et al._ (2011);
Verhoef  _et al._ (2010); Babushkin  _et al._ (2017); Jürgens  _et al._
(2020); Li  _et al._ (2020); Campi  _et al._ (2021) in the frequency domain,
which arises from the nearly stepwise increase of the plasma density
$\rho_{\rm{SFI}}(t)$ twice per optical cycle of the driving laser field Siders
_et al._ (2001); Serebryannikov  _et al._ (2009). However, the relative role
of the involved physical processes in the associated frequency conversion
process remained elusive.
In a recent publication we unequivocally demonstrated SFI-driven, non-
perturbative wave-mixing in strongly exited fused silica below the damage
threshold and identified the previously disregarded injection current as a
dominant source of low-order harmonic generation Jürgens  _et al._ (2020).
Here, we exploit the optical signature of this injection current to
selectively extract the plasma formation induced by SFI. We demonstrate that
the cycle-averaged conduction band electron density increase can be obtained
from an analysis of the temporal signature of an arbitrary harmonic order. We
further combine the time-resolved measurement of low-order, injection
harmonics with time-domain ptychography Spangenberg  _et al._ (2015); Witting
_et al._ (2016) to access the sub-cycle dynamics of $\rho_{\rm{SFI}}(t)$.

##  II Experimental Setup

The experimental apparatus is described in detail in Ref. Jürgens  _et al._
(2020). In brief, an electron-hole plasma was produced in a
$0.5\text{\,}\mathrm{mm}$-thick a-SiO2 sample (UV grade fused silica, Corning
7980) by focusing a short-wavelength infrared (SWIR,
$\lambda_{\textrm{pump}}=$ $2.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$) laser
pulse along with a weak, time-delayed $790\text{\,}\mathrm{nm}$ (near-
infrared, NIR) probe laser pulse using a cross-polarized, close-to-collinear
pump-probe geometry. Employing cross-polarized pump and probe laser beams
ensures that the weak probe laser pulse does not significantly modify the
ionization yield (see Fig. S1). Both laser pulses were obtained using a
Ti:Sapphire regenerative amplifier ($45\text{\,}\mathrm{fs}$,
$3.3\text{\,}\mathrm{mJ}$ pulse energy, $1\text{\,}\mathrm{kHz}$ repetition
rate, linearly polarized) that pumped an optical parametric amplifier (TOPAS
C, Light Conversion). An off-axis parabolic mirror (focal length of
$50\text{\,}\mathrm{mm}$, gold-coated) focused the SWIR beam down to a spot
size of $\sim$$35\text{\,}\mathrm{\SIUnitSymbolMicro m}$ ($1/e^{2}$ beam
diameter in air) within the sample. These focusing conditions induced a local
intensity in the focal region of up to
$15\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. The NIR probe beam was
directed through a central hole in the parabolic mirror at a small angle
($<\,$$5\text{\,}\mathrm{\SIUnitSymbolDegree}$). The spectral properties of
the light emerging from the sample were measured as a function of the pump-
probe delay $\tau$ with the help of a commercial visible / ultraviolet
(VIS/UV) spectrometer (Avantes AvaSpec-HS1024x58/122TEC). Due to the
noncollinear pump-probe geometry the two-color and the single-color wave-
mixing signals propagated in different directions. The spectrometer was
aligned such, that within its spectral range only the two-color, sum-frequency
mixing signals were detected.
Simultaneously, the transmitted NIR probe laser pulse was analyzed with the
help of a photodiode (Ophir PD-10) and an NIR spectrometer (Avantes AvaSpec-
LS2048), in order to quantify the absorption by the plasma and modifications
of the spectral properties, respectively.

##  III Results

A spectrogram displaying the emitted radiation as a function of the relative
delay $\tau$ between both laser pulses is shown in Fig. 1(a). An intense,
transient, harmonic emission was observed at frequencies

$\omega_{n}=2n\times\omega_{\textrm{pump}}+\omega_{\textrm{probe}}\quad\text{with}\quad
n=1,2,3,4$ (1)

(indicated by the vertical dashed lines).

Figure 1: Two-color wave-mixing in fused silica. (a) Time-frequency map of the
two-color wave-mixing signal obtained by exciting the fused silica sample with
a pump intensity of $13\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. The
vertical dashed lines indicate the positions of the first four harmonic orders
predicted by Eq. 1. (b) Yield of the $n=1$ wave-mixing signal as a function of
the pump laser intensity. Figure 1(a) was measured at the point where the
$n=1$ yield significantly deviates from the quadratic fit (dashed line).

Figure 1(b) displays the yield of the $n=1$ wave-mixing signal as a function
of the pump laser intensity. At intensities
$\leq$$12\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ the nearly quadratic
dependence of the measured yield on the pump laser intensity indicates that
the harmonic formation is dominated by Kerr-type four-wave-mixing, using two
pump photons and one probe photon. The deviation from the parabolic dependence
at pump laser intensities approaching the damage threshold of the SiO2
substrate (i.e. at $\geq$ $14\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$)
results from a change of the prevailing wave-mixing mechanism from a Kerr-type
nonlinearity at low to intermediate intensity to SFI at high intensity.

Following Ref. Jürgens  _et al._ (2019) the time-derivative of the current
density associated with the strong-field-induced population transfer from the
valence to the conduction band and the subsequent field-driven classical
carrier dynamics in the conduction band is given by Geissler  _et al._ (1999)

$\frac{\partial\textbf{j}}{\partial
t}=q_{e}n_{0}\left[\underbrace{\frac{q_{e}}{m_{e}}\textbf{E}\rho_{\rm{SFI}}}_{\rm{Brunel}}+\underbrace{\textbf{v}_{0}\dot{\rho}_{\rm{SFI}}}_{\rm{velocity}}+\underbrace{\frac{\partial}{\partial
t}(\textbf{x}_{0}\dot{\rho}_{\rm{SFI}})}_{\rm{injection}}\right]$ (2)

where $q_{e}=-e$ is the electron charge, $n_{0}$ is the SiO2 molecular
density, $m_{e}$ is the electronic mass, E, the driving laser field,
$\rho_{\rm{SFI}}$ the density of electrons in the conduction band and
$\dot{\rho_{\rm{SFI}}}$ the SFI rate. The parameters $\textbf{x}_{0}$ and
$\textbf{v}_{0}$ denote the spatial displacement and the velocity of the
electrons after injection into the conduction band. Regarding the associated
contributions to nonlinear wave-mixing, it has been demonstrated that the
spatial displacement of electrons during the transition from the valence to
the conduction band [in other words the injection current $j_{\rm{inj}}(t)$]
dominates over the Brunel current and the velocity-current over a wide range
of parameters. In a cross-polarized pump-probe geometry where the polarization
of the harmonic emission follows the polarization of the probe laser pulse and
only a single probe photon is used in the frequency conversion process, the
injection current can be simplified to

$\displaystyle j_{\rm{inj}}(t,\tau)$
$\displaystyle=\frac{n_{0}E_{g}\dot{\rho}_{\rm{SFI}}(t)}{E_{\rm{pump}}(t)^{2}}E_{\rm{probe}}(t-\tau)$
$\displaystyle\equiv\sigma_{\rm{inj}}(t)E_{\rm{probe}}(t-\tau).$ (3)

For convenience we have introduced the effective conductivity

$\displaystyle\sigma_{\rm{inj}}(t)$
$\displaystyle=\frac{n_{0}E_{g}\dot{\rho}_{\rm{SFI}}(t)}{E_{\rm{pump}}(t)^{2}}$
$\displaystyle\propto\frac{\dot{\rho}_{\rm{SFI}}(t)}{E_{\rm{pump}}(t)^{2}}$
$\displaystyle\approx\frac{\dot{\rho}_{\rm{SFI}}(t)}{I_{\rm{pump}}(t)}$ (4)

as the quantity that will be retrieved from a measurement of the injection
harmonics as a function of the pump-probe delay. It should be emphasized that
replacing $E_{\rm{pump}}(t)^{2}$ by the cycle-averaged intensity of the SWIR
pump laser pulse $I_{\rm{pump}}(t)$ generally constitutes a strong
approximation whose justification has to be checked carefully. Here, however,
it is applicable due to the strongly localized character and the high degree
of nonlinearity of $\dot{\rho}_{\rm{SFI}}(t)$ (see supplemental document). As
the effective, strong-field-induced conductivity is nonzero only near the
crests of individual laser half cycles, the field in the denominator can
effectively be replaced by its amplitude without introducing a substantial
error.
In the next step, the SFI rate $\dot{\rho}_{\rm{SFI}}(t)$ can be decomposed
into the sum of a smooth, cycle-averaged component -
$\dot{\rho}_{\rm{SFI}}^{\rm{ca}}(t)$ \- and an oscillating, sub-cycle
component - $\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)$ (for details see
supplemental document). Only the sub-cycle plasma formation due to SFI leads
to the harmonic emission that is reported in Fig. 1(a). Finally, the
experimentally measured harmonic yield $Y_{\rm{inj}}(\omega,\tau)$ can be
related to the sub-cycle ionization dynamics according to (for a derivation
see supplemental document):

$\displaystyle Y_{\rm{inj}}(\omega,\tau)$
$\displaystyle\propto\left|\int\frac{\partial
j_{\rm{inj}}^{\rm{sub}}(t,\tau)}{\partial t}e^{i\omega t}dt\right|^{2}$
$\displaystyle=\omega^{2}\left|\int\sigma_{\rm{inj}}^{\rm{sub}}(t)E_{\rm{probe}}(t-\tau)e^{i\omega
t}dt\right|^{2}$ (5)

###  III.1 Cycle-averaged ionization dynamics

In order to demonstrate the link of the emission of injection harmonics to the
strong-field-induced plasma formation we focus on the $n=1$ harmonic. Figure
2(a) and (b) show the experimental intensity of the $n=1$ harmonic,
$Y_{n=1}(\tau)$, as a function of the pump-probe delay, compared to the
convolution of the pump and probe laser pulses and the expected signature of a
$\chi^{(3)}$ wave-mixing process ($I_{\rm{probe}}(t-\tau)\otimes
I_{\rm{pump}}(t)^{2}$, with linear convolution operator $\otimes$) at
different pump laser intensities. At $I_{\rm{pump}}=$
$2\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, where four-wave-mixing is
assumed to dominate [see Fig. 1(b)], the spectrally integrated $n=1$ harmonic
signal agrees very well with the $\chi^{(3)}$ signature while at
$I_{\rm{pump}}=$ $13\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$ the
emission is temporally more confined, as a consequence of the significantly
higher order nonlinearity of the generation process.

Figure 2: Linking injection harmonics to cycle-averaged carrier dynamics by a
direct deconvolution. (a) Intensity $Y_{n=1}(\tau)$ of the $n=1$ harmonic as a
function of the pump-probe delay obtained at a pump intensity of
$2\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, together with the pump-
probe convolution and the corresponding temporal profile of a
$\chi^{(3)}$-driven wave-mixing process. (b) Same as in (a) for a pump
intensity of $13\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. (c)
Effective, strong-field-induced conductivity $\sigma_{\rm{inj}}(t)$, derived
from (b), together with the intensity envelope of the pump laser pulse. (d)
Cycle-averaged SFI rate obtained by multiplying $\sigma_{\rm{inj}}(t)$ with
$I_{\rm{pump}}(t)$. (e) Plasma density induced by SFI, resulting from a direct
integration of the SFI rate shown in (d), compared to a numerical result for
$\rho_{\rm{SFI}}^{\rm{ca}}(t)$.

At near-threshold intensities, the dependence of the harmonic intensity on the
pump-probe delay is related to a convolution of the effective conductivity
with the probe laser pulse. Analogous to the plasma density and the SFI rate,
the $\sigma_{\rm{inj}}$ can be divided into a cycle-averaged and a sub-cycle
component. Only the cycle-averaged dynamics of the sub-cycle modulated,
effective conductivity can be extracted from the delay dependence of the
harmonic intensity that is given by
$Y_{\rm{inj}}^{n=1}(\tau)\propto\left|\sigma_{\rm{inj}}^{\rm{ca}}(t)\otimes
E_{\rm{probe}}(t-\tau)\right|^{2}$. By means of a numerical deconvolution
using the temporal profile of the NIR probe laser pulse, the cycle-averaged
contribution of $\sigma_{\rm{inj}}(t)$ can be retrieved [see Fig. 2(c)].
Multiplication of this conductivity with the intensity envelope of the pump
laser pulse $I_{\rm{pump}}(t)$ (see Eq. III) reveals the cycle-averaged SFI
rate - $\dot{\rho}_{\rm{SFI}}^{\rm{ca}}(t)$ \- shown in Fig. 2(d). The cycle-
averaged, strong-field-induced plasma density buildup can be computed by a
direct integration of $\dot{\rho}_{\rm{SFI}}^{\rm{ca}}(t)$ and is compared to
a numerical simulation under our experimental conditions in Fig. 2(e) (the
used numerical model is described in detail in Ref. Jürgens  _et al._ (2020)).
The plasma density obtained from the experiments starts to increase at
slightly earlier times with a slightly lower slope than the numerically
determined one. A possible explanation for the deviation between the numerical
prediction and the experimentally determined plasma formation is an
uncertainty of the temporal profile of the probe laser pulse in the
interaction region. The retrieval of $\rho_{\rm{SFI}}^{\rm{ca}}(t)$ from a
direct deconvolution of $Y_{\rm{inj}}^{n=1}(\tau)$ requires the knowledge of
the temporal envelope of $E_{\rm{probe}}(t)$. Due to propagation through the
fused silica sample and the electron-hole plasma the temporal properties of
the probe laser pulse may be modified leading to an additional uncertainty in
the retrieval of the plasma density buildup. Note that the described
characterization of the cycle-averaged plasma dynamics can be performed using
an arbitrary harmonic order, as all injection harmonics exhibit an identical
temporal signature (for details see supplemental document). The analysis of
higher orders is particularly attractive as they are quasi background-free
with respect to Kerr-type contributions.

###  III.2 Sub-cycle ionization dynamics

The previous analysis of the time-resolved injection harmonics connected the
delay-dependent intensity of the $n=1$ harmonic to the cycle-averaged plasma
dynamics induced by SFI. In what follows, a more robust numerical method is
introduced that does not require an a priori knowledge of the probe laser
field and furthermore gives insight into the sub-cycle dynamics of the SFI-
induced plasma formation.
Equation III contains the product of two complex-valued functions whose phase
information is lost as the measured quantity $Y_{\rm{inj}}(\omega,\tau)$ is a
purely real-valued object. Hence, Eq. III resembles a classical phase-
retrieval problem that is well-known in nonlinear optics from many FROG/X-FROG
applications (see e.g. Trebino and Kane (1993); Trebino  _et al._ (1997);
Linde  _et al._ (1998)).
In order to deduce information on the sub-cycle dynamics of the underlying SFI
process we numerically retrieve $\sigma_{\rm{inj}}^{\rm{sub}}(t)$ by a
ptychographic analysis of $Y_{\rm{inj}}(\omega,\tau)/\omega^{2}$.
Multiplication of the retrieved effective conductivity with the intensity
envelope of the pump laser pulse $I_{\rm{pump}}(t)$ reveals the oscillatory
dynamics of the SFI rate [i.e. $\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)$]. To
reconstruct the full SFI rate, we analyze the temporal envelope of
$\rho_{\rm{SFI}}^{\rm{sub}}(t)$. As detailed in the supplemental document, the
envelope of the retrieved sub-cycle component of the SFI rate corresponds to
the cycle-averaged SFI rate $\dot{\rho}_{\rm{SFI}}^{\rm{ca}}(t)$. Adding up
the two components while fulfilling the condition
$\dot{\rho}_{\rm{SFI}}(t)\geq 0\,\forall\,t$ (i.e.
$\dot{\rho}_{\rm{SFI}}(t)=\alpha\dot{\rho}_{\rm{SFI}}^{\rm{ca}}(t)+\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)$
with $\alpha=|\textrm{min}\\{\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)\\}|$) yields
the full SFI rate that can be directly integrated to retrieve the plasma
formation induced by SFI from the ptychographic analysis of time-resolved,
low-order harmonic spectra.

Figure 3: Relating the time-resolved measurement of injection harmonics to
sub-cycle ionization dynamics . (a). Measured time-frequency map of the two-
color wave-mixing signal at an SWIR excitation intensity of
$13\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$. (b). Retrieved
spectrogram by the ptychographic iterative engine after 1000 iterations. (c).
Effective, strong-field-induced conductivity $\sigma_{\rm{inj}}(t)$ as
retrieved by time-domain ptychography. (d). Reconstructed fast and slow
component of the SFI rate - $\dot{\rho}_{\rm{SFI}}(t)$. (e). Resulting total
ionization rate due to SFI obtained from experimental, two-color, harmonic
spectra. (f). Strong-field-induced plasma formation as retrieved from the
time-resolved measurement of the $n=1$ harmonic compared to numerical
simulations under our experimental conditions.

A combined analysis of the slow (femtosecond timescale in the current
experiment) dynamics that are encoded in the delay dependence of the harmonic
intensity and the spectral properties of the emitted harmonics in principle
allows a full reconstruction of the sub-cycle ionization dynamics induced by
the strong SWIR laser field. However, doing so requires a full (amplitude &
phase) retrieval of the harmonics. Since the neighbouring harmonic orders do
not spectrally overlap [see Fig. 1(a)], the relative phase between two
consecutive harmonics remains unknown and we continue our analysis with the
$n=1$ harmonic only.
Results of the ptychographic analysis of the measured $n=1$ harmonic are
displayed in Fig. 3. Figure 3(a) and (b) show the measured and retrieved time-
frequency map containing the $n=1$ harmonic. The ptychographic algorithm is
able to reconstruct the measured wave-mixing signal with high accuracy down to
the noise level of the experiment. The sub-cycle-modulated, effective
conductivity $\sigma_{\rm{inj}}^{\rm{sub}}(t)$ retrieved by the ptychographic
algorithm is shown in Fig. 3(c). The resulting sub-cycle component of the SFI
rate - $\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)$ is depicted together with the
cycle-averaged, low-frequency contribution in Fig. 3(d) (Note that due to the
decomposition
$\dot{\rho}_{\rm{SFI}}(t)=\dot{\rho}_{\rm{SFI}}^{\rm{ca}}(t)+\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)$
the sub-cycle component of the SFI rate can have negative values, while the
cycle-averaged and the total SFI rate are always $\geq 0$). Combining the fast
and the slow component leads to the total SFI rate depicted in Fig. 3(e).
Direct integration of $\dot{\rho}_{\rm{SFI}}(t)$ yields the total plasma
formation due to SFI reconstructed from the time-resolved measurement of the
$n=1$ harmonic. The resulting $\rho_{\rm{SFI}}(t)$ is shown together with
results obtained by numerical simulations under identical conditions in Fig.
3(f). A nearly stepwise increase of the resulting plasma density with two
steps per optical cycle of the SWIR driving field is found. As the carrier-
envelope phase (CEP) of the pump laser pulses is not actively stabilized in
our experiments and the ptychographic algorithm is not able to retrieve the
absolute phase of $\sigma_{\rm{inj}}^{\rm{sub}}(t)$ our results exhibit a
temporal uncertainty of one optical cycle [see shaded area in Fig. 3(f)].
When compared to the numerically obtained plasma density buildup, a slightly
later liftoff of the $\rho_{\rm{SFI}}(t)$ obtained from our measurement is
observed. The overall slope of the reconstructed $\rho_{\rm{SFI}}(t)$ is found
to be marginally steeper than the numerical prediction. Discrepancies between
the reconstructed and numerically predicted $\rho_{\rm{SFI}}(t)$ might be due
to dynamic modifications of the bandgap Winkler  _et al._ (2017); Lucchini
_et al._ (2016) leading to transiently enhanced or diminished ionization
probabilities.

##  IV Discussion and Outlook

While the cycle-averaged dynamics can be extracted from the temporal signature
of an arbitrary harmonic order (here, we chose the $n=1$ harmonic), the
stepwise increase of the plasma density can be reconstructed by a combined
analysis of the injection harmonics in the time and frequency domain. The
$n=1$ harmonic suffices to characterize the cycle-averaged plasma dynamics and
additionally provides valuable insights into the sub-cycle dynamics. As a
consequence, it is possible to apply the ptychographic reconstruction of
strong-field-induced carrier dynamics on a wide range of materials including
semiconductors with bandgaps in the visible spectral range where only
$\hbar\omega_{n=1}\leq E_{g}$.
In addition to the retrieval of $\sigma_{\rm{inj}}^{\rm{sub}}(t)$ time-domain
ptychography retrieves the temporal profile of the NIR probe laser pulse in
the interaction region. The retrieved duration
($\sim$$65\text{\,}\mathrm{fs}$) differs significantly from the initial
duration ($45\text{\,}\mathrm{fs}$) indicating that the propagation through
the fused silica sample and the electron-hole plasma strongly affect the
temporal structure of the probe pulse. Indeed, using the probe laser pulse
retrieved by the ptychographic iterative engine for the deconvolution
considerably improves the agreement between the retrieved and the simulated
$\rho_{\rm{SFI}}^{\rm{ca}}(t)$ [see Fig. 2(e)].
Finally, we are able to disentangle the plasma dynamics induced by SFI from
the competing ionization and relaxation mechanisms. Together with the
simultaneously measured transmission of the NIR probe laser that carries
information on the total plasma density (see Ref. Jürgens  _et al._ (2019))
this enables the determination of the relative contribution of SFI and IMP.
Figure 4 shows a comparison between the total plasma density obtained from a
time-resolved transmission measurement and the SFI-induced plasma formation
extracted from the time-resolved detection of the $n=1$ injection harmonic.
Here, both densities are normalized to unity as the absolute number of excited
carriers cannot be determined by both presented methods. When comparing the
temporal profile of $\rho_{\rm{tot}}(t)$ and $\rho_{\rm{SFI}}(t)$ it becomes
evident that $\rho_{\rm{SFI}}(t)$ starts to increase at much earlier times (at
$\tau\approx$ $-50\text{\,}\mathrm{fs}$). $\rho_{\rm{tot}}(t)$ lifts off at
$\tau\approx$ $25\text{\,}\mathrm{fs}$ where $\rho_{\rm{SFI}}(t)$ already
reaches $90\%$ of its maximum value. The fact that the experimentally
determined $\rho_{\rm{SFI}}(t)$ reaches its maximum value where no significant
probe absorption occurs indicates that IMP is the dominating excitation
mechanism under our experimental parameters. This is further supported by the
delay between the peak of the pulse and the maximum of $\rho_{\rm{tot}}(t)$ at
$\approx$ $110\text{\,}\mathrm{fs}$ and is in very good agreement with
numerical simulations that predict a complete dominance of IMP for our
experimental configuration. At pump-probe delays $\geq$
$110\text{\,}\mathrm{fs}$, $\rho_{\rm{tot}}(t)$ exhibits an exponential decay
due to ultrafast relaxation processes that are not accounted for in the
determination of $\rho_{\rm{SFI}}(t)$. Therefore, $\rho_{\rm{SFI}}(t)$ stays
at the maximum level for longer time delays.

Figure 4: Comparison of the SFI-induced plasma dynamics extracted from the
ptychographic analysis of the $n=1$ harmonic and the total plasma dynamics
obtained from a time-resolved absorption measurement. The dashed lines show
numerical results for $\rho_{\rm{SFI}}(t)$ and $\rho_{\rm{tot}}(t)$ at a pump
intensity of $17\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$.

In conclusion, we have presented time-resolved experiments on two-color, low-
order harmonic generation in fused silica with the aim of isolating the plasma
density buildup induced by SFI. Our analysis demonstrates a possible
application of the recently identified injection harmonics for ultrafast
metrology in solid state systems. The fact that the detected signal is not
screened by concomitant ionization mechanisms such as IMP, allows a
background-free investigation of strong-field phenomena in bulk solids.
We demonstrated the possibility to extract the cycle-averaged plasma density
buildup due to SFI from a simple numerical analysis of the intensity of an
arbitrary harmonic order. By applying an iterative phase retrieval algorithm
on the time-resolved measurement of the $n=1$ harmonic we revealed the
stepwise character of the field-induced plasma formation. The exact shape of
these steps can be obtained by including higher harmonic orders with a well-
determined relative phase while the exact temporal position of the ionization
bursts is determined by the relative phase between the driving laser field and
the emitted harmonics. Comparing the SFI-induced plasma formation to the total
evolution of the electron-hole plasma provides intimate insights on the
relative contribution of field-driven and impact-driven ionization channels to
the total plasma density buildup. Importantly, SFI-induced signals manifest
well-below the threshold for a permanent modification and thus hold the
potential for an in-situ detection of the onset of optical damage before a
catastrophic damage occurs.

### Funding

Deutsche Forschungsgemeinschaft (IDs DFG ME4427/1-1 and DFG ME4427/1-2,
Heisenberg Grant ID 398382624, SPP1840: QUTIF ID 281272685), Leibniz
Association grant SAW‐LAPTON (K266/2019), Ministry of Education, Science and
Culture of Mecklenburg-Vorpommern, Germany via the project ”NEISS”.

### Acknowledgements

We thank M. Jupé (Laser Zentrum Hannover) for measuring the bandgap of the
samples used in this study.

### Supplemental document

See supplemental document for supporting content.

## References

  * Glezer  _et al._ (1995) E. N. Glezer, M. Milosavljevic, L. Huang, R. J. Finlay, T.-H. Her, J. P. Callan, and E. Mazur, Optics Letters 21, 2023 (1995).
  * Davis  _et al._ (1996) K. M. Davis, K. Miura, N. Sugimoto, and K. Hirao, Optics Letters 21, 1729 (1996).
  * Miura  _et al._ (1997) K. Miura, J. Qiu, H. Inouye, T. Mitsuyu, and K. Hirao, Applied Physics Letters 71, 3329 (1997).
  * Mittholiya  _et al._ (2017) K. Mittholiya, P. K. Anshad, A. K. Mallik, S. Bhardwaj, A. Hedge, A. Bhatnagar, R. Bernard, J. A. Dharmadhikari, D. Mathur, and A. K. Dharmadhikari, Journal of Optics 46, 304 (2017).
  * Nisar  _et al._ (2013) S. Nisar, L. Li, and A. Sheikh, Journal of Laser Applications 25 (2013).
  * Ahmed  _et al._ (2008) F. Ahmed, M. S. Lee, H. Sekita, T. Sumiyoshi, and M. Kamata, Applied Physics A 93, 189 (2008).
  * von der Linde and Schüler (1996) D. von der Linde and H. Schüler, J. Opt. Soc. Am. B 13 (1996).
  * Schaffer  _et al._ (2001) C. B. Schaffer, A. Brodeur, and E. Mazur, Measurement Science and Technology 12 (2001).
  * Gattass and Mazur (2008) R. R. Gattass and E. Mazur, Nature Photonics 2, 219 (2008).
  * Mottay  _et al._ (2016) E. Mottay, X. Liu, H. Zhang, E. Mazur, R. Sanatinia, and W. Pfleging, MRS Bulletin 41, 984 (2016).
  * Sugioka and Cheng (2014) K. Sugioka and Y. Cheng, Light Sci Appl 3 (2014).
  * Phillips  _et al._ (2015) K. C. Phillips, J. H. Gandhi, E. Mazur, and S. K. Sundaram, Advances in Optics and Photonics 7, 684 (2015).
  * Malinauskas  _et al._ (2016) M. Malinauskas, A. Žukauskas, S. Hasegawa, Y. Hayasaki, V. Mizeikis, R. Buividas, and S. Juodkazis, Light: Science & Applications 5 (2016).
  * Furch  _et al._ (2019) F. J. Furch, D. W. Engel, T. Witting, A. Perez-Leija, M. J. J. Vrakking, and A. Mermillod-Blondin, Optics Letters 44, 4267 (2019).
  * Garg  _et al._ (2016) M. Garg, M. Zhan, T. Luu, H. Lakhotia, T. Klostermann, A. Guggenmos, and E. Goulielmakis, Nature 538, 359 (2016).
  * Yang  _et al._ (2020) Y. Yang, M. Tuchetti, P. Vasireddy, W. P. Putnam, O. Karnbach, A. Nardi, F. X. Kärtner, K. K. Berggren, and P. D. Keathley, Nature Communications 11 (2020).
  * Goulielmakis  _et al._ (2007) E. Goulielmakis, V. S. Yakovlev, A. L. Cavalieri, M. Uiberacker, V. Pervak, A. Apolonski, R. Kienberger, U. Kleineberg, and F. Krausz, Science 317, 769 (2007).
  * Schiffrin  _et al._ (2013) A. Schiffrin, T. Paasch-Colberg, N. Karpowicz, V. Apalkov, D. Gerster, S. Mühlbrandt, M. Korbman, J. Reichert, M. Schultze, S. Holzner, J. V. Barth, R. Kienberger, R. Ernstorfer, V. S. Yakovlev, M. I. Stockman, and F. Krausz, Nature 493, 70 (2013).
  * Schultze  _et al._ (2013) M. Schultze, E. M. Bothschafter, A. Sommer, S. Holzner, W. Schweinberger, M. Fiess, M. Hofstetter, R. Kienberger, V. Apalkov, V. S. Yakovlev, M. I. Stockman, and F. Krausz, Nature 493, 75 (2013).
  * Sederberg  _et al._ (2020) S. Sederberg, D. Zimin, S. Keiber, F. Siegrist, M. S. Wismer, V. S. Yakovlev, I. Floss, C. Lemell, J. Burgdörfer, M. Schultze, F. Krausz, and N. Karpowicz, Nature Communications 11 (2020).
  * Greenland  _et al._ (2010) P. T. Greenland, S. A. Lynch, A. F. G. van der Merr, B. N. Murdin, C. R. Pidgeon, B. Redlich, N. Q. Vinh, and G. Aeppli, Nature 465, 1057 (2010).
  * Horstmann  _et al._ (2020) J. G. Horstmann, H. Böckmann, B. Wit, F. Kurtz, G. Storeck, and C. Ropers, Nature 583, 232 (2020).
  * Lebugle, M. and Sanner, N. and Varkentina, N. and Sentis, M. and Utéza, O. (2014) Lebugle, M. and Sanner, N. and Varkentina, N. and Sentis, M. and Utéza, O., J. Appl. Phys. 116 (2014).
  * Quéré  _et al._ (2001) F. Quéré, S. Guizard, and P. Martin, Europhys. Lett. 56, 138 (2001).
  * Mao, S. S. and Quéré, F. and Guizard, S. and Mao, X. and Russo, R.E., and Petite, G. and Martin, P. (2004) Mao, S. S. and Quéré, F. and Guizard, S. and Mao, X. and Russo, R.E., and Petite, G. and Martin, P., Applied Physics A 79, 1695 (2004).
  * Winkler  _et al._ (2006) S. Winkler, I. Burakov, R. Stoian, N. Bulgakova, A. Husakou, A. Mermillod-Blondin, A. Rosenfeld, D. Ashkenasi, and I. Hertel, Applied Physics A 84, 413 (2006).
  * Kaiser  _et al._ (2000) A. Kaiser, B. Rethfeld, M. Vicanek, and G. Simon, Physical Review B 61, 11437 (2000).
  * Jupé  _et al._ (2009) M. Jupé, L. Jensen, A. Melnikaitis, V. Sirutkaitis, and D. Ristau, Optics Express 17 (2009).
  * Gertsvolf  _et al._ (2010) M. Gertsvolf, M. Spanner, D. M. Rayner, and P. B. Corkum, J. Phys. B: At. Mol. Opt. Phys. 43, 1 (2010).
  * Yudin and Ivanov (2001) G. L. Yudin and M. Y. Ivanov, Physical Review A 64 (2001).
  * Zhokhov and Zheltikov (2018) P. A. Zhokhov and A. M. Zheltikov, Scientific Reports 8 (2018).
  * Audebert  _et al._ (1994) P. Audebert, P. Daguzan, A. Dos Santos, J. C. Gauthier, J. Geindre, S. Guizard, G. Hamoniaux, K. Krastev, P. Martin, G. Petite, and A. Antonetti, Physical Review Letters 73, 1990 (1994).
  * Martin  _et al._ (1997) P. Martin, S. Guizard, P. Daguzan, G. Petite, P. D’Oliveira, P. Meaynadier, and M. Perdrix, Physical Review B 55 (1997).
  * Mero  _et al._ (2003) M. Mero, J. Liu, A. Sabbah, J. Jasapara, K. Starke, D. Ristau, J. K. McIver, and W. Rudolph, in  _Proc. SPIE 4932_ (2003) pp. 202 – 212.
  * Garcia-Lechuga  _et al._ (2017) M. Garcia-Lechuga, L. Haahr-Lillevang, J. Siegel, P. Balling, S. Guizard, and J. Solis, Physical Review B 95, 214114 (2017).
  * Winkler  _et al._ (2017) T. Winkler, L. Haahr-Lillevang, C. Sarpe, B. Zielinski, N. Götte, A. Senftleben, P. Balling, and T. Baumert, Nature Physics 14 (2017).
  * Møller  _et al._ (2020) S. H. Møller, S. T. Andersen, and P. Balling, Phys. Rev. Research 2, 043010 (2020).
  * Schultze  _et al._ (2014) M. Schultze, K. Ramasesha, C. D. Pemmaraju, S. A. Sato, D. Whitmore, A. Gandman, J. S. Prell, L. J. Borja, D. Prendergast, K. Yabana, D. M. Neumark, and S. R. Leone, Science 346, 1348 (2014).
  * Schlaepfer  _et al._ (2018) F. Schlaepfer, M. Lucchini, S. A. Sato, M. Volkov, L. Kasmi, N. Hartmann, A. Rubio, L. Gallmann, and U. Keller, Nature Physics 14, 560 (2018).
  * Sommer  _et al._ (2016) A. Sommer, E. M. Bothschafter, S. A. Sato, C. Jakubeit, T. Latka, O. Razskazovskaya, H. Fattahi, M. Jobst, W. Schweinberger, V. Shirvanyan, V. Yakovlev, R. Kienberger, K. Yabana, N. Karpowicz, M. Schultze, and F. Krausz, Nature 534, 86 (2016).
  * Lucchini  _et al._ (2016) M. Lucchini, S. A. Sato, A. Ludwig, J. Herrmann, M. Volkov, L. Kasmi, Y. Shinohara, K. Yabana, L. Gallmann, and U. Keller, Science 353, 916 (2016).
  * Otobe  _et al._ (2016) T. Otobe, Y. Shinohara, S. Sato, and K. Yabana, Physical Review B 93, 045124 (2016).
  * Brunel (1990) F. Brunel, J. Opt. Soc. Am. B 7, 4 (1990).
  * Mitrofanov  _et al._ (2011) A. V. Mitrofanov, A. J. Verhoef, E. E. Serebryannikov, J. Lumeau, L. Glebov, A. M. Zheltikov, and A. Baltuška, Phys. Rev. Lett. 106, 147401 (2011).
  * Verhoef  _et al._ (2010) A. J. Verhoef, A. V. Mitrofanov, E. E. Serebryannikov, D. V. Kartashov, A. M. Zheltikov, and A. Baltuška, Phys. Rev. Lett. 104, 163904 (2010).
  * Babushkin  _et al._ (2017) I. Babushkin, C. Brée, C. M. Dietrich, A. Demircan, U. Morgner, and A. Husakou, J. Mod. Opt. 64, 10 (2017).
  * Jürgens  _et al._ (2020) P. Jürgens, B. Liewehr, B. Kruse, C. Peltz, D. Engel, A. Husakou, T. Witting, M. Ivanov, M. Vrakking, T. Fennel, and A. Mermillod-Blondin, Nat. Phys. 16, 1035 (2020).
  * Li  _et al._ (2020) L. Li, T. Huang, P. Lan, Y. Zhang, J. Li, X. Zhu, L. He, W. Cao, and P. Lu, “Second-order nonlinear polarization in the strong-field regime,” (2020), arXiv:2010.13298 [physics.optics] .
  * Campi  _et al._ (2021) F. Campi, S. D. R. Abbing, Z.-Y. Zhang, M. van der Geest, and P. M. Kraus, in  _International Conference on X-Ray Lasers 2020_, Vol. 11886, edited by D. Bleiner, International Society for Optics and Photonics (SPIE, 2021) pp. 117 – 124.
  * Siders  _et al._ (2001) C. W. Siders, G. Rodriguez, J. L. W. Siders, F. G. Omenetto, and A. J. Taylor, Phys. Rev. Lett. 87, 263002 (2001).
  * Serebryannikov  _et al._ (2009) E. E. Serebryannikov, A. J. Verhoef, A. Mitrofanov, A. Baltuška, and A. M. Zheltikov, Phys. Rev. A 80, 053809 (2009).
  * Spangenberg  _et al._ (2015) D. Spangenberg, P. Neethling, E. Rohwer, M. H. Brügmann, and T. Feurer, Phys. Rev. A 91, 021803 (2015).
  * Witting  _et al._ (2016) T. Witting, D. Greening, D. Walke, Matia-Hernando, B. T., J. P. Marangos, and J. W. G. Tisch, Optics Letters 41, 4218 (2016).
  * Jürgens  _et al._ (2019) P. Jürgens, M. Vrakking, A. Husakou, R. Stoian, and A. Mermillod-Blondin, Applied Physics Letters 115 (2019).
  * Geissler  _et al._ (1999) M. Geissler, G. Tempea, A. Scrinzi, M. Schnürer, F. Krausz, and T. Brabec, Phys. Rev. Lett. 83, 2930 (1999).
  * Trebino and Kane (1993) R. Trebino and D. J. Kane, Journal of the Optical Society of America A 10, 1101 (1993).
  * Trebino  _et al._ (1997) R. Trebino, K. W. DeLong, D. N. Fittinghoff, J. N. Sweetser, M. A. Krumbügel, B. A. Richman, and D. J. Kane, Review of Scientific Instruments 68, 3277 (1997).
  * Linde  _et al._ (1998) S. Linde, H. Giesse, and J. Kuhl, phys. stat. sol. (b) 206, 119 (1998).
  * Mermillod-Blondin  _et al._ (2019) A. Mermillod-Blondin, P. Jürgens, B. Liewehr, B. Kruse, C. Peltz, T. Witting, A. Husakou, M. Ivanov, T. Fennel, and M. J. J. Vrakking, in  _Laser Applications in Microelectronic and Optoelectronic Manufacturing (LAMOM) XXIV_ , Vol. 10905 (International Society for Optics and Photonics, 2019) p. 109050O.
  * Ammosov  _et al._ (1986) M. V. Ammosov, N. B. Delone, and V. P. Krainov, Sov. Phys. JETP 64, 1191 (1986).

## Supplemental Document

###  IV.1 Two-color plasma formation using parallel and perpendicular
polarizations

Numerical simulations of the plasma formation due to SFI using the model
described in Jürgens  _et al._ (2020) illustrate the importance of the
perpendicular polarizations of the pump and probe laser fields. Figure S1
shows the number of electrons promoted into the conduction band by the pump
laser field only, the pump laser field in combination with a weak, cross-
polarized probe laser field, the pump laser field and a weak parallel-
polarized probe laser field ($I_{\rm{pump}}+I_{\rm{probe}}$ and the pump field
and a probe laser field of thrice the initial intensity
($I_{\rm{pump}}+I_{\rm{probe}}$). While the ionization yield is essentially
the same in the first two cases, in a parallel polarization configuration even
a weak probe laser field is sufficient to substantially modify the total
number of conduction band electrons generated during pump-probe overlap.
Slight changes of the probe laser intensity considerably enhance the cross-
excitation efficiency leading to a clear departure of the resulting plasma
dynamics from the pump-field-induced plasma formation.

Figure S1: Calculated strong-field-induced plasma dynamics due to two-color
excitation of fused silica in different polarization configurations normalized
to the yield obtained from irradiation with the pump laser pulse only.
Simulation parameters: $I_{\rm{pump}}=$
$12\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$, $\lambda_{\rm{pump}}=$
$2100\text{\,}\mathrm{nm}$, $\tau_{\rm{pump}}=$ $140\text{\,}\mathrm{fs}$,
$I_{\rm{probe}}=$ $0.015\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$,
$\lambda_{\rm{probe}}=$ $800\text{\,}\mathrm{nm}$, $\tau_{\rm{probe}}=$
$45\text{\,}\mathrm{fs}$.

###  IV.2 Derivation of $Y_{\rm{inj}}(\omega,\tau)$

Only the sub-cycle dynamics of the SFI rate contribute to the formation of the
harmonics. The sub-cycle component of the injection current relevant for the
emission of the $n=1$ wave-mixing signal, assuming cross-polarized pump and
probe laser fields and a weak probe pulse can be expressed as

$j_{\rm{inj}}^{\rm{sub}}(t,\tau)=\frac{n_{0}E_{g}\dot{\rho}_{\rm{SFI}}^{\rm{sub}}(t)}{E_{\rm{pump}}(t)^{2}}E_{\rm{probe}}(t-\tau)$
(6)

The expression for the measured quantity $Y_{\rm{inj}}(\omega,\tau)$ can be
derived by

$\displaystyle Y_{\rm{inj}}(\omega,\tau)$
$\displaystyle\propto\left|\int\frac{\partial
j_{\rm{inj}}^{\rm{sub}}(t,\tau)}{\partial t}e^{i\omega t}dt\right|^{2}$
$\displaystyle=\omega^{2}\left|\int j_{\rm{inj}}^{\rm{sub}}(t,\tau)e^{i\omega
t}dt\right|^{2}$
$\displaystyle=\omega^{2}\left|\int\frac{n_{0}E_{g}\dot{\rho}^{\rm{sub}}_{\rm{SFI}}(t)}{E_{\rm{pump}}(t)^{2}}E_{\rm{probe}}(t-\tau)e^{i\omega
t}dt\right|^{2}$
$\displaystyle=\omega^{2}\left|\int\sigma_{\rm{inj}}^{\rm{sub}}(t)E_{\rm{probe}}(t-\tau)e^{i\omega
t}dt\right|^{2}$ (7)

All variables are defined in the main text.

###  IV.3 Simplifications in order to analyze the retrieved quantity

The calculation of $\dot{\rho}_{\rm{SFI}}(t)$ from Eq. III requires the exact
knowledge of the pump laser field together with the relative phase between
$E_{\rm{pump}}(t)$ and $\sigma_{\rm{inj}}(t)$, which is not accessible in the
presented experiment. However, we can exploit the fact that the pump laser
field varies slowly compared to the SFI rate due to its high nonlinearity.
Since SFI takes place mainly close to the extrema of $E_{\rm{pump}}(t)$ it is
possible to replace $E_{\rm{pump}}(t)^{2}$ by the intensity envelope
$I_{\rm{pump}}(t)$ of the SWIR pump laser pulse. This approximation is however
only valid for certain experimental parameters. In our case, the high
nonlinearity of the excitation process due to the small photon energy of the
SWIR pump laser pulse ($0.59\text{\,}\mathrm{eV}$) in combination with the
wide bandgap of the fused silica samples ($7.7\text{\,}\mathrm{eV}$) justifies
the simplification. For shorter wavelengths (larger photon energies) or
materials with a smaller bandgap this approximation tends to fail as the
excitation process is temporarily less confined due to the lower nonlinearity
of the SFI process.

Figure S2: Numerical analysis of the assumptions made for
$\sigma_{\rm{inj}}(t)$ in the main text. (c) Comparison of the SFI rate and
$E(t)^{2}$. (b). Normalized quantity
$\dot{\rho}_{\rm{SFI}}(t)/E_{\rm{pump}}(t)^{2}$ compared to the discussed
approximation $\dot{\rho}_{\rm{SFI}}(t)/I_{\rm{pump}}(t)$ and the ionization
rate. (c) Error $\Delta$ resulting from the discussed simplification as a
function of the wavelength.

Figure S2 illustrates this by comparing the SFI rate to the instantaneous
electric field [see Fig. S2(a)]. Due to the high nonlinearity of the SFI
process, ionization takes place only in the vicinity of the peaks of
$E(t)^{2}$ where $E(t)^{2}\approx I(t)$. Figure S2(b) compares the quantities
$\sigma_{\rm{inj}}(t)=\dot{\rho}_{\rm{SFI}}(t)/E_{\rm{pump}}(t)^{2}$ to the
simplification $\dot{\rho}_{\rm{SFI}}(t)/I_{\rm{pump}}(t)$ for the
experimental parameters
($\lambda_{\rm{pump}}=\,$$2.1\text{\,}\mathrm{\SIUnitSymbolMicro m}$,
$I_{\rm{pump}}=\,$ $13\text{\,}\mathrm{TW}\text{\,}{\mathrm{cm}}^{-2}$). No
difference between the exact conductivity and the approximation is observed.
In order to analyze the influence of the wavelength on the validity of the
approximation we can assume a pure multi-photon ionization scenario and
approximate the ionization rate by $\dot{\rho}_{\rm{SFI}}(t)\propto
E_{\rm{pump}}(t)^{2N}$ where $N$ denotes the multi-photon order defined by
$N=]E_{\rm{g}}/\hbar\omega+1[$. Here $E_{\rm{g}}$ is the bandgap of the
material and $]x[$ denotes the integer part of $x$. The mean absolute error
$\Delta$ shown in Fig. S2(b) is calculated as follows:

$\Delta=\frac{1}{M}\sum_{i=0}^{M}\left|\frac{E_{\rm{pump}}^{2N}(t_{i})}{E_{\rm{pump}}^{2}(t_{i})}-\frac{E_{\rm{pump}}^{2N}(t_{i})}{I_{\rm{pump}}(t_{i})}\right|.$
(8)

The observed dependence of $\Delta$ on the excitation wavelength indicates
that the replacement of $E_{\rm{pump}}(t)^{2}$ by the intensity envelope
$I_{\rm{pump}}(t)$ is only valid for long wavelengths (high multi-photon
order) while for shorter wavelengths the approximation is expected to
substantially modify the result.

###  IV.4 Decomposition of the SFI rate

Assuming an intensity profile that is given by the product of an envelope and
an oscillating part $I(t)=A(t)\cos(\omega_{0}t)^{2}$, with amplitude $A(t)$
and central angular frequency $\omega_{0}$, the strong-field ionization rate
$\dot{\rho}_{\rm{SFI}}(t)$ can be decomposed as a Fourier cosine series with a
fundamental frequency of $2\omega_{0}$ (see also Mermillod-Blondin  _et al._
(2019)):

Figure S3: Decomposition of the strong-field-induced plasma density, the SFI
rate and the effective conductivity into a cycle-averaged and a sub-cycle
component. (a)-(c) Strong-field-induced plasma formation -
$\rho_{\rm{SFI}}(t)$. (d)-(f) SFI rate - $\dot{\rho}_{\rm{SFI}}(t)$. (g)-(i)
Strong-field-induced conductivity - $\sigma_{\rm{inj}}(t)$.

$\displaystyle\frac{\partial\rho_{\rm{SFI}}(t)}{\partial t}$
$\displaystyle=B(t)\left[\frac{1}{2}a_{0}+\sum_{n=1}^{\infty}a_{n}(t)\cos(2n\omega_{0}t)\right]$
$\displaystyle=\underbrace{\frac{a_{0}B(t)}{2}}_{\rm{cycle-
averaged}}+\underbrace{B(t)\sum_{n=1}^{\infty}a_{n}(t)\cos(2n\omega_{0}t)}_{\rm{sub-
cycle}}$ (9)

with $B(t)$ being a nonlinear function in $A(t)$ taking into account the
intensity dependence of the SFI rate. Here $n$ denotes the order of the wave-
mixing process and $a_{i}$ the $i$-th Fourier coefficient. Hence, the total
plasma density induced by SFI is given by:

$\displaystyle\rho_{\rm{SFI}}(t)$
$\displaystyle=a_{0}\int_{-\infty}^{t}\frac{B(t^{\prime})}{2}dt^{\prime}+\int_{-\infty}^{t}\frac{B(t^{\prime})}{2}\sum_{n=1}^{\infty}a_{n}(t)\cos(2n\omega_{0}t^{\prime})dt^{\prime}$
$\displaystyle\approx\underbrace{\int_{-\infty}^{t}\frac{B(t^{\prime})}{2}dt^{\prime}}_{\rm{cycle-
averaged}}+\underbrace{\frac{B(t^{\prime})}{2}\sum_{n=1}^{\infty}\int_{-\infty}^{t}a_{n}(t)\cos(2n\omega_{0}t^{\prime})dt^{\prime}}_{\rm{sub-
cycle}}.$ (10)

The decomposition of the plasma density, the SFI rate and the strong-field-
induced conductivity into a cycle-averaged and a sub-cycle component is
depicted in Fig. S3. While the cycle-averaged contributions have only positive
values, the sub-cycle components can be both positive or negative.
As can be seen from Eq. IV.4, all harmonic orders of the two-color wave-mixing
signal exhibit the same temporal signature given by $B(t)$. This is supported
by the identical temporal profile of the detected harmonic peaks (see Fig. S4)
and allows the extraction of the cycle-averaged, strong-field-induced plasma
dynamics from an arbitrary order of the injection harmonics.

Figure S4: Projection of the harmonic peaks from Fig. 1(b) on the pump-probe
delay axis compared to the pump-probe convolution.

###  IV.5 Test of the reconstruction procedure with the help of numerical
simulations

Figure S5: Numerical test of the reconstruction procedure. (a). Simulated
time-frequency map of the two-color wave-mixing signal generated by the
injection current. (b). Retrieved spectrogram by the ptychographic iterative
engine after 1000 iterations. (c) Direct output of the phase retrieval
algorithm for $\sigma_{\rm{inj}}(t)$. (d). Corrected effective strong-field-
induced conductivity after levelling the spectral phase. (e). Cycle-averaged
and sub-cycle ionization rate computed from the corrected
$\sigma_{\rm{inj}}(t)$. (f). Total SFI rate $\dot{\rho}_{\rm{SFI}}(t)$
obtained from a ptychographic analysis of simulated, two-color injection
harmonics. (g). Resulting plasma dynamics induced by SFI compared to the
numerical input data.

Figure S5 shows a numerical test of the reconstruction procedure using
numerical simulations of the optical signature of the injection current
according to Eq. III. In the numerical model an instantaneous Ammosov-Delone-
Krainov (ADK) tunneling rate Ammosov  _et al._ (1986) is employed to calculate
$\dot{\rho}_{\rm{SFI}}(t)$. The resulting time-resolved spectrum
$Y_{\rm{inj}}(\omega,\tau)$ containing the first eight wave-mixing orders is
displayed in Fig. S5(a). Figure S5(b) shows the reconstructed spectrogram of
the ptychographic iterative engine after 1000 iterations. The spectral
features of the input data are accurately reproduced and lead to the retrieved
effective conductivity $\sigma_{\rm{inj}}(t)$ shown in Fig. S5(c). Analogous
to the case of the experimental harmonic spectra the harmonic orders do not
spectrally overlap and therefore the phase retrieval algorithm is not able to
determine the relative phase of the harmonics. In the case of simulated data
where the absolute and relative phases of all fields are well-known and the
SFI process is assumed to be instantaneous we can manually set the spectral
phase over the whole frequency range to 0. This leads to the sub-cycle
modulated effective conductivity displayed in Fig. S5(d).
The resulting $\sigma_{\rm{inj}}^{\rm{sub}}(t)$ is used to reconstruct the
sub-cycle and the cycle-averaged component of the SFI rate shown in Fig. S5(e)
while the corresponding total SFI rate is depicted in Fig. S5(f). Direct
integration of the obtained $\dot{\rho}_{\rm{SFI}}(t)$ leads to the expected
stepwise increase of the plasma density [$\rho_{\rm{SFI}}(t)$ in Fig. S5(g)]
that is compared to the cycle-averaged plasma formation and the numerical
input of the simulation. The numerical input of the simulation is reproduced
with a very high precision demonstrating that the reconstruction of
$\rho_{\rm{SFI}}(t)$ from a time-resolved measurement of low-order injection
harmonics is indeed feasible. Note that the shift between the simulated and
reconstructed curve results from the CEP-insensitivity of the retrieval
algorithm that is accounted for by the uncertainty of one optical cycle
indicated by the shaded area in Fig. S5(g). Hence, the presented experimental
strategy enables the reconstruction of sub-cycle ionization dynamics with
attosecond precision if the relative phase of the harmonic orders is known.

Figure S6: Influence of the maximum harmonic order $n$ on the reconstruction
of $\rho_{\rm{SFI}}(t)$ from simulated, time-resolved spectra (results for
$n\geq 2$ are vertically shifted for visibility). Inset: Standard deviation
between the individual reconstructed densities and the numerical input data.

Analysis of the influence of the maximum harmonic order $n_{\rm{max}}$ on the
reconstruction of the SFI-induced plasma density is analyzed in Fig. S6.
Herein, the plasma densities reconstructed from more than the first harmonic
order are shifted vertically for better visibility. The inset shows the
standard deviation of the individual reconstructed plasma densities with
respect to the numerical input data. These results indicate that the
reconstruction of the sub-cycle dynamics of SFI benefits from an increased
number of harmonics inserted into the retrieval algorithm. However, orders
with $n>1$ only lead to minor modifications of the reconstructed ionization
dynamics. Hence, analyzing the $n=1$ harmonic only already gives valuable
insights into the strong-field-driven plasma dynamics while the exact shape of
the nearly-stepwise increase of $\rho_{\rm{SFI}}(t)$ is encoded in the higher
harmonic orders.

###  IV.6 Case of dominating Brunel harmonics

If in contrast to the above discussion the Brunel current instead of the
injection current becomes the dominant source of low-order harmonic
generation, it is still possible to reconstruct the strong-field-induced
plasma formation from a time-resolved measurement of the two-color harmonics.
The expression for the current density given in Eq. III needs to be replaced
by the Brunel current $j_{\rm{Br}}(t)$ Brunel (1990); Babushkin  _et al._
(2017):

$\frac{\partial j_{\rm{Br}}(t,\tau)}{\partial
t}=\frac{q_{e}^{2}}{m_{e}}\rho_{\rm{SFI}}(t)E_{\rm{probe}}(t-\tau)$ (11)

Again, the above equation is only valid for cross-polarized pump and probe
laser fields where the ionization is solely driven by the strong pump laser
field. In analogy to the derivation of Eq. III the resulting spectrogram can
be written as

$\displaystyle Y_{\rm{Br}}(\omega,\tau)$
$\displaystyle=\left|\int\frac{q_{e}^{2}}{m_{e}}\rho_{\rm{SFI}}^{\rm{sub}}(t)E_{\rm{probe}}(t-\tau)e^{i\omega
t}dt\right|^{2}$ (12)

Hence, phase-retrieval analysis of $Y_{\rm{Br}}(\omega,\tau)$ gives direct
access to the strong-field-induced, sub-cycle variations of the plasma density
$\rho_{\rm{SFI}}^{\rm{sub}}(t)$. The total plasma density due to SFI can be
obtained in an analogous fashion to the previously described case of
dominating injection harmonics.


# Do We Run Large-scale Multi-Robot Systems on the Edge?
More Evidence for Two-Phase Performance in System Size Scaling

Jonas Kuckling,1,2 Robin Luckey,3 Viktor Avrutin,4 Andrew Vardy,5
Andreagiovanni Reina,2,6 Heiko Hamann1,2 The work of VA was supported by the
German Research Foundation within the scope of the project “Generic
bifurcation structures in piecewise-smooth maps with extremely high number of
borders in theory and applications for power converter systems – 2”. AR
acknowledges support from the Belgian F.R.S.-FNRS, of which he is Chargé de
Recherches.1JK and HH are with the Department of Computer and Information
Science, University of Konstanz, 78457 Konstanz, Germany. {jonas.kuckling,
<EMAIL_ADDRESS>JK, AR, and HH are also with the Centre for
the Advanced Study of Collective Behaviour, University of Konstanz,
Germany.3RL is with the Institute of Computer Engineering, University of
Lübeck, Lübeck, Germany.4VA is with the Institute for Systems Theory and
Automatic Control, University of Stuttgart, Stuttgart, Germany.5AV is with the
Department of Computer Science and the Department of Electrical and Computer
Engineering, Memorial University of Newfoundland, St. John’s, Canada.6AR is
with IRIDIA, Université Libre de Bruxelles, Brussels, Belgium.

###### Abstract

With increasing numbers of mobile robots arriving in real-world applications,
more robots coexist in the same space, interact, and possibly collaborate.
Methods to provide such systems with system size scalability are known, for
example, from swarm robotics. Example strategies are self-organizing behavior,
a strict decentralized approach, and limiting the robot-robot communication.
Despite applying such strategies, any multi-robot system breaks above a
certain critical system size (i.e., number of robots) as too many robots share
a resource (e.g., space, communication channel). We provide additional
evidence based on simulations, that at these critical system sizes, the system
performance separates into two phases: nearly optimal and minimal performance.
We speculate that in real-world applications that are configured for optimal
system size, the supposedly high-performing system may actually live on
borrowed time as it is on a transient to breakdown. We provide two modeling
options (based on queueing theory and a population model) that may help to
support this reasoning.

##  I Introduction

To scale up to larger and larger problems, robot tasks in industry,
environmental/infrastructure monitoring, and our everyday lives will require
larger and larger multi-robot systems, acting in parallel and cooperating [1,
2]. With an increasing number of robots sharing resources, such as physical
space and radio channels, scalability in system size will become an increasing
challenge [3]. Methods for designing robot systems that scale (at least for
fixed swarm density) are developed and studied in swarm robotics [4]. Even
robots that were not designed to interact and cooperate may be required to do
so in the future—further increasing the necessity to develop scalable systems.
[5].

Swarm robotics achieves scalability by enforcing a decentralized approach
based on self-organizing behaviors, local robot-robot interactions, and simple
neighbor-based communication [4]. Several studies have measured a certain
intuitive swarm performance (e.g., number of collected items during a defined
time interval in a foraging task [6]), usually over swarm density $\rho=N/A$,
for swarm size $N$ and a bounded operation area of size $A$ [7, 8, 9, 3]. If
area size $A$ is fixed, scaling swarm size $N$ effectively modulates the swarm
density. Almost all studies that measure mean swarm performance over swarm
density (e.g., see [10]) show similar curves characterized by three stages: an
initial increase with system size, reaching peak performance (at a critical
system size $N_{c}$), and for $N>N_{c}$ decreasing system performance.
Interestingly, a diversity of systems from different fields of research show
similar scalability results (e.g., see [11, 12, 13]). While measuring mean
swarm performance is a good approach for an initial performance analysis, it
is likely to hide an important property. The measured swarm performance may be
bimodal showing two phases: an upper phase with high performance (i.e., each
robot operates efficiently most of the time), a lower phase with minimal
performance (i.e., most robots trapped in congestion etc.), and almost nothing
in between [14]. This is, however, a potentially dangerous condition as swarm
robotics systems often show linear or even super-linear performance increases
for low swarm densities. By optimizing for performance, one may unknowingly
put the system on the edge, where model predictions and real-world experience
may catastrophically diverge. In addition, the robotic system is likely to
have a long transient at the critical system size which may have two
implications. First, the system can actually be run on ‘borrowed time’ if the
transient is longer than what is relevant in applications (i.e., potentially
contradicting mathematical predictions built on assumed convergence). Second,
the system at runtime initially shows good performance but then suddenly
breaks down completely.

As our main contribution, we show simulation results for a large real-world
multi-robot system in a warehouse and two more, very different, academic and
smaller scenarios (swarm robotics object clustering and emergent taxis)
indicating clearly the two-phase situation. In addition, we provide two
modeling options that are able to qualitatively describe this crucial two-
phase system behavior. Which model may prove to be best for modeling two-phase
swarm performance remains an open question at this point. It remains
encouraging that different modeling techniques are able to represent the two-
phase situation and yield quite similar predictions. Besides the obvious
conclusion that a swarm engineer wants to avoid putting their robotic systems
on the edge between performing maximally and breaking down completely, there
are other immediate implications. For example, hysteresis effects are possible
that may have extreme effects on the swarm performance at runtime. Known
approaches of ‘robust scalability,’ such as the online voluntary retreat of
robots in order to reduce swarm density [15] or of adapting algorithm
parameters autonomously in reaction to a currently measured swarm density
[16], may trigger sudden performance jumps if the system operates close to a
bifurcation point (see Sec. IV-B).

##  II Three Scenarios

###  II-A Warehouse

We consider an actual modern warehouse solution as required, for example, by
today’s online retailers of fashion to maximize economy of scale. The system
is based on robots called ‘AFLE’ (German abbreviation for Autonomous Driving
Storage Units) that are currently under development by the company EMHS GmbH
(see Fig. 1(a)). The idea is to build a large swarm (in the order of $10^{5}$)
of low-cost mobile autonomous ‘hangers’ as a solution to traditional pocket
sorters instead of a dragging chain. Each AFLE carries a bag forming a storage
unit. The robots hang on a rail grid that also provides power. By dividing the
warehouse into independent segments the system has high potential to scale. To
further enforce scalability, the system could combine central coordination by
global messaging with decentralized approaches to coordinate pairs or small
groups of AFLEs, for example, at crossings of rails. Hence, the system at
least partially qualifies as a swarm robotics system.

The simulation software is focused on large-scale simulations of warehouses
with more than 3000 AFLEs. The warehouse is simplified to a grid world with
grid cells that can have an edge per cardinal direction and contain not more
than one robot. All traffic is one-way and organized in ‘highways’; the actual
storage rail sections are not modeled. Robots emerge on highways modeled by
Poison distributions (i.e., mimicking the rather unpredictable influx of
orders). The AFLE robots are routed with $A^{*}$. No physics is simulated but
we penalize robots doing turns as the real robots would need to slow down.

###  II-B Object Clustering

In the object clustering task a swarm of robots manipulate objects, either
bringing them to a pre-defined goal position, or aggregating them in a self-
organized manner (see Fig. 1(b)). If the goal is pre-defined, then this task
becomes synonymous with foraging [17]. We test the _lasso method_ for object
clustering as its performance as a function of swarm size has already been
partially explored [18]. This method relies on an external scalar field which
acts to guide the robots around the environment while nudging objects towards
the global minimum of the field, located at the goal.

The control algorithm operates independently on each robot. The most outlying
visible puck is identified and the value of this puck’s position on the scalar
field is determined. A control law drives the robot to reach and follow the
contour line defined by this value. We observe that the robots tend to occupy
the same contour line and travel clockwise around it, nudging pucks inwards
toward the goal. Thus, a ‘lasso’ is formed that tightens around the objects
until they are forced to the goal position. Performance was found to increase
with swarm size up to a point, then decrease due to spatial interference [18].
The results depicted in Fig. 1(e) were collected in a custom simulator which
models two-dimensional physics111open-source code repository:
https://github.com/avardy/cwaggle_lasso, video: https://youtu.be/_KOU5SzpQBg..

###  II-C Emergent Taxis

In the emergent taxis task [14], the objective is to move the swarm towards a
light beacon (see Fig. 1(c)). However, the robots only have an omnidirectional
beacon sensor to identify the beacon. That is they can only perceive if their
line of sight to the beacon is unoccupied but no bearing information is
available. If a robot has a line of sight to the beacon, it is referred to as
_illuminated_ otherwise it is considered _shadowed_. The robots move straight,
until they approach a defined distance to another robot and trigger an
obstacle avoidance behavior. This _avoid radius_ depends on the state of the
robot: if they are shadowed, the radius is smaller than if they are
illuminated. Thus, robots are biased to move towards the light beacon.
Additionally, a parameter $\alpha$ controls the coherence of the swarm.

With increasing swarm size, the effect of increased number of encounters
influences the swarm performance. Initially, the increased number of robots
increases the density with which shadowed robots form behind and approach
illuminated ones. This in turn increases the number of times that an
illuminated robot avoids a shadowed one and thus the movement bias towards the
light beacon becomes stronger. However, if the required coherence forces the
robots to aggregate too tightly, the obstacle avoidance behavior might be
continuously triggered thus slowing down any movement. Previous results showed
that performance increased with swarm size but a second phase of low
performance emerged at higher swarm densities [14].

##  III Results

Next, we present simulation results for each of the three scenarios. We focus
on swarm performance over swarm size $N$. Swarm performance is throughput for
the warehouse scenario (i.e., how many autonomous hangers reach their
destination within a given total experiment time), average distance between
objects for the object clustering scenario, and the speed of the swarm’s
barycenter towards the light for the emergent taxis scenario. As we keep the
operation area size $A$ constant in all three cases, we effectively change the
swarm density $\rho$ by varying swarm size $N$. In the case of the warehouse,
we vary the swarm size $N$ of coexisting robots indirectly via an arrival rate
of how many robots are requested to approach the packing station per hour.

###  III-A Warehouse

(a) Warehouse prototype view (b) Object clustering overhead view (c) Emergent
taxis schematic

(d) Warehouse results (e) Object clustering results (f) Emergent taxis results

Figure 1: Upper row (a-c): Representative illustration of the three scenarios.
Lower row (d-f): histograms of the swarm performance in simulation over swarm
size (or arrival rate in the case of the warehouse scenario).

Fig. 1(d) shows a 2d-histogram of the performance in simulation in the
warehouse scenario. We investigated different arrival rates with which new
AFLEs are added to the warehouse. The performance is measured as throughput,
that is the number of AFLEs that reach their destination. As AFLEs arrive and
leave the warehouse during the experiment, the swarm size fluctuates over
time. However, assuming a constant time that the AFLEs spend inside the
warehouse, the arrival rate can be an approximation for the number of AFLEs in
the warehouse. For each investigated value of the arrival rate, 96 repetitions
were performed. The color value of a cell in the histogram indicates how often
a given arrival rate resulted in the corresponding throughput.

At low arrival rates, the system does not suffer from congestion and the
throughput increases linearly with the arrival rate, as expected for a system
with low density. At an arrival rate of about $10^{4}$ units per hour
interference effects lead to some experiments in which the throughput drops
considerably. This drop of throughput is caused by a congestion in the
warehouse, when too many AFLEs interfere with each other. When two AFLEs are
trying to traverse the same rail segment, they need to negotiate the traversal
order and AFLEs with lower priority wait until AFLEs with higher priorities
are passed. With increased arrival rate, the number of experiments with
congestion increases until all repetitions suffer from congestion. It is
important to note that there is a range of arrival rate values for which part
of the experiments operated properly, with maximum performance, and the
remaining experiments collapsed into complete congestion with performance
close to zero. As a result, the distribution of the throughput splits into two
_phases_ , one high performing and one low performing.

###  III-B Object Clustering

Figure 2: The performance plot from Fig. 1(e) augmented with snapshots of the
final arena configuration. The robots are displayed in red. Pucks are shown in
green. The goal position is depicted at top-left.

Fig. 1(e) shows the performance of simulated swarms performing the object
clustering task for increasing swarm sizes, from one to twenty robots. For
each swarm size, we evaluated the performance $1000$ times. Initially, all
robots and objects were randomly distributed in the arena. The performance is
the normalized average distance of all objects to the target location. A
performance measure of $0$ is the ideal, but is not attainable due to physical
interference between objects.

When the number of robots is too small, a small number of objects is nudged
towards the goal, resulting in poor collective performance. Performance
increases with increasing swarm size, as predicted by theory with optimal
performance around 6 robots. For swarms with 10 or more robots, we observed
that the performance splits into two phases—roughly, successful and
unsuccessful. Fig. 2 provides examples of the final arena configuration with
16 robots for different levels of performance. The top example shows a
successful outcome in which all objects have been delivered to the goal. The
remaining examples represent situations where congestion has occurred and
severely reduced performance.

###  III-C Emergent Taxis

Fig. 1(f) shows the performance in simulation of the swarm in the emergent
taxis task. This result was published before and data was taken directly from
the previous publication [14]. Each swarm size of $N\in\\{2,3,\dots,40\\}$
robots was evaluated 800 times. At the beginning of the experiment, robots
were distributed randomly on the far side of the arena. The score of each
evaluation was the speed with which the barycenter of the robot swarm moved
towards the light during the experiment. For negative speeds, the barycenter
of the swarm moved away from the light.

For small swarm sizes, we did not observe any emergent taxis behavior.
However, with increasing swarm size, the swarm started moving towards the
light and the speed of this taxis behavior also increased. At around a swarm
size of 15 robots, the maximum speed was reached and no further improvement
was found. At the same swarm size, we observed the emergence of a second
phase, in which only the observed speed of the taxis behavior was low or
sometimes nonexistent. In the first phase, the swarm is still able to function
properly, whereas in the second case it becomes _pinned_ , with many obstacle
avoidance maneuvers triggered resulting in only slow movements.

##  IV Models

Next, we propose two modeling approaches that capture the two-phase character
of swarm performance as observed in real-world systems. The first is a simple
method based on queueing theory as an example of low model complexity from
computer science that is still capable of describing key features. As second
method we use our previously published population model [3] that we introduced
as a generic approach to model multi-robot systems but also large parallel
computing systems. It comes with a lot of complexity but also proves to be
capable of catching key properties of the swarm systems including predictions
for transient times, that is for example, how long a real-world robot system
would perform high despite sitting on the edge. We provide a full discussion
of the resulting nonlinear system in terms of nonlinear dynamics. We believe
that especially the model predictions of transient times may prove to be
relevant for real-world applications. Hence, a deep understanding of all the
mathematical implications is essential for a complete (future) understanding
of scalability in robot swarms.

###  IV-A Queues

A simple approach of modeling the two-phase situation of swarm performance can
be based on queueing theory. An intuitive and sophisticated model could be
based on queueing networks to explicitly represent the shared resources, such
as space and radio bandwidth (possibly also truncated queues). However, a
simpler approach is possible by an unorthodox change of the standard M/M/1
queue with exponential interarrival-time distribution $\lambda\exp{(-\lambda
t)}$, exponential service time $\mu\exp{(-\mu t)}$, one channel, infinite
capacity. Key change to the standard model is that we make the service time
dependent on the current queue length $N_{q}(t)$ (i.e., a cost on keeping the
queue) by setting mean service time
$1/\mu(t)=1/(c_{1}\log(c_{2}N_{q}(t))+\mu_{c})$ for arbitrary constants
$c_{1}=0.3$ and $c_{2}=0.2$. For the following example we set $\mu_{c}=1$. We
simulate the queue for mean arrival rates $\lambda\in[0.4,0.9]$ (500 samples
each) for $2\times 10^{6}$ time steps and exclude data from the first $2\times
10^{6}$ time steps as transient. Resulting mean throughput and mean queue
length are shown in Fig. 3(a) and a histogram of the throughput in Fig. 3(b).
It is to note that the curves shown in Fig. 3(a) resemble those predicted by
the population model (see Fig. 5a). This queuing model is able to represent
the linear increase and the breakdown of the system at a critical arrival rate
$\mu_{c}\approx 0.65$.

(a) throughput and queue length $N_{q}$

(b) throughput histogram (frequency color-code log.)

Figure 3: Queuing model: throughput and queue length $N_{q}$ over arrival rate
for modified M/M/1 queue (service time dependent on queue length).

###  IV-B Population model and nonlinear dynamics analysis

In previous work, we have presented a population model that grasps most of the
relevant system features concerning scalability [3]. In this alternative
approach to model scalability we assume that each robot can be assigned to one
of three abstract robot states relevant for scalability analysis. The system
dynamics is based on $N$ robots that transit between these three operational
states. We assume robots can be in state _solo_ (operating alone without
sharing communication and/or space resources with others), in state _grupo_
(sharing resources with other robots), or in state _fermo_ (interfering with
others in the use of resources). The model describes through seven transition
rates ($k_{1},\dots,k_{7}$) the probability that robots change their state as
a consequence of independent actions or interactions with other robots. As in
any population model we make assumptions, for example, that it is more likely
for a _solo_ robot to approach a _solo_ robot if the _solo_ state currently
has a high percentage in the swarm. We assume also that two _solo_ robots
approaching each other trigger transitions to _grupo_ robots and similar for
any other interaction of two robots (see Fig. 4).

Figure 4: Underlying schema for population model: transitions for robots in
state _solo_ , _grupo_ , and _fermo_ [3].

The model is expressed as a system of ordinary differential equations (ODEs):

$\left\\{\begin{aligned} \frac{ds}{dt}&=-2k_{1}s^{2}-k_{2}sg-k_{3}sf+k_{4}g\\\
\frac{df}{dt}&=2k_{5}g^{2}+k_{6}gf-k_{7}f\,,\end{aligned}\right.$ (1)

with $g=N-s-f$ due to the robot conservation assumption. The system variables
$s$, $g$, and $f$ represent the proportion of robots in states _solo_ ,
_grupo_ , and _fermo_ , respectively. $(s^{*},g^{*},f^{*})$, the stable fixed
point of system (1), describes the long-term distribution of the robots among
the three states.

(a) (b) (c)

Figure 5: Stationary states of Eq. (1) under variation of $N$ for (a)
$k_{1}=0.005$ and (b) $k=0.001$, other parameters as in Eq. (2). (b) Domain
between the saddle-node bifurcation points $N^{(1)}_{\text{\tiny SN}}$,
$N^{(2)}_{\text{\tiny SN}}$ associated with bi-stability, the non-physical
domain outside of $0\leq s,f,g\leq N$ is shown in gray; (c) transient time
$T_{\text{trans}}$ for fixed initial state $(s_{0},f_{0})=(0,0)$ and
convergence to stationary state $(s_{1}^{*},f_{1}^{*})$ with accuracy
$10^{-12}$; $k_{1}=0.005$. Insets show the indicated rectangles magnified.
Figure 6: _fermo_ -_solo_ state space in the bi-stable domain ($N=40$); color-
coded logarithmic transient time $T_{\text{trans}}$, that is, convergence
towards stable stationary states $(s_{1}^{*},f_{1}^{*})$ and
$(s_{3}^{*},f_{3}^{*})$ with the same accuracy of $10^{-12}$. Also the saddle
point $(s_{2}^{*},f_{2}^{*})$ and its stable and unstable manifolds
($W^{s}_{\pm}$ and $W^{u}_{\pm}$, respectively) are shown. The inset shows the
transient times for initial values located at the path indicated with A.

In the numerical simulations presented below, we use the following set of
parameters, unless stated otherwise:

$\begin{split}&k_{1}=0.005,\>k_{2}=0.1,\>k_{3}=0.06,\>k_{4}=10,\\\
&k_{5}=0.15,\>k_{6}=0.3,\>k_{7}=0.8\end{split}$ (2)

For a quick understanding of this model, we have analyzed the stationary
system behavior. Fig. 5a shows the stationary states $(s^{*},f^{*})$ of system
(1) at the parameter set specified in Eq. (2) under variation of swarm size
$N$. For completeness, the resulting value $g^{*}=N-s^{*}-f^{*}$ due to robot
conservation is shown, too. Fig. 5a shows, how for small swarm sizes $N$, the
number of _solo_ robots $s^{*}$ grows almost linearly (most robots remain in
_solo_ , few in _fermo_ and _grupo_). For increasing swarm size $N$, a sudden
change of the system’s configuration occurs. The number of _solo_ robots
reaches a maximum (here, $N=31$) and then decrease rapidly. Simultaneously, we
get more _fermo_ robots, resulting eventually in their predominance. From the
perspective of nonlinear dynamics, the shifted location of the asymptotically
stable fixed point is not a bifurcation as the topological structure of the
state space remains unchanged under parameter variation. However, in any
practical application the sudden decrease of _solo_ robots would have a
significant impact on system performance.

For other parameter values we may observe more sophisticated dynamics. As an
example, Fig. 5b shows the stationary states of the system for the value
$k_{1}=0.001$ instead of $k_{1}=0.005$ (i.e., slower transition of _solo_
robots to _grupo_ robots). As before, we observe an almost linear growth of
_solo_ robots until it reaches a maximal value. However, soon after, the
stable solution $(s_{1}^{*},f_{1}^{*})$ merges with an unstable one
$(s_{2}^{*},f_{2}^{*})$ and disappears in a saddle-node bifurcation at
$N=N^{(2)}_{\text{\tiny SN}}$. The unstable stationary state
$(s_{2}^{*},f_{2}^{*})$ involved in this bifurcation appears already before,
via a different saddle-node bifurcation which occurs at
$N=N^{(1)}_{\text{\tiny SN}}<N^{(2)}_{\text{\tiny SN}}$ and gives rise to
another asymptotically stable stationary state $(s_{3}^{*},f_{3}^{*})$. Thus,
in the parameter interval $N^{(1)}_{\text{\tiny SN}}<N<N^{(2)}_{\text{\tiny
SN}}$, two stable stationary states $(s_{1}^{*},f_{1}^{*})$,
$(s_{3}^{*},f_{3}^{*})$ coexist and the system exhibits a bi-stable behavior,
so that it may converge either to the stationary state where _solo_ robots
prevail, or to the one characterized by the prevalence of _fermo_ robots,
depending on initial values $s_{0}$ and $f_{0}$.

In Fig. 6, we illustrate the system’s state space. As two saddle-node
bifurcations occurring at $N=N^{(1)}_{\text{\tiny SN}}$ and
$N=N^{(2)}_{\text{\tiny SN}}$ are connected by the unstable (saddle) point
$(s_{2}^{*},f_{2}^{*})$, the basins of attraction of the stable stationary
states are separated from each other by the branches $W^{s}_{\pm}$ of the
stable manifold of the saddle $(s_{2}^{*},f_{2}^{*})$. Here, the initial
values located above or below $W^{s}_{\pm}$ converge towards
$(s_{1}^{*},f_{1}^{*})$ or $(s_{3}^{*},f_{3}^{*})$, respectively.

Importantly for any application, the bi-stable dynamics lead necessarily to a
hysteresis effect. If we start from small swarm sizes $N$ and increase swarm
size incrementally, then the system follows the upper branch of the solutions,
that is, the stationary state $(s_{1}^{*},f_{1}^{*})$ (see Fig. 5b). At the
bifurcation point $N=N^{(2)}_{\text{\tiny SN}}$, the system undergoes a sudden
jump to the lower branch $(s_{3}^{*},f_{3}^{*})$ (the saddle-node bifurcation
acts as a so-called hard bifurcation, or a hard loss of stability [19]). If
one additional robot is introduced ($N+1$), the system undergoes a sudden
change, which corresponds to a sudden catastrophic loss of productivity. This
roughly resembles the Braess paradox from game theory [20] where the overall
performance of a system may decrease as an additional resource is provided.
Instead, if we start from a large swarm size $N$ and we switch off robots one
by one, then the system follows the lower branch of solutions (stationary
state $(s_{3}^{*},f_{3}^{*})$) operating at low productivity until the
bifurcation point $N=N^{(1)}_{\text{\tiny SN}}$. Then, the system jumps to the
upper branch $(s_{1}^{*},f_{1}^{*})$ and the productivity is suddenly
increased.

Important for all applications is the transient time (i.e., time to
converge—with a given accuracy—to the stationary state). As an example, Fig.
5c shows the transient time on logarithmic scale corresponding to the
bifurcation diagram in Fig. 5a. The transient time grows over-exponentially as
the system approaches the critical state $N_{\text{\tiny C}}$. This behavior,
known as _critical slowing down_ , is typically observed in a neighborhood of
a bifurcation point. This may be surprising, since, as already mentioned,
there is no bifurcation (in the rigorous sense) in Fig. 5a. However, in Fig.
5b (different parameter $k_{1}$) we observe two saddle-node bifurcations.
Therefore, in the 2D parameter space $(N,k_{1})$ there exists a codimension-2
_cusp_ -bifurcation point [21] from which two saddle-node bifurcation points
emerge, leading to bistability and hysteresis. Hence, the long transient time
close to the critical state (see Fig. 5a) is indeed caused by the critical
slowing down in a neighborhood of a cusp bifurcation in 2D parameter space
even if this bifurcation is not visible in the bifurcation diagram shown in
Fig. 5a. In addition, the transient time $T_{\text{trans}}$ depends not
exclusively on the parameters but also on the initial values $s_{0}$, and
$f_{0}$. The transient time grows (over-exponentially, tending to infinity)
for initial values close to the boundary between the basins of attraction of
the coexisting stable fixed points (see Fig. 6).

##  V Discussion and Conclusion

We have presented three multi-robot systems that show the characteristic two-
phase performance signature where, for a critical swarm size, the phases of
high and low performance coexist. We have presented that for simulations of
the AFLE robot system by EMHS GmbH. Hence, our findings clearly show that not
only real-world implementations of robot swarms but also real-world
applications of actual productive systems in existing industries show this
same property in their system size scaling (or in this case scaling with
arrival rates). We have also presented simulation results of an object
clustering scenario that was previously implemented on real robots, showing
the validity of the simulation [18]. Also in this more academic example from
swarm robotics, we observe the same two-phase performance signature. As a last
example, we have reminded the reader of simulation results for the well-known
emergent taxis scenario that were published before [14]. This original finding
inspired the presented research. To point to a possible interpretation of
these results, we have presented two modeling options based on queueing theory
and population models. The adapted M/M/1 queue is able to catch most of the
discussed scalability features despite its simplicity. Our key contribution in
terms of modeling, however, is the population model and especially its
interpretation based on nonlinear dynamics. These interpretations point to
three important findings and issues. (a) There is indeed a bifurcation in the
system and rapid changes in performance with small changes in swarm size
should be expected. (b) In the vicinity of bifurcation points we should indeed
expect extremely increased transient times. (c) We should expect hysteresis
effects depending on swarm size that is increased or decreased online while
operating. Swarm robotics systems that are optimized for swarm performance may
unknowingly be put on the edge with the effect of unexpected, catastrophic
performance breakdowns. Long transients may either obscure this effect or may
buy the robot swarm time that is possibly long enough to survive in high-
performing states for the whole mission time. Future research will hopefully
provide more evidence of whether the two-phase performance signature in
scalability of swarm robotics is real and relevant even in industrial
applications.

## References

  * [1] G.-Z. Yang, J. Bellingham, P. E. Dupont, P. Fischer, L. Floridi, R. Full, N. Jacobstein, V. Kumar, M. McNutt, R. Merrifield, B. J. Nelson, B. Scassellati, M. Taddeo, R. Taylor, M. Veloso, Z. L. Wang, and R. Wood, “The grand challenges of Science Robotics,” _Science Robotics_ , vol. 3, no. 14, p. eaar7650, 2018.
  * [2] M. Dorigo, G. Theraulaz, and V. Trianni, “Swarm robotics: Past, present, and future [point of view],” _Proceedings of the IEEE_ , vol. 109, no. 7, pp. 1152–1165, 2021.
  * [3] H. Hamann and A. Reina, “Scalability in computing and robotics,” _IEEE Transactions on Computers_ , vol. 71, no. 6, pp. 1453–1465, 2022.
  * [4] H. Hamann, _Swarm Robotics: A Formal Approach_. Springer, 2018.
  * [5] I. Rahwan, M. Cebrian, N. Obradovich, _et al._ , “Machine behaviour,” _Nature_ , vol. 568, p. 477–486, 2019.
  * [6] Q. Lu, G. Fricke, J. Ericksen, and M. Moses, “Swarm foraging review: Closing the gap between proof and practice,” _Current Robotics Reports_ , vol. 1, pp. 215–225, 2020.
  * [7] H. Kwa, J. Philippot, and R. Bouffanais, “Effect of swarm density on collective tracking performance,” _Swarm Intelligence_ , vol. 17, p. 253–281, 2023.
  * [8] E. R. Hunt, S. Jones, and S. Hauert, “Testing the limits of pheromone stigmergy in high-density robot swarms,” _Royal Society Open Science_ , vol. 6, p. 190225, 2023.
  * [9] Y. Khaluf, C. Pinciroli, G. Valentini, and H. Hamann, “The impact of agent density on scalability in collective systems: noise-induced versus majority-based bistability,” _Swarm Intelligence_ , vol. 11, no. 2, pp. 155–179, Jun 2017. [Online]. Available: https://doi.org/10.1007/s11721-017-0137-6
  * [10] A. Rosenfeld, G. A. Kaminka, and S. Kraus, “A study of marginal performance properties in robotic groups,” in _Proceedings of the Third International Joint Conference on Autonomous Agents and Multiagent Systems - Volume 3_ , ser. AAMAS ’04. USA: IEEE Computer Society, 2004, p. 1536–1537.
  * [11] J. Ward and R. T. Compton, “High throughput slotted ALOHA packet radio networks with adaptive arrays,” _IEEE Transactions on Communications_ , vol. 41, no. 3, pp. 460–470, 1993.
  * [12] D. Mateo, Y. Kuan, and R. Bouffanais, “Effect of correlations in swarms on collective response,” _Scientific Reports_ , vol. 7, p. 10388, 2020.
  * [13] B. S. Kerner, _Introduction to modern traffic flow theory and control_. Springer, 2009.
  * [14] H. Hamann, “Towards swarm calculus: Urn models of collective decisions and universal properties of swarm performance,” _Swarm Intelligence_ , vol. 7, no. 2-3, pp. 145–172, 2013. [Online]. Available: http://dx.doi.org/10.1007/s11721-013-0080-0
  * [15] S. Mayya, P. Pierpaoli, and M. Egerstedt, “Voluntary retreat for decentralized interference reduction in robot swarms,” in _International Conference on Robotics and Automation_. IEEE, 2019, pp. 9667–9673.
  * [16] M. Wahby, J. Petzold, C. Eschke, T. Schmickl, and H. Hamann, “Collective change detection: Adaptivity to dynamic swarm densities and light conditions in robot swarms,” _Artificial Life Conference Proceedings_ , no. 31, pp. 642–649, 2019.
  * [17] O. Zedadra, N. Jouandeau, H. Seridi, and G. Fortino, “Multi-agent foraging: state-of-the-art and research challenges,” _Complex Adaptive Systems Modeling_ , vol. 5, no. 1, pp. 1–24, 2017.
  * [18] A. Vardy, “The lasso method for multi-robot foraging,” in _19th Conference on Robots and Vision (CRV)_ , 2022, pp. 106–113.
  * [19] V. Arnol’d, _Catastrophe theory_. Springer Science & Business Media, 2003.
  * [20] D. Braess, “Über ein Paradoxon aus der Verkehrsplanung,” _Unternehmensforschung_ , vol. 12, pp. 258–268, 1968.
  * [21] Y. Kuznetsov, _Elements of Applied Bifurcation Theory_. Springer, 2004, 3rd edition.