# Collisions enhance tracer diffusivity in an active Janus bath

Maximilian R. Bailey  Laboratory for Soft Materials and Interfaces, Department
of Materials, ETH Zürich, Vladimir-Prelog-Weg 5, 8093 Zürich, Switzerland
Dmitry A. Fedosov  Theoretical Physics of Living Matter, Institute of
Biological Information Processing and Institute for Advanced Simulation,
Forschungszentrum Jülich, 52425 Jülich, Germany  Federico Paratore  Laboratory
for Soft Materials and Interfaces, Department of Materials, ETH Zürich,
Vladimir-Prelog-Weg 5, 8093 Zürich, Switzerland  Fabio Grillo  Laboratory for
Soft Materials and Interfaces, Department of Materials, ETH Zürich, Vladimir-
Prelog-Weg 5, 8093 Zürich, Switzerland  Gerhard Gompper  Theoretical Physics
of Living Matter, Institute of Biological Information Processing and Institute
for Advanced Simulation, Forschungszentrum Jülich, 52425 Jülich, Germany
Lucio Isa  E-mail<EMAIL_ADDRESS>Laboratory for Soft Materials and
Interfaces, Department of Materials, ETH Zürich, Vladimir-Prelog-Weg 5, 8093
Zürich, Switzerland

(August 0 d , 2024)

###### Abstract

The generation of fluid flows by autophoretic microswimmers has been proposed
as a mechanism to enhance mass transport at the micro- and nanoscale. Here, we
experimentally investigate the ability of model “active baths” of
photocatalytic silica-titania Janus microspheres to enhance the diffusivity of
tracer particles at different microswimmer densities excluding collective
effects. Inspired by the similarities between our experimental findings and
previous results for biological microorganisms, we also model our Janus
microswimmers using a general squirmer framework. The numerical simulations
faithfully capture our observations, offer an insight into the microscopic
mechanism underpinning tracer transport, and allow us to expand the parameter
space beyond our experimental system. We find strong evidence for a collision-
based mechanism for increased tracer diffusivity, which leads to the
identification of an operating window for enhanced tracer transport by
chemical microswimmers based on dimensionless scaling arguments.

###### pacs:

???

The ability of active matter to transduce available energy into directed
motion has inspired research across various disciplines Ramaswamy (2010);
Gompper  _et al._ (2020). Amongst the systems studied at the micron-scale,
where Brownian fluctuations and viscous forces are significant Bechinger  _et
al._ (2016), the “autophoretic” motion of chemical microswimmers has been
proposed to induce micro-mixing in fluids Karshalev  _et al._ (2018); Dai  _et
al._ (2022), with potential applications ranging from environmental
remediation Wang  _et al._ (2019a) to synthetic chemistry Wittmann  _et al._
(2022). While biological microswimmers have been shown to enhance the
diffusivity of passive tracers Wu and Libchaber (2000); Leptos  _et al._
(2009); Valeriani  _et al._ (2011); Wilson  _et al._ (2011); Jeanneret  _et
al._ (2016); Zhao  _et al._ (2017); Lagarde  _et al._ (2020), analogous
systematic studies of this effect are lacking for synthetic, catalytic
microswimmers Miño  _et al._ (2011), except for systems propelled by bubble
generation Wang  _et al._ (2019b); Orozco  _et al._ (2014). As most synthetic
and biological active matter systems do not involve bubble propulsion, an
opportunity arises to address this open question with model active colloids.

Here, we investigate the mixing properties of a widely studied synthetic
active matter model system consisting of a suspension of photocatalytically
active Janus particles with a titania cap under UV illumination. We overcome
the limited scale of conventional Janus particle fabrication techniques via
‘Toposelective Nanoparticle Attachment” (TNA), which exploits an emulsion
template and multi-functional polymers to facilitate the selective binding of
a broad range of nanoparticles, including the photocatalytic titania
nanoparticles used in this study, onto silica particle supports Bailey  _et
al._ (2022a, 2023a). While experimentally probing the mechanisms leading to
the enhanced transport of passive tracers, we observe unexpected parallels
between their dynamics in active baths of our synthetic active colloids and
previous findings for biological microswimmers, despite their vastly different
underlying propulsion schemes. Employing a simple scaling argument, we
attribute the increased diffusivity of the colloidal tracers to a collision-
based interaction with the Janus microswimmers. We then use a coarse-grained
hydrodynamic swimmer model - which generates a force dipole and source dipole
in the far-field - to capture the generic aspects of tracer mixing Götze and
Gompper (2010); Zöttl and Stark (2023). Finding good agreement between
numerical simulations and experiments, we proceed to study the mixing
properties of our hydrodynamic “squirmers” under conditions beyond those
experimentally achievable, and highlight key considerations for the
realisation of the promise of catalytic ”micro-stirrers”.

Specifically, our microswimmers are spherical SiO2 colloids (diameter
$d_{S}=2.16\mu$m) half-decorated with photocatalytic TiO2 nanoparticles, as
described in detail in Refs. Bailey  _et al._ (2022a, 2023a). Under UV
illumination, the particles asymmetrically degrade H2O2 (fuel) and thus
establish local chemical gradients, leading to self-propulsion by autophoresis
Golestanian (2009); Popescu  _et al._ (2010); Dey  _et al._ (2016); Bailey
_et al._ (2022b). The density mismatch of the particles, coupled with the
experimental conditions used (see Supporting Information, Materials & Methods)
minimises their out-of-plane motion Bailey  _et al._ (2022c, 2023b) and
ensures that particles move in 2D (xy plane). Under these conditions, we
systematically study the behavior of the microswimmers as their effective area
fraction ($\Phi=N_{p}A_{p}/A_{T}$, where $N_{p}$ is the number of
microswimmers, $A_{p}=\pi d_{S}^{2}/4$, and $A_{T}$ is the area of the system)
is varied. In agreement with previous reports for both synthetic and
biological systems Theurkauff  _et al._ (2012); Thakur and Kapral (2012);
Palacci  _et al._ (2013a); Mognetti  _et al._ (2013); Liebchen  _et al._
(2017); Ginot  _et al._ (2018); Chen  _et al._ (2012, 2015); Theers  _et al._
(2018), we observe the formation of “dynamic clusters” above a threshold
swimmer density, (see Supporting Information, Figure S1, and compare Videos S1
and S2). The emergence of collective behaviours, although of significant
fundamental interest, not only complicates the analysis of the mixing
properties of microswimmers, but also likely hinders their overall mixing
efficiency Mou  _et al._ (2016), as we outline below. We therefore determine
the limiting area fraction of microswimmers before the onset of clustering
events using the methodology outlined by Theurkauff et al (see Supporting
Information, Figure S1a Theurkauff  _et al._ (2012)). In agreement with their
findings, we determine an onset of clustering at $\Phi\approx 0.05$ (see
Supporting Information, Figure S1c), which we thus set as the upper bound for
our experiments to focus on single-particle effects. Similar to the seminal
work by Wu and Libchaber on the mixing properties of bacterial baths Wu and
Libchaber (2000), the velocity distribution of our synthetic “free” particles
obeys Maxwell’s statistics, however, in its 2D form as particle motion is
confined above the substrate:
$P(v_{0})=\frac{v_{0}}{v_{M}^{2}}\text{exp}(-\frac{v_{0}^{2}}{2v_{M}^{2}})$,
where $v_{M}$ is the modal velocity. Under the experimental conditions used,
our active colloids have a modal velocity $v_{M}\approx 3.1\mu$ms-1 (see
Figure 1b).

Figure 1: Janus microswimmer 2D motion above a substrate in the absence of
clustering. a. Micrograph of “free” particles at low area fraction $\Phi$ =
0.02 (1s trajectory). Inset: HR-SEM image of the studied Janus particles.
Scale bar represents 1 $\mu$m. b. Distribution of the free microswimmer in-
plane velocities measured over 40000 trajectories.

Following previous work on biological microswimmers, we study the effect of
microswimmers on the dynamics of micron-sized tracers (diameter
$d_{T}=1.70\mu$m) Wu and Libchaber (2000); Leptos  _et al._ (2009); Miño  _et
al._ (2011); Valeriani  _et al._ (2011); Jeanneret  _et al._ (2016); Zhao  _et
al._ (2017); Lagarde  _et al._ (2020). We find that the tracers display
extensive displacements when a microswimmer passes in their vicinity (Figure
2a, Supporting Information Video S3) - which we henceforth refer to as
“jumps”. We stress that these enhancements in the motion of the tracers are
distinct to previous studies on “colloidal transport”, where the microswimmer
actively transports the tracer by a “pick-up and place” type mechanism Baraban
_et al._ (2012); Palacci  _et al._ (2013b); Mou  _et al._ (2019) (also see
Supporting Information Figure S2 on the over-weighted effect such transported
tracers can have on the ensemble mean-squared-displacement (EMSD)). By simple
scaling arguments, we now demonstrate that these “jumps” are indeed likely
caused by collision events between the tracers and the microswimmers Lagarde
_et al._ (2020).

Figure 2: Emergence of large displacements in the trajectories of tracers with
the introduction of microswimmers. a. Tracer trajectories for 3 different
experimental cases corresponding to various $\hat{\tau_{c}}$. Scale bars
represent 5 $\mu$m. Insets are representative snapshots of the experimental
system (bright tracers and dark microswimmers, 16.25 x 16.25 $\mu$m). b.
Evolution of the kurtosis $\gamma$ of the distribution of displacements with
lag time $\Delta\tau$ for the 3 experimental conditions depicted in a. Each
data-set contains at least 23 trajectories. The dashed lines indicate
$\hat{\tau_{c}}$ for the given experimental conditions. Inset: Correlation
between the estimated collision time $\hat{\tau_{c}}$ and the experimentally
observed lag time at which the kurtosis is maximum
($\tau_{\gamma=\operatorname*{argmax}_{\tau}(\gamma)}$ \- see Supporting
Information Figure S3). c. Ensemble mean-squared-displacement (EMSD) of the 3
experimental systems depicted in a.

We start by assuming that the tracers and swimmers are equidistantly
distributed. If the translational and rotational diffusion of the swimmers are
negligible, the MSD of N swimmers moving with speed $v_{0}$ at time $\tau$ can
be approximated as $\langle\Delta\bm{\mathrm{r}}^{\,2}(\tau)\rangle\sim
N_{p}v_{0}^{2}\tau^{2}$ at short times Bailey  _et al._ (2022d). Thus, the
time taken for a swimmer to (on average) travel the inter-tracer distance and
collide with the tracers can be approximated as
$\hat{\tau_{c}}\sim\frac{\sqrt{A_{T}}}{v_{0}\sqrt{N_{p}N_{trac}}}$, where
$N_{trac}$ is the number of tracer particles (see Supporting Information,
Materials and Methods). By plotting the kurtosis $\gamma$ of tracer
displacements for different lag times $\Delta\tau$ and fitting a smoothing
curve Fernandez-Rodriguez  _et al._ (2020), we observe a Pearson’s correlation
coefficient R = 0.8 between the scaling time $\hat{\tau_{c}}$ for collisions
estimated by our argumentation, and the experimentally determined lag time
$\tau_{\gamma=\operatorname*{argmax}_{\tau}(\gamma)}$ at which the maximum in
the kurtosis of the tracer displacements $\gamma(\Delta\tau)$ occurs (see
Figure 2b, inset).

$\gamma(\Delta\tau)$ is a measure of non-Gaussianity ($\gamma>0$ corresponding
to heavy tails, and $\gamma=0$ to Gaussian distributions), and can be expected
to capture the presence of tracer jumps. In previous theoretical and
experimental work on biological systems Pushkin and Yeomans (2013); Jeanneret
_et al._ (2016), the presence of extended tracer displacements upon near-field
interactions with the microorganisms was referred to as “entrainment”, arising
from no-slip surface interactions and enhancing tracer diffusivity Pushkin and
Yeomans (2013); Pushkin  _et al._ (2013); Mathijssen  _et al._ (2018). In our
synthetic analogues, tracers in active baths with low $\hat{\tau_{c}}$ also
exhibit a larger MSD and thus an increased diffusion coefficient (see Figure
2c). However, we expect that our Janus particles have a diffusiophoretic
surface slip velocity Golestanian  _et al._ (2007), precluding the mechanism
of entrainment Mathijssen  _et al._ (2018).

Figure 3: Influence of the jumps on the tracer distribution of displacements
at different lag times $\Delta\tau$. Top row: Tracer trajectory in a
microswimmer bath. The non-Gaussian, exponential tails in the distribution of
displacements becomes prominent at longer lag-times. We note that these tails
are not observed for “passive” microswimmers (see Supporting Information,
Figure S4). Bottom row: Our numerical neutral squirmer model maps
qualitatively well onto the experimental findings, with the presence of jumps
at longer $\Delta\tau$ seen in the tails of the distribution of displacements
(blue for tracers in the presence of neutral squirmers, black for the passive
case).

Further evidence for the enhancement in tracer diffusivity as a result of
tracer jumps is presented in Figure 3 (top row). Following Jeanneret et al.,
we distinguish between standard diffusive motion during a fixed time and
jumps, using the autocorrelation of subsequent displacements Jeanneret  _et
al._ (2016). Although the significantly lower speeds of the studied synthetic
microswimmers reduces the effective size of the tracer jumps, it is clear that
long-time exponential tails emerge in the distribution of tracer displacements
when jumps are included, replicating the findings for microalgae Jeanneret
_et al._ (2016). This qualitative agreement was not expected a-priori, but
further motivates our modelling of the active colloids as squirmers.

The squirmer model Götze and Gompper (2010), first proposed for microorganisms
Lighthill (1952); Blake (1971), is frequently invoked to describe the flow
fields surrounding chemical Janus microswimmers Michelin and Lauga (2014);
Zöttl and Stark (2023); Bailey  _et al._ (2023b), and generates flows which
are consistent with a coarse-grained framework incorporating hydrodynamic
effects in the far field via a multipole expansion, together with appropriate
near-field fluid flows. We note that a more complete description of chemically
active colloids requires the analysis of the phoretic flows resulting from the
chemical gradients produced by the microswimmers Popescu  _et al._ (2018);
Katuri  _et al._ (2021) (see Supporting Information, Figure S5). Still, a
complete description of the phoretic and osmotic flows requires the relevant
mobility parameters - reducing generalisability to other mixing systems. For a
coarse-grained model to describe the mixing properties of our microswimmers,
we therefore neglect the phoretic contribution to tracer displacements,
following Ref. Campbell  _et al._ (2019) (also see Supporting Information,
Figure S6). In particular, we model our synthetic microswimmers as neutral
squirmers (squirmer parameter $\beta$ = 0 - see Supporting Information,
Materials and Methods), informed by previous work on catalytic Janus systems
Yang  _et al._ (2014). The neutral squirmer model furthermore agrees with the
absence of transport behaviours associated with flow vortices Madden  _et al._
(2022), the dominant contribution of tracer “jumps” to their diffusivity, and
the absence of long-range flow effects Pushkin and Yeomans (2013). Fluid flow
is modelled using dissipative particle dynamics (DPD), a mesoscopic
hydrodynamics simulation technique Hoogerbrugge and Koelman (1992); Espanol
and Warren (1995). Both passive and active colloids are represented by a
triangulated network of particles placed homogeneously on a spherical surface
to ensure a near-rigid shape (see Supporting Information, Materials and
Methods). Propulsion is modelled by imposing a slip velocity at the surface of
the active particles. For passive particles, no-slip boundary conditions (BCs)
are imposed. The simulation domain corresponds to 34 x 34 x 6.8 $\mu$m with
periodic BCs in the xy plane, and walls at z = 0 and z = 6.8 $\mu$m. The
diameter $d_{Sq}$ of the squirmers is taken to be 2.16 $\mu$m, as in
experiments. Both passive and active particles are confined within a layer of
2.45 $\mu$m at the lower wall with no-slip BCs, to mimic the quasi 2D setup in
the experiments. The wall at z = 6.8 $\mu$m assumes slip BCs, to better
capture experimental conditions. The area fraction of active particles
corresponds to $\Phi$ = 0.05, identified as the onset of clustering in both
experiments and simulations. We underline that our simulation of squirmers
with a surface velocity precludes tracer entrainment, which requires a no-slip
microswimmer surface Mathijssen  _et al._ (2018).

We find that the dynamics of our simulated “tracer” particles in the presence
of neutral squirmers are very similar to our experimental findings at similar
velocities. Specifically, we identify the presence of “jumps” in the tracer
trajectories, which emerge at longer lag times in the distribution of tracer
displacements (see Figure 3, bottom row). We therefore proceed to a more
quantitative characterisation of the mixing efficiencies of different
squirmers, not only for neutral squirmers ($\beta$ = 0), but also for pullers
($\beta$ = 5), and pushers ($\beta$ = -5). As a control, we also study the
dynamics of tracer particles in the presence of particles that do not swim
(i.e. passive particles). We observe that pushers are significantly superior
mixers to other squirmer types in our experimental conditions, as reflected in
the EMSD of the tracers (Figure 4a). In fact, we determine a linear trend in
Deff/D0 with $\langle$V${}_{0}\rangle$ (Figure 4b), as expected for a
predominantly collision-based mechanism for mixing. The superior mixing
performance of pushers is attributed to their higher in-plane swimming speeds,
which results from the preferential alignment of their swimming axis parallel
to the substrate, whereas pullers are often oriented perpendicular to the
surface (see Figure 4c and inset). As an aside, we note the similarity in the
distribution of swimming speeds of our synthetic microswimmers and the neutral
squirmers (compare Figures 1b and 4c), underlining the role of orientation
with respect to the substrate.

Figure 4: Dynamics of simulated active and tracer particles (black for passive
particles, red for pullers, blue for neutrals, and magenta for pushers.) a.
Mean-squared displacements of tracers for different squirmer types for
velocities approximating experimental values. b. Scaling of normalised tracer
diffusion Deff/D0 with the swimming velocity V0 of the squirmers. c.
Distributions of squirmer in-plane (2D) velocities for fixed squirmer
strengths (all types would have the same swimming velocity under unbounded
conditions). Inset: Orientation distributions of active and tracer particles
with respect to the z axis (i.e. perpendicular to the wall). d. Deff/D0 for
different sizes of passive particles.

The tracers’ diffusivity enhancement Deff/D0 also scales approximately
linearly with their size in the resolution of our DPD simulations (for a given
squirmer velocity - see Figure 4d), indicating the robustness of the
“collision” mechanism to not only particle speed but also to the size of the
tracers. We thus hypothesise that the mixing capacity of microswimmers is
limited by the relative size and frequency of the imparted “jumps” via
collisions compared to the underlying diffusivity of the component to be mixed
(i.e. $l_{j}\gg\sqrt{\text{D}_{0}\tau_{c}}$ for enhanced mixing, where $l_{j}$
is the jump length after a collision occurring with frequency $1/\tau_{c}$).
This mirrors the experimental findings using 230 nm nanoparticle tracers where
no enhancement in $D_{eff}$ over $D_{0}$ was observed (see Supporting
Information Figures S7-8), suggesting the unsuitability of this type of
microswimmer systems to mix e.g. species with very high diffusion
coefficients.

Despite the best laid schemes for microswimmers and mixing, there are still
significant hurdles to the “micro-stirrer” concept. This can perhaps be simply
illustrated by taking a leaf from Purcell Purcell (1977), and evaluating the
relevant dimensionless numbers involved in enhanced transport for chemical
microswimmers. Presuming low Reynolds number conditions $Re=vl/\nu\ll 1$,
where $\nu$ is the kinematic viscosity and $vl$ relates to the fluid flow at a
characteristic length-scale $l$, for microswimmers swimming in bulk with
velocity $v$ (i.e. motion in the viscous regime and ignoring substrate
effects), the first relevant ratio is the Péclet “transport” number,
$Pe_{T}=vl/D$ \- where $D$ is the translational diffusion coefficient of the
transported species. If $Pe_{T}\ll 1$, the microswimmer is essentially static
with respect to the transported species, and we do not expect any diffusivity
enhancement. Considering the special case of “chemistry-on-the-fly”, where
molecules are “mixed” to enhance the chemical reaction rate Karshalev  _et
al._ (2018); Dai  _et al._ (2022); Wittmann  _et al._ (2022),
$Pe_{T}\rightarrow 0$ for most experimental systems. However, $Pe_{T}$ can be
increased by scaling the particle size and speed appropriately (e.g.
$d_{S}\sim 10^{2}$ $\mu$m moving at some $d_{S}$ per second). In this case,
and staying on the “chemistry-on-the-fly” approach (while ignoring phoretic
effects Michelin and Lauga (2014)), we can then evaluate the Damköhler $Da$
and Sherwood $Sh$ numbers. $Da=\kappa l/D$ details the balance between the
reaction rate $\kappa$ and the (diffusive) mass transfer process. Notably, the
classical continuum phoretic framework to describe chemical microswimmer
motion presumes reaction-rate limited kinetics, i.e. $Da\rightarrow 0$ Howse
_et al._ (2007); Kreissl  _et al._ (2016); Popescu  _et al._ (2018); Domínguez
and Popescu (2022), with finite values of $Da$ leading to decreasing
microswimmer speeds Michelin and Lauga (2014) \- calling into question the
goal of overcoming mass-transfer limitations in reactions using microswimmer
motility Karshalev  _et al._ (2018); Dai  _et al._ (2022); Wittmann  _et al._
(2022). Nevertheless, assuming transport-limited reactions, we examine the
Sherwood number for mass transfer to a sphere in Stokes flow,
$Sh=2+0.991Re^{1/3}Sc^{1/3}$ Friedlander (1961); Armenante and Kirwan (1989),
where $Sc=\nu/D$ is the Schmidt number (the ratio of momentum to mass
diffusivity), the first term is for the asymptotic limit of molecular
diffusivity, and $0.991Re^{1/3}Sc^{1/3}$ is the advective contribution to
mass-transport Armenante and Kirwan (1989). In creeping flow, we thus expect
the enhancement in mass transfer from fluid advection to be small even when
$Pe_{T}>1$ for most physical parameters. In highly viscous reactor settings,
the micro-stirrer strategy could be effective for e.g. macromolecules with a
low diffusivity (assuming the microswimmers do not slow down), however, we
underline the poor scaling with the exponent $1/3$.

Basic dimensionless arguments based on convection-reaction-diffusion theory
Squires  _et al._ (2008) thus provide some additional insight into the
challenges of using autophoretic microswimmers as micro-stirrers, besides the
limitations posed by neglected effects such as clustering. Still, it is worth
noting that collective behaviours could also enhance fluid intermixing and
thereby improve reaction-diffusion limitations Wang  _et al._ (2021), while
appropriate material selection can modify generated flow fields Katuri  _et
al._ (2018); Sharan  _et al._ (2022) and thereby prevent undesirable
aggregation phenomena (i.e. for pushers Theers  _et al._ (2018); Qi  _et al._
(2022)). Other promising avenues to increase microswimmer mixing could be
achieved via modifying their shape away from classical spheres Xiong  _et al._
(2023), or exploiting bubble generation Mundaca-Uribe  _et al._ (2023). Our
results therefore underline the need for a more considered approach to the
design of synthetic microswimmers to enhance micro- and nanoscale mass
transport, and provide evidence that deviations from model (spherical)
swimmers may be required for such applications.

Acknowledgements

The authors thank M. Popescu and V. Niggel for helpful discussions on theory
and experiment respectively. The authors also thank S. Vasudevan for his
assistance with the performed DDM analysis. The authors gratefully acknowledge
computing time on the supercomputer JURECA Thörnig (2021) at Forschungszentrum
Jülich under grant no. actsys.

Author Contribution Statement

Author contributions are defined based on the CRediT (Contributor Roles
Taxonomy). Conceptualisation: M.R.B., D.A.F., F.P., F.G., G.G., L.I. Funding
acquisition: L.I. Investigation and analysis: M.R.B., D.A.F. Methodology:
M.R.B., D.A.F., F.P. Supervision: L.I. Visualisation and writing: M.R.B.,
D.A.F., F.P., F.G., G.G., L.I.

## References

  * Ramaswamy (2010) S. Ramaswamy, Annual Review of Condensed Matter Physics 1, 323 (2010).
  * Gompper  _et al._ (2020) G. Gompper, R. G. Winkler, T. Speck, A. Solon, C. Nardini, F. Peruani, H. Löwen, R. Golestanian, U. B. Kaupp, _et al._ , Journal of Physics: Condensed Matter 32, 193001 (2020).
  * Bechinger  _et al._ (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Reviews of Modern Physics 88, 045006 (2016).
  * Karshalev  _et al._ (2018) E. Karshalev, B. Esteban-Fernández De Ávila, and J. Wang, Journal of the American Chemical Society 140, 3810 (2018).
  * Dai  _et al._ (2022) B. Dai, Y. Zhou, X. Xiao, Y. Chen, J. Guo, C. Gao, Y. Xie, and J. Chen, Advanced Science 9 (2022).
  * Wang  _et al._ (2019a) L. Wang, A. Kaeppler, D. Fischer, and J. Simmchen, ACS Applied Materials and Interfaces 11, 32937 (2019a).
  * Wittmann  _et al._ (2022) M. Wittmann, S. Heckel, F. Wurl, Z. Xiao, T. Gemming, T. Strassner, and J. Simmchen, Chemical Communications 58, 4052 (2022).
  * Wu and Libchaber (2000) X. L. Wu and A. Libchaber, Physical Review Letters 84, 3017 (2000).
  * Leptos  _et al._ (2009) K. C. Leptos, J. S. Guasto, J. P. Gollub, A. I. Pesci, and R. E. Goldstein, Physical Review Letters 103, 1 (2009).
  * Valeriani  _et al._ (2011) C. Valeriani, M. Li, J. Novosel, J. Arlt, and D. Marenduzzo, Soft Matter 7, 5228 (2011).
  * Wilson  _et al._ (2011) L. G. Wilson, V. A. Martinez, J. Schwarz-Linek, J. Tailleur, G. Bryant, P. N. Pusey, and W. C. Poon, Physical Review Letters 106, 7 (2011).
  * Jeanneret  _et al._ (2016) R. Jeanneret, D. O. Pushkin, V. Kantsler, and M. Polin, Nature Communications 7, 1 (2016).
  * Zhao  _et al._ (2017) X. Zhao, K. K. Dey, S. Jeganathan, P. J. Butler, U. M. Córdova-Figueroa, and A. Sen, Nano Letters 17, 4807 (2017).
  * Lagarde  _et al._ (2020) A. Lagarde, N. Dagès, T. Nemoto, V. Démery, D. Bartolo, and T. Gibaud, Soft Matter 16, 7503 (2020).
  * Miño  _et al._ (2011) G. Miño, T. E. Mallouk, T. Darnige, M. Hoyos, J. Dauchet, J. Dunstan, R. Soto, Y. Wang, A. Rousselet, and E. Clement, Physical Review Letters 106, 48102 (2011).
  * Wang  _et al._ (2019b) B. Wang, F. Ji, J. Yu, L. Yang, Q. Wang, and L. Zhang, iScience 19, 760 (2019b).
  * Orozco  _et al._ (2014) J. Orozco, B. Jurado-Sánchez, G. Wagner, W. Gao, R. Vazquez-Duhalt, S. Sattayasamitsathit, M. Galarnyk, A. Cortés, D. Saintillan, and J. Wang, Langmuir 30, 5082 (2014).
  * Bailey  _et al._ (2022a) M. R. Bailey, F. Grillo, N. D. Spencer, and L. Isa, Advanced Functional Materials 32, 2109175 (2022a).
  * Bailey  _et al._ (2023a) M. R. Bailey, T. A. Gmür, F. Grillo, and L. Isa, Chemistry of Materials 35, 3731–3741 (2023a).
  * Götze and Gompper (2010) I. O. Götze and G. Gompper, Physical Review E 82, 041921 (2010).
  * Zöttl and Stark (2023) A. Zöttl and H. Stark, Annual Review of Condensed Matter Physics 14, 109 (2023).
  * Golestanian (2009) R. Golestanian, Physical Review Letters 102, 188305 (2009).
  * Popescu  _et al._ (2010) M. N. Popescu, S. Dietrich, M. Tasinkevych, and J. Ralston, European Physical Journal E 31, 351 (2010).
  * Dey  _et al._ (2016) K. K. Dey, F. Wong, A. Altemose, and A. Sen, Current Opinion in Colloid and Interface Science 21, 4 (2016).
  * Bailey  _et al._ (2022b) M. R. Bailey, N. Reichholf, A. Flechsig, F. Grillo, and L. Isa, Particle and Particle Systems Characterization 39, 2100200 (2022b).
  * Bailey  _et al._ (2022c) M. R. Bailey, F. Grillo, and L. Isa, Soft Matter 18, 7291–7300 (2022c).
  * Bailey  _et al._ (2023b) M. R. Bailey, C. M. B. Gutiérrez, J. Martín-Roca, V. Niggel, V. Carrasco-Fadanelli, I. Buttinoni, I. Pagonabarraga, L. Isa, and C. Valeriani, arXiv (2023b), 10.48550/arXiv.2307.16193.
  * Theurkauff  _et al._ (2012) I. Theurkauff, C. Cottin-Bizonne, J. Palacci, C. Ybert, and L. Bocquet, Physical Review Letters 108, 268303 (2012).
  * Thakur and Kapral (2012) S. Thakur and R. Kapral, Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 85, 026121 (2012).
  * Palacci  _et al._ (2013a) J. Palacci, S. Sacanna, A. P. Steinberg, D. J. Pine, and P. M. Chaikin, Science 339, 936 (2013a).
  * Mognetti  _et al._ (2013) B. M. Mognetti, A. Šarić, S. Angioletti-Uberti, A. Cacciuto, C. Valeriani, and D. Frenkel, Physical Review Letters 111, 245702 (2013).
  * Liebchen  _et al._ (2017) B. Liebchen, D. Marenduzzo, and M. E. Cates, Physical Review Letters 118, 268001 (2017).
  * Ginot  _et al._ (2018) F. Ginot, I. Theurkauff, F. Detcheverry, C. Ybert, and C. Cottin-Bizonne, Nature Communications 9, 1 (2018).
  * Chen  _et al._ (2012) X. Chen, X. Dong, A. Be’er, H. L. Swinney, and H. P. Zhang, Physical Review Letters 108, 148101 (2012).
  * Chen  _et al._ (2015) X. Chen, X. Yang, M. Yang, and H. P. Zhang, EPL (Europhysics Letters) 111, 54002 (2015).
  * Theers  _et al._ (2018) M. Theers, E. Westphal, K. Qi, R. G. Winkler, and G. Gompper, Soft matter 14, 8590 (2018).
  * Mou  _et al._ (2016) F. Mou, L. Kong, C. Chen, Z. Chen, L. Xu, and J. Guan, Nanoscale 8, 4976 (2016).
  * Baraban  _et al._ (2012) L. Baraban, M. Tasinkevych, M. N. Popescu, S. Sanchez, S. Dietrich, and O. G. Schmidt, Soft Matter 8, 48 (2012).
  * Palacci  _et al._ (2013b) J. Palacci, S. Sacanna, A. Vatchinsky, P. M. Chaikin, and D. J. Pine, Journal of the American Chemical Society 135, 15978 (2013b).
  * Mou  _et al._ (2019) F. Mou, J. Zhang, Z. Wu, S. Du, Z. Zhang, L. Xu, and J. Guan, iScience 19, 415 (2019).
  * Bailey  _et al._ (2022d) M. R. Bailey, A. R. Sprenger, F. Grillo, H. Löwen, and L. Isa, Physical Review E 106, L052602 (2022d).
  * Fernandez-Rodriguez  _et al._ (2020) M. A. Fernandez-Rodriguez, F. Grillo, L. Alvarez, M. Rathlef, I. Buttinoni, G. Volpe, and L. Isa, Nature Communications 11, 1 (2020).
  * Pushkin and Yeomans (2013) D. O. Pushkin and J. M. Yeomans, Physical Review Letters 111, 1 (2013).
  * Pushkin  _et al._ (2013) D. O. Pushkin, H. Shum, and J. M. Yeomans, Journal of Fluid Mechanics 726, 5 (2013).
  * Mathijssen  _et al._ (2018) A. J. Mathijssen, R. Jeanneret, and M. Polin, Physical Review Fluids 3, 033103 (2018).
  * Golestanian  _et al._ (2007) R. Golestanian, T. B. Liverpool, and A. Ajdari, New Journal of Physics 9 (2007), 10.1088/1367-2630/9/5/126, arXiv:0701168 [cond-mat] .
  * Lighthill (1952) M. J. Lighthill, Communications on Pure and Applied Mathematics 5, 109 (1952).
  * Blake (1971) J. R. Blake, Mathematical Proceedings of the Cambridge Philosophical Society 70, 303 (1971).
  * Michelin and Lauga (2014) S. Michelin and E. Lauga, Journal of fluid mechanics 747, 572 (2014).
  * Popescu  _et al._ (2018) M. N. Popescu, W. E. Uspal, Z. Eskandari, M. Tasinkevych, and S. Dietrich, European Physical Journal E 41, 1 (2018).
  * Katuri  _et al._ (2021) J. Katuri, W. E. Uspal, M. N. Popescu, and S. Sánchez, Science Advances 7, eabd0719 (2021).
  * Campbell  _et al._ (2019) A. I. Campbell, S. J. Ebbens, P. Illien, and R. Golestanian, Nature Communications 10, 1 (2019).
  * Yang  _et al._ (2014) M. Yang, A. Wysocki, and M. Ripoll, Soft Matter 10, 6208 (2014).
  * Madden  _et al._ (2022) I. P. Madden, L. Wang, J. Simmchen, and E. Luijten, Small 18, 2107023 (2022).
  * Hoogerbrugge and Koelman (1992) P. J. Hoogerbrugge and J. M. Koelman, Epl 19, 155 (1992).
  * Espanol and Warren (1995) P. Espanol and P. Warren, Epl 30, 191 (1995).
  * Purcell (1977) E. M. Purcell, American Journal of Physics 45, 3 (1977).
  * Howse  _et al._ (2007) J. R. Howse, R. A. Jones, A. J. Ryan, T. Gough, R. Vafabakhsh, and R. Golestanian, Physical Review Letters 99, 048102 (2007).
  * Kreissl  _et al._ (2016) P. Kreissl, C. Holm, and J. De Graaf, The Journal of Chemical Physics 144 (2016).
  * Domínguez and Popescu (2022) A. Domínguez and M. N. Popescu, Current Opinion in Colloid and Interface Science 61, 101610 (2022).
  * Friedlander (1961) S. Friedlander, AIChE Journal 7, 347 (1961).
  * Armenante and Kirwan (1989) P. M. Armenante and D. J. Kirwan, Chemical Engineering Science 44, 2781 (1989).
  * Squires  _et al._ (2008) T. M. Squires, R. J. Messinger, and S. R. Manalis, Nature biotechnology 26, 417 (2008).
  * Wang  _et al._ (2021) Q. Wang, F. Ji, S. Wang, and L. Zhang, ChemNanoMat 7, 600 (2021).
  * Katuri  _et al._ (2018) J. Katuri, W. E. Uspal, J. Simmchen, A. Miguel-López, and S. Sánchez, Science Advances 4, eaao1755 (2018).
  * Sharan  _et al._ (2022) P. Sharan, Z. Xiao, V. Mancuso, W. E. Uspal, and J. Simmchen, ACS Nano , 4599 (2022).
  * Qi  _et al._ (2022) K. Qi, E. Westphal, G. Gompper, and R. G. Winkler, Communications Physics 5, 49 (2022).
  * Xiong  _et al._ (2023) K. Xiong, J. Lin, Q. Chen, T. Gao, L. Xu, and J. Guan, Matter 6, 1–18 (2023).
  * Mundaca-Uribe  _et al._ (2023) R. Mundaca-Uribe, M. Holay, A. Abbas, N. Askarinam, J. S. Sage-Sepulveda, L. Kubiatowicz, R. H. Fang, L. Zhang, and J. Wang, ACS Nano 17, 9272 (2023).
  * Thörnig (2021) P. Thörnig, Journal of large-scale research facilities JLSRF 7, A182 (2021).
  * Fedosov  _et al._ (2010) D. A. Fedosov, B. Caswell, and G. E. Karniadakis, Biophys. J. 98, 2215 (2010).
  * Helfrich (1973) W. Helfrich, Z. Naturforsch. 28, 693 (1973).
  * Alarcón  _et al._ (2017) F. Alarcón, C. Valeriani, and I. Pagonabarraga, Soft Matter 13, 814 (2017).
  * Singh  _et al._ (2017) D. P. Singh, U. Choudhury, P. Fischer, and A. G. Mark, Advanced Materials 29, 1701328 (2017).
  * Kepten  _et al._ (2013) E. Kepten, I. Bronshtein, and Y. Garini, Physical Review E - Statistical, Nonlinear, and Soft Matter Physics 87, 1 (2013), arXiv:1212.0793 .
  * Bailey  _et al._ (2022e) M. R. Bailey, A. R. Sprenger, F. Grillo, H. Löwen, and L. Isa, Physical Review E 106 (2022e), 10.1103/physreve.106.l052602.
  * Cerbino and Trappe (2008) R. Cerbino and V. Trappe, Physical Review Letters 100 (2008), 10.1103/physrevlett.100.188102.

## The Supporting Information includes:

  * Text S1: Materials and Methods

  * Fig. S1: Definition of clustering and the different velocity distributions

  * Fig. S2: Influence of actively transported tracers on ensemble averaged values

  * Fig. S3: Fitting the lag-time at which the maximal kurtosis for the distribution of displacements of tracers is observed

  * Fig. S4: Distribution of tracer displacements in the presence of non-swimming microswimmers (UV light off)

  * Figure S5: Evidence for phoretic interactions between the Janus particles (pumps) and passive nanoparticle tracers

  * Figure S6: Influence of Janus pumps on nanoparticle tracer dynamics in the near- and far-field

  * Figure S7: Lack of mixing effect observed for nanoparticle tracers in a microswimmer bath

  * Figure S8: Comparison of different approaches to evaluate nanoparticle tracer dynamics to evaluate their accuracies

## Materials and Methods

Experimental setup and materials. Particles were imaged on an inverted
microscope (Nikon Eclipse Ti2e) under Köhler illumination with white light
using a 40× objective (CFI S Plan Fluor ELWD 40XC) with adjustable collar (set
to 1 mm) unless otherwise specified. Due to a density mismatch, the particles
sediment to the glass substrate, where they remain under the experimental
conditions used. Images were acquired at 50 FPS using a Hamamatsu C14440-20UP
digital camera. The Janus microswimmers were activated with a UV LED (365 nm,
Prizmatix UHP-F-365), at a power density of 44 mW/mm2, in the presence of a
chemical fuel (H2O2 3 % v/v, Acros Organics). The tracer particles (230 nm and
1700 nm diameter, microParticle GmbH) were fluorescent under illumination with
the UV LED source, and did not require an additional excitation source (Ex/Em
502/518 nm). To ensure a near-constant number density (area fraction $\Phi$)
of studied microswimmers, it was necessary to create a confined system to
prevent their escape to non-illuminated regimes. This was achieved by creating
wells (diameter = 250 $\mu$m) in glass slides by HF etching. To avoid boundary
effects, only the region bounded by a box (83.2 x 83.2 $\mu$m) centred in this
well was studied. The cells were sealed using an upper glass slide and PDMS
spacer to prevent convective drift.

The expression for $\hat{\tau_{c}}$ was derived as follows: Assuming that N
ballistic swimmers (ignoring diffusion) will explore an area given by:

$\displaystyle\langle\Delta\bm{\mathrm{r}}^{\,2}(\tau)\rangle\sim
N_{p}v_{0}^{2}\tau^{2}.$ (1)

then it follows that they will travel a distance $l_{N}$ in
$\tau=l_{N}/(v_{0}\sqrt{N_{p}})$ seconds.

Likewise, assuming $N_{t}$ equidistantly spaced tracers in the system area
$A$, we approximate the mean inter-tracer distance to be $\langle
l_{t}\rangle=\sqrt{A}/\sqrt{N_{t}}$. Therefore, the average predicted time
$\hat{\tau_{c}}$ for N ballistic swimmers to travel the mean inter-tracer
distance can be obtained by equating $l_{N}$ with $\langle l_{t}\rangle$,
giving the following collision time between swimmers and tracers:

$\displaystyle\hat{\tau_{c}}\sim\frac{\sqrt{A_{T}}}{v_{0}\sqrt{N_{p}N_{trac}}}$
(2)

Modeling fluid flow in simulations. Fluid flow is modeled by the dissipative
particle dynamics (DPD) method Hoogerbrugge and Koelman (1992); Espanol and
Warren (1995), which is a mesoscopic hydrodynamics simulation technique. In
DPD, fluid is represented by a collection of $N_{s}$ solvent particles. DPD
particle interactions are pairwise and soft, while particle velocities ${\bf
v}_{i}$ and positions ${\bf r}_{i}$ follow the Newton’s second law of motion

$\displaystyle\frac{d{\bf r}_{i}}{dt}$ $\displaystyle=$ $\displaystyle{\bf
v}_{i},$ $\displaystyle\frac{d{\bf v}_{i}}{dt}$ $\displaystyle=$
$\displaystyle\frac{1}{m_{i}}\sum_{j\neq i}({\bf F}^{C}_{ij}+{\bf
F}^{D}_{ij}+{\bf F}^{R}_{ij}),$ (3)

where $m_{i}$ is the particle mass and the forces ${\bf F}^{C}_{ij}$, ${\bf
F}^{D}_{ij}$, and ${\bf F}^{R}_{ij}$ define conservative, dissipative, and
random interactions, respectively. The above equations of motion are
integrated using the velocity-Verlet algorithm. The total force exerted on a
particle $i$ by particle $j$ is given by

$\displaystyle{\bf F}^{C}_{ij}$ $\displaystyle=$ $\displaystyle
b\left(1-\frac{r_{ij}}{r_{c}}\right)\hat{\bf r}_{ij},$ $\displaystyle{\bf
F}^{D}_{ij}$ $\displaystyle=$ $\displaystyle-\gamma\omega^{D}(r_{ij})({\bf
v}_{ij}\cdot\hat{\bf r}_{ij})\hat{\bf r}_{ij},$ $\displaystyle{\bf
F}^{R}_{ij}$ $\displaystyle=$
$\displaystyle\sigma_{dpd}\omega^{R}(r_{ij})\frac{\xi_{ij}}{\sqrt{dt}}\hat{\bf
r}_{ij},$ (4)

where ${\bf r}_{ij}={\bf r}_{i}-{\bf r}_{j}$, $\hat{\bf r}_{ij}={\bf
r}_{ij}/r_{ij}$, and ${\bf v}_{ij}={\bf v}_{i}-{\bf v}_{j}$. The coefficients
$b$, $\gamma$, and $\sigma_{dpd}$ determine the strength of conservative,
dissipative, and random forces, respectively. $\omega^{D}$ and $\omega^{R}$
are weight functions and $\xi_{ij}$ is a symmetric normally-distributed random
variable with zero mean and unit variance. All forces act within a sphere of
cutoff radius $r_{c}$. The random and dissipative forces form a thermostat and
have to satisfy the fluctuation-dissipation theorem Espanol and Warren (1995)
to maintain the equilibrium temperature $T$ in DPD system. This results in two
conditions,

$\omega^{D}(r_{ij})=\left[\omega^{R}(r_{ij})\right]^{2},\quad\quad\sigma_{dpd}^{2}=2\gamma
k_{B}T.$ (5)

The typical choice for the weight function is
$\omega^{R}(r_{ij})=(1-r_{ij}/r_{c})^{s}$, where $s$ is a selected exponent.

Rigid walls are simulated by immobilized DPD particles having identical number
density and interactions as the fluid particles. The immobile wall particles
enforce no-slip boundary conditions (BCs) at the walls through dissipative
interactions with the fluid particles. In order to prevent wall penetration by
fluid particles, reflective surfaces are introduced at the interface between
the fluid and the walls. Slip BCs at the walls can be modeled by turning off
the dissipative interactions between fluid and wall particles.

Model of active particles and passive tracers in simulations. Both active
swimmers and passive tracers are represented by an arrangement of particles
uniformly distributed on a spherical surface and forming a triangulated
network (i.e. membrane model). To ensure the structural integrity of the
triangulated network, it is supplied with shear elasticity, bending rigidity,
and constraints related to the surface area and enclosed volume Fedosov  _et
al._ (2010). The bonds connecting adjacent vertices within the network provide
shear elasticity defined by a potential function $U_{e}=U_{wlc}+U_{pow}$ as
Fedosov  _et al._ (2010).

$U_{wlc}=\frac{k_{B}Tl_{m}}{4p}\frac{3x^{2}-2x^{3}}{1-x},\quad\quad
U_{pow}=\frac{k_{p}}{l},$ (6)

where $l$ is the spring length, $l_{m}$ is the maximum spring extension,
$x=l/l_{m}$, $p$ is the persistence length, and $k_{p}$ is the repulsive force
coefficient. The bending rigidity is described by the Helfrich bending energy
Helfrich (1973)

$U_{b}=\frac{\kappa}{2}\sum_{i}\sigma_{i}(H_{i}-H_{0}^{i})^{2},$ (7)

where $\kappa$ is the bending rigidity, $H_{i}={\bf
n}_{i}\cdot\sum_{j(i)}\sigma_{ij}{\bf r}_{ij}/(\sigma_{i}r_{ij})$ is the
discretized mean curvature at vertex $i$, ${\bf n}_{i}$ is the unit normal at
the membrane vertex $i$, $\sigma_{i}=\sum_{j(i)}\sigma_{ij}r_{ij}/4$ is the
area corresponding to vertex $i$ (the area of the dual cell), $j(i)$
corresponds to all vertices linked to vertex $i$, and
$\sigma_{ij}=r_{ij}(\cot\theta_{1}+\cot\theta_{2})/2$ is the length of the
bond in the dual lattice, where $\theta_{1}$ and $\theta_{2}$ are the angles
at the two vertices opposite to the edge $ij$ in the dihedral. Finally,
$H_{0}^{i}$ is the spontaneous curvature at vertex $i$.

The area and volume conservation constraints are given by Fedosov  _et al._
(2010)

$U_{a}=\frac{k_{a}(A-A_{0}^{tot})^{2}}{2A_{0}^{tot}}+\sum_{m\in
1...N_{t}}{\frac{k_{d}(A_{m}-A^{m}_{0})^{2}}{2A^{m}_{0}}},$ (8)
$U_{v}=\frac{k_{v}(V-V_{0}^{tot})^{2}}{2V_{0}^{tot}},$ (9)

where $A$ is the current area of the membrane, $A_{0}^{tot}$ is the targeted
global area, $A_{m}$ is the area of the m-th triangle from $N_{t}$ triangles
on the triangulated surface, $A^{m}_{0}$ is the targeted area of the m-th
triangle determined by the triangular mesh on the particle surface, $V$ is the
instantaneous membrane volume, and $V_{0}^{tot}$ is the targeted volume.
Furthermore, $k_{a}$, $k_{d}$, and $k_{v}$ are the global area, local area and
volume constraint coefficients, respectively.

Both swimmers and tracers are submerged within a DPD fluid, and also filled by
fluid DPD particles due to their membrane-like representation. The membrane
surfaces of all suspended active particles and tracers serve as a boundary
separating DPD particles inside and outside of the membranes. This is achieved
through the reflection of solvent particles at the membrane surfaces from
inside and outside. Note that the dissipative and random forces between the
internal and external fluid particles are deactivated, and only the
conservative force is employed to maintain uniform fluid pressure across the
membranes. No slip BCs between the fluid and suspended colloidal particles are
implemented through the proper adjustment of the friction coefficient $\gamma$
between fluid and membrane particles Fedosov  _et al._ (2010).

Propulsion of active swimmers is realized through a slip velocity at their
surface using the spherical squirmer model Alarcón  _et al._ (2017)

${\bf u}_{s}=\left(B_{1}\sin(\phi)+B_{2}\sin(\phi)\cos(\phi)\right){\bf
t_{\phi}},$ (10)

where ${\bf t_{\phi}}$ represents the local tangential vector in the $\phi$
direction at the swimmer surface, and $\phi$ corresponds to the angle between
the orientation of the squirmer and the vector extending from its center of
mass to a local position on its surface. The coefficient $B_{1}$ governs the
propulsion strength of the squirmer, while $B_{2}$ introduces a front-back
asymmetry into the slip-velocity field. In bulk, the spherical squirmer swims
with a velocity $v_{0}=2B_{1}/3$. We employ the dimensionless ratio
$\beta=\frac{B_{2}}{B_{1}}$ to characterize the mode of propulsion, where
$\beta<0$ represents a pusher, $\beta=0$ corresponds to a neutral swimmer, and
$\beta>0$ models a puller. To enforce the slip velocity at the squirmer
surface, the dissipative interaction between the squirmer particles and those
of the surrounding fluid is adjusted as follows

${\bf F}_{ij}^{D}=-\gamma\omega^{D}(r_{ij})(\hat{\bf r}_{ij}\cdot{\bf
v}_{ij}^{\ast})\hat{\bf r}_{ij},\quad\quad{\bf v}_{ij}^{\ast}={\bf v}_{i}-{\bf
v}_{j}+{\bf u}_{s}^{i},$ (11)

where ${\bf u}_{s}^{i}$ is the slip velocity at the position of squirmer
particle $i$, while $j$ corresponds to an outer-fluid DPD particle.

Excluded volume interactions between active squirmers and tracers are imposed
through the 12-6 Lennard-Jones (LJ) potential as

$U_{LJ}(r)=4\epsilon\left[\left(\frac{\sigma_{LJ}}{r}\right)^{12}-\left(\frac{\sigma_{LJ}}{r}\right)^{6}\right],$
(12)

where $\sigma_{LJ}$ is the distance at which $U_{LJ}(\sigma_{LJ})=0$ and
$\epsilon$ is the potential strength. Note that only the repulsive part of the
LJ potential is employed by setting the cutoff radius to
$r_{c}^{LJ}=2^{1/6}\sigma_{LJ}$.

Simulation setup and parameters. The simulation domain corresponds to a slit-
like geometry with a size of $100r_{c}\times 100r_{c}\times 22r_{c}$ and
periodic BCs in the xy plane. In the $z$ axis, the fluid layer of thickness
$20r_{c}$ is confined between two walls, each with a thickness of $r_{c}$. The
fluid modeled by a collection of DPD particles has a number density of
$n=5/r_{c}^{3}$ ($r_{c}=1.0$ is selected in simulations). The DPD force
coefficients for interactions between fluid particles are $b=50k_{B}T/r_{c}$,
$\gamma=14.15\sqrt{mk_{B}T}/r_{c}$, and $s=0.15$, where $m=1$ is the particle
mass and $k_{B}T=2$ is the unit of energy. These parameters lead to a fluid
dynamic viscosity of $\eta=16.7\sqrt{mk_{B}T}/r_{c}^{2}$.

Suspended squirmers are constructed from $1000$ particles placed on a
spherical surface with a diameter $d_{Sq}=6.4r_{c}$. The largest investigated
tracers have the same diameter. The area fraction of $\Phi=0.05$ corresponds
to $18$ simulated squirmers, to which $10$ tracers were added. Both the
squirmers and tracers are placed at the lower wall ($z=0$) with no-slip BCs,
and restricted to a fluid layer of thickness $7.3r_{c}$, in order to properly
capture the quasi-2D environment in experiments. At the upper wall
($z=20r_{c}$) slip BCs are imposed, which partially accounts for a larger
domain in experiments. Coupling between the squirmer/tracer particles and
fluid particles assumes $b=0$ and $\gamma=12\sqrt{mk_{B}T}/r_{c}$ with $s=0.1$
and $r^{\prime}_{c}=0.7r_{c}$. The value of $\gamma$ here is computed such
that the imposed slip BCs at the squirmer surface are obtained Fedosov  _et
al._ (2010). The time step in simulations is $\Delta
t=0.007r_{c}\sqrt{m/k_{B}T}$.

To relate simulation and physical units, we introduce characteristic length
and time scales. The length scale is based on the squirmer diameter
$d_{Sq}=6.4r_{c}$, resulting in $r_{c}$ to represent approximately $0.34$ $\mu
m$ ($d_{S}=2.16$ $\mu m$ in experiments). Therefore, the simulation domain
corresponds to $34\times 34\times 6.8$ $\mu m$. The time scale is
$\tau=d_{Sq}/v_{M}$, where $v_{M}$ is the modal velocity of the velocity
distribution of active particles in the xy plane, see Fig. 1(b) in the main
text.

Figure. S 1: Separating “clustered” particles from “free” particles. a.
Identification of “clustered” particles from micrographs. Clusters are defined
as particles with 2 or more nearest neighbours within 2.5R for at least 12
frames (0.24 s, Theurkauff  _et al._ (2012)). b. Mean velocity of free (open
circle) and clustered (closed diamond) microswimmers. c. Onset of clustering
at an area fraction $\Phi\sim 0.05$. d. Instantaneous velocity distributions
of free and clustered particles (symbols as before). The red curves (solid and
dashed for the free and clustered particles respectively) show the best fits
of the Maxwell distribution Wu and Libchaber (2000). We clearly see that
clustered particles cannot be well explained by this distribution, while the
motion of the non-interacting microswimmers is well captured. Notably, the
distribution of velocities follows that for neutral squirmers (see main text,
Figure 4). Figure. S 2: Effect of “actively transported” colloids on the
ensemble-mean-squared-displacement (EMSD) of all tracers. Left: The individual
MSDs (IMSDs) of multiple tracer particles (here, $1\mu m$) in a dilute
suspension of synthetic microswimmers. We see one extreme outlier, which
corresponds to a singular tracer particle being actively transported by a
microswimmer (via phoretic attraction). Right: The single actively transported
tracer significantly contributes to the overall EMSD (red - due to averaging
of all values). Removing this single IMSD from the ensemble (blue) collapses
the EMSD of the tracers onto the control EMSD (green) of passive tracers in
the absence of active microswimmers. This highlights the issues of blindly
using averaged values alone when evaluating systems with large outliers.
Figure. S 3: Extracting the experimental
$\tau_{\gamma=\operatorname*{argmax}_{\tau}(\gamma)}$ from experiments
performed under different conditions (see Figure 2 b. of the main text).
$\gamma(\tau)$ is extracted from the displacements of the tracers in different
experiments (varying tracer and swimmer concentration), and a smoothing spline
is then fitted to extract $\tau$ at which the value of $\gamma$ is maximised
Fernandez-Rodriguez  _et al._ (2020).  Figure. S 4: Tracer displacements in
the presence of microswimmers which are passively diffusing (no UV
illumination). At both short and long lag times $\Delta\tau$, the distribution
is Gaussian. Contrast this to the enhanced diffusivity of tracers in the
presence of motile microswimmers, shown in Figure 3 (top row) of the main
text. Note that the autocorrelation approach to extracting jumps also
identifies jumps in the passive case, but that the resulting distributions are
unaffected.  Figure. S 5: Demonstration of the phoretic interactions between
the Janus microswimmers and the PS tracer nanoparticles. By immobilizing our
microswimmers and modulating the pH, we also observe a phoretic contribution
to the generated flow fields, in agreement with previous reports for TiO2-SiO2
systems Singh  _et al._ (2017). The values are averaged over at least 4 Janus
pumps (stuck to the substrate via drying) in each case. Specifically, we note
the change in the accumulation of the tracer nanoparticles as
tetramethylammonium hydroxide (TMAH) is added to modulate the pH, following
Singh  _et al._ (2017). a. Accumulation of PS nano-tracers around a fixed
Janus particle at the SiO2 hemisphere under standard experimental conditions
in the presence of H2O2 with UV illumination (pH $\sim 4$, orientation shown
in inset, top left). Colour coding indicates the normalised number density, as
well as the averaged profiles in X and Y. b. Distribution of displacements of
the nanoparticle tracers in the regions marked I and II in a. We see that near
the Janus particle, tracer displacements decrease due to accumulation. Far
from the particle, we observe the expected Gaussian distribution (also see
Figure S5). c. Accumulation of PS nano-tracers around a fixed Janus particle
at the TiO2 hemisphere under modulation of the pH (pH $\sim 7$) via addition
of TMAH, in the presence of H2O2 with UV illumination (orientation shown in
inset, top left). For these experiments, tracers were visualised using a 100x
oil objective. Figure. S 6: Average instantaneous displacement of tracer
nanoparticles in rings evaluated around the Janus pumps evaluated in Figure S4
a. We see that outside of the inner rings (where the tracers accumulate due to
phoretic interactions), the average displacements of the tracers rapidly
saturate, indicating little-to-no effect of the hydrodynamic flows created by
the microswimmer on tracers at a distance.  Figure. S 7: Fitted translational
diffusion coefficients $D_{T}$ of tracer nanoparticles in the presence of
different concentrations (represented by area fraction, $\Phi$) of
microswimmers, normalised by the control at $\Phi=0$. The mean value and error
bars are obtained from the distribution of $D_{T}$ fitted from the mean
logarithmic squared displacementKepten  _et al._ (2013); Bailey  _et al._
(2022e) (MLSD, given by $MLSD=log(4D_{T})+\alpha log(\tau)$), assuming a log-
normal distribution (also see Figure S7 c.). We see that there is no clear
trend in the values of $D_{T}$ with the presence of microswimmers, indicating
the lack of an enhancement in their diffusivity (via mixing), which we
attribute to their high initial diffusion coefficient (arising from their
small size).  Figure. S 8: Comparison of the fitting approaches used to
extract $D_{T}$ values of the (control, i.e. in the absence of microswimmers)
PS nanoparticle tracers, to confirm the accuracy of various approaches. a-b.
Fitting of $D_{T}$ via differential dynamic microscopy (DDM Cerbino and Trappe
(2008)). c. Fitting a log-normal distribution to the binned values of $D_{T}$
obtained using the MLSD of individual tracers (also see Figure S6) d. Fitting
the EMSD of the tracers directly. e. Comparison of the fits obtained using
different approaches (in $\mu m^{2}s^{-1}$). We see good agreement between
using DDM and the distribution of the MLSDs of individual nanoparticle tracers
in the passive case. As the presence of microswimmers (additional scattering
sources) affects the interpretation of DDM results, we therefore use the
values obtained from the distribution of the IMSDs.


# Mapping Housing Stock Characteristics from Drone Images for Climate
Resilience in the Caribbean

Isabelle Tingzon, Nuala Margaret Cowan, and Pierre Chrzanowski
The World Bank Group, GFDRR
{tisabelle, ncowan<EMAIL_ADDRESS>
###### Abstract

Comprehensive information on housing stock is crucial for climate adaptation
initiatives aiming to reduce the adverse impacts of climate-extreme hazards in
high-risk regions like the Caribbean. In this study, we propose a workflow for
rapidly generating critical baseline housing stock data using very high-
resolution drone images and deep learning techniques. Specifically, our work
leverages the Segment Anything Model and convolutional neural networks for the
automated generation of building footprints and roof classification maps. By
strengthening local capacity within government agencies to leverage AI and
Earth Observation-based solutions, this work seeks to improve the climate
resilience of the housing sector in small island developing states in the
Caribbean.

##  1 Introduction

The Caribbean is among the world’s most climate-vulnerable regions due to the
prevalence and intensity of extreme climate hazards such as storms, floods,
and landslides. Category 5 hurricanes like Dorian, Irma, and Maria have
devastated many small island developing states (SIDS) in recent years, leaving
widespread trails of loss and destruction across the region. SIDS bear
substantial economic costs from climate-extreme events, with the highest
degree of damages often sustained in the housing sector [19]. Hurricane Maria,
for example, destroyed over 90% of Dominica’s housing stock, accumulating
costs over 200% of the nation’s GDP [18]. As global temperatures continue to
rise, extreme weather events will only grow in severity, putting many more
vulnerable homes and shelters at risk.

In response to these challenges, governments and international organizations
have developed ambitious climate resilience programs to reduce the adverse
effects of extreme climatic hazards in the housing sector [19, 9]. Climate
resilience initiatives generally require comprehensive housing stock
information to inform better retrofitting, reconstruction, and relocation
plans. However, the traditional house-to-house approach to identifying high-
risk buildings can be extremely expensive and time-consuming to implement,
prompting the need for more timely and cost-efficient alternatives.

Recent years have seen a growing interest in using drones for disaster risk
reduction and recovery [9, 28]. Previous works have successfully applied deep
learning (DL) techniques to very high-resolution (VHR) aerial imagery to
extract baseline housing information such as building footprints and rooftop
attributes [27, 10, 22, 14, 24]. However, despite the evident advantages DL
and Earth observation (EO)-based technologies, the widespread adoption of
these solutions is often hindered by gaps in local capacity to develop and
maintain systems for generating baseline exposure datasets.

This work aims to bridge these gaps by providing government agencies with an
end-to-end workflow for rapidly generating critical baseline housing
information using VHR drone images. Specifically, we leverage the Segment
Anything Model (SAM) for building footprint delineation and convolutional
neural networks (CNNs) for roof type and roof material classification. We also
evaluate the cross-country generalizability of roof classification models
across SIDS to determine the extent to which models trained in one country can
be adapted to another. This work is developed under the Digital Earth Project
for Resilient Housing and Infrastructure in partnership with the Government of
the Commonwealth of Dominica (GoCD) and the Government of Saint Lucia (GoSL)
to enhance the climate resilience of the housing sector of SIDS in the
Caribbean.

(a) Healthy Metal

(b) Irregular Metal

(c) Blue Tarpaulin

(d) Concrete/cement

(e) Incomplete

(a) Hip

(b) Gable

(c) Flat

(d) No roof

Figure 1: Examples of VHR drone-derived roof image tiles for each of the roof
material categories (top row) and roof type categories (bottom row).

##  2 Data

To generate our ground truth datasets, we use (1) VHR aerial images in the
form of aircraft- and drone-derived optical imagery and (2) building
footprints in vector polygons, with Dominica and Saint Lucia as our primary
regions of interest.

VHR Aerial Images. We acquired the following VHR aerial images from partner
government agencies GoCD and GoSL, the World Bank Global Program for Resilient
Housing (GPRH) [8], and the open data platform OpenAerialMap [5]: (1)
aircraft-derived post-disaster orthophotos of Dominica; (2) aircraft-derived
orthophotos of Saint Lucia; (3) pre- and post-disaster drone images of 10
villages and cities across Dominica; (4) drone images of 3 districts in Saint
Lucia. The spatial resolutions of aircraft-derived images range from 10 to 20
cm/px, whereas the resolutions of drone-derived images range from 2 to 5
cm/px. For more information on the spatial resolution, coverage, and year of
acquisition of the aerial images, see Table 3 in the Appendix.

Building Footprints Data. We obtained nationwide building footprints
delineated from the aircraft-derived orthophotos of Dominica and Saint Lucia
from the World Bank. For drone images with no corresponding building
footprints, we initially looked to alternative data sources such as
OpenStreetMap [6], Microsoft [4], and Google [2]. However, we found these
publicly available building footprint datasets to be nonviable due to
significant misalignment with the underlying drone images, as illustrated in
Figure 2. To address this challenge, we used the Segment Anything Model (SAM)
to delineate building instances directly from drone images [15]. Additional
information on the SAM configuration used in this study is detailed in Section
3.

(a) Microsoft

(b) Google

(c) OpenStreetMap

(d) SAM

Figure 2: Building footprints from (a) Microsoft, (b) Google, (c)
OpenStreetMap, and (d) SAM overlaid on a drone image taken in Salisbury,
Dominica from OpenAerialMap [5].

Rooftop Image Tiles. For roof classification, we selected approximately 15,000
building footprints across Dominica and Saint Lucia. For each building
footprint, we extract the minimal bounding rectangle of the building polygon,
scaled by a factor of 1.5, from the corresponding aerial image. We then
proceeded to annotate a total of 9,150 buildings in Dominica and 5,175
buildings in Saint Lucia via visual interpretation of VHR aerial images. The
buildings are classified based on two attributes: roof type and roof material,
the class distributions of which are presented in Table 1. We note that Saint
Lucia does not contain any samples for the "Blue Tarpaulin" category. Figure 1
illustrates examples of drone image tiles for each roof type and roof material
category.

Table 1: The distribution of roof type and roof material labels across Dominica and Saint Lucia. |  | Dominica | Saint Lucia
---|---|---|---
|  | Train | Test | Total | Train | Test | Total
Roof Type |  Gable | 2,669 | 653 | 3,322 | 2,347 | 585 | 2,932
Hip | 1,579 | 393 | 1,972 | 1,089 | 271 | 1,360
Flat | 1,894 | 475 | 2,369 | 456 | 106 | 562
No Roof | 1,190 | 297 | 1,487 | 269 | 52 | 321
Roof Material |  Healthy metal | 1,934 | 482 | 2,416 | 2,396 | 598 | 2,994
Irregular metal | 1,733 | 432 | 2,165 | 1,113 | 276 | 1,389
Concrete/cement | 1,240 | 312 | 1,552 | 328 | 75 | 403
Blue tarpaulin | 1,094 | 260 | 1,354 | 0 | 0 | 0
Incomplete | 1,331 | 332 | 1,663 | 324 | 65 | 389
Source |  Aircraft | 5,936 | 0 | 5,936 | 2,485 | 0 | 2,485
|  Drone | 1,396 | 1,818 | 3,214 | 1,676 | 1,014 | 2,690
|  Total | 7,332 | 1,818 | 9,150 | 4,161 | 1,014 | 5,175

##  3 Methods

This section outlines our workflow for generating critical housing stock
information in the Caribbean using drone images, as summarized in Figure 3.

Figure 3: Proposed workflow for the automatic generation of housing stock
information from drone images using DL models.

To extract building footprints from drone images, we use Segment-Geospatial
for segmenting raster images based on the Segment Anything Model (SAM) [29,
20, 21, 15]. We also leverage Language Segment Anything [16] to combine
instance segmentation with text prompts to generate masks of specific objects
in the drone images. We set our text prompt to "house", the box threshold
(i.e. the threshold value used for object detection in the image) to 0.30, and
the text threshold (i.e. the threshold value used to associate the detected
objects with the provided text prompt) to 0.30. As a post-processing step, we
apply the Douglas-Peucker algorithm to simplify the generated building
polygons with a tolerance of $5e^{-6}$ [12].

Given the VHR aerial images and corresponding building footprint polygons, we
proceed with developing our roof classification models. We begin by fine-
tuning CNN models pre-trained on the ImageNet dataset [11] with architectures
ResNet50 [13], VGG-16 [23], Inceptionv3 [25], and EfficientNet-B0 [26] using
cross-entropy loss for both roof type and roof material classification tasks.
The input rooftop image tiles are zero-padded to a square of size 224 x 224 px
for ResNet50, EfficientNet-B0, and VGG-16 and 299 x 299 px for InceptionV3. We
set the batch size to 32 and the maximum number of epochs to 60, and we use an
Adam optimizer with an initial learning rate of $1e^{-}5$, which decays by a
factor of 0.1 after every 7 epochs with no improvement. For data augmentation,
we implement horizontal and vertical image flips with a probability of 0.50
and random rotations ranging from $-90^{\circ}$ to $90^{\circ}$. Given that
our data is imbalanced (see Table 1), we implement random oversampling for the
minority classes. To prevent overconfident predictions, we apply label
smoothing as a regularization technique, with smoothing set to 0.1 [17].

##  4 Results and Discussion

For each country in our dataset, we set aside 80% of the data for training and
the remaining 20% for testing, using stratified random sampling for the test
set to preserve the percentage of samples per class as shown in Table 1. We
note that the test sets for both countries are comprised entirely of drone
images. Additionally, to test whether geographically diverse training data
improves the prediction, we combine the training sets across Dominica and
Saint Lucia (henceforth referred to as the "combined" dataset). We report the
standard performance metrics F1-score, precision, recall, and accuracy.

Our results indicate that for Dominica, the best performance is achieved by an
EfficientNet-B0 model for roof type classification (F1-score: 87.1%) and a
ResNet50 model for roof material classification (F1-score: 89.5%). Likewise,
for Saint Lucia, the best F1-score is attained by a ResNet50 model for roof
type classification (F1-score: 89.5%) and an EfficientNet-B0 model for roof
material classification (F1-score: 91.71%). For models trained on the combined
dataset (i.e., "combined model"), the best results are obtained by an
EfficientNet-B0 model for roof type classification (F1-score: 90.0%) and an
Inceptionv3 model for roof material classification (F1-score: 90.4%). For the
complete results and sample outputs, see Table 4 and Figure 4, respectively.

(a) 2017 drone image from GoCD

(b) 2018 drone image from OpenAerialMap

(c) 2017 roof material classification map

(d) 2018 roof material classification map

Figure 4: Drone images (top) and corresponding roof material classification
maps (bottom) of Coulibistrie, Dominica taken before (left) and after (right)
Hurricane Maria.

Next, we investigate the cross-country generalizability of the best-performing
models by evaluating their performance on the designated test sets of each
country. As shown in Table 2, our results indicate that for roof type
classification, the combined model performs marginally better than models
trained using only local data (i.e. data from the same country); however, for
roof material classification, we find that locally trained models consistently
outperform the combined model. In general, we find that local models
demonstrate high levels of variability in performance in out-of-distribution
countries, indicating the importance of collecting localized training data and
the need for further studies on domain adaptation to reduce performance
degradation in the face of geographical distribution shifts.

Table 2: Cross-country generalization of roof type and roof material
classification models.

(a) Roof Type |  |
---|---|---
Training Data | Test Data | F1 score | Precision | Recall | Accuracy
Dominica | Dominica | 87.08 | 87.06 | 87.37 | 86.41
Saint Lucia | 88.39 | 88.93 | 87.92 | 93.04
Saint Lucia | Dominica | 82.59 | 84.72 | 81.39 | 82.31
Saint Lucia | 89.49 | 94.11 | 85.97 | 94.28
Combined | Dominica | 88.02 | 88.51 | 88.01 | 87.39
Saint Lucia | 91.92 | 93.63 | 90.36 | 95.67

(a) Roof Material |  |
---|---|---
Training Data | Test Data | F1 score | Precision | Recall | Accuracy
Dominica | Dominica | 89.50 | 90.49 | 88.84 | 89.44
Saint Lucia | 90.83 | 92.04 | 89.93 | 92.94
Saint Lucia | Dominica | 64.15 | 63.14 | 67.32 | 72.99
Saint Lucia | 91.71 | 93.81 | 90.00 | 93.79
Combined | Dominica | 88.17 | 89.66 | 87.16 | 87.78
Saint Lucia | 90.37 | 93.10 | 88.58 | 93.55

##  5 Local Capacity Building

The overarching goal of this project is to strengthen local capacity in SIDS
to leverage AI and EO-based solutions for resilient housing operations. To
this end, our team is assisting government agencies in establishing Geographic
Information Systems (GIS) units capable of generating, managing, and
maintaining large-scale disaster risk datasets. We have partnered with the
Humanitarian OpenStreetMap Team (HOT) [3] to design a training program on how
to operate drones, coordinate pilots, process collected imagery, and manage
the resulting geospatial datasets. Additionally, we have developed educational
resources on geospatial data processing aimed at local government staff,
community mappers, and disaster responders, including Google Colaboratory and
Jupyter Notebook tutorials demonstrating how to run the SAM and CNN models to
quickly generate baseline housing stock data from locally collected drone
images [1, 7].

##  6 Conclusion

This work proposes a workflow for filling baseline exposure data gaps in the
Caribbean using VHR drone images and DL. We demonstrate how computer vision
and drone-based technologies can be used to rapidly generate housing stock
information, including building footprints and roof classification maps, for
disaster risk reduction and recovery. Based on our evaluation of the cross-
country generalizability of DL models, we urge caution in applying locally-
trained models off the shelf to new geographic regions and emphasize the
importance of collecting local, highly contextualized training data. We also
emphasize the importance of local capacity building, skills development, and
co-creation of geospatial datasets in deploying sustainable AI-for-climate
solutions, especially in Global South contexts. We conclude this study by
urging governments in SIDS to invest in the digital infrastructure and local
capabilities needed to sustainably generate and maintain EO-derived housing
stock data for climate resilience.

## Acknowledgments

This work is developed under Digital Earth for Resilient Housing and
Infrastructure in the Caribbean, a World Bank project funded by the Global
Facility for Disaster Reduction and Recovery (GFDRR), in partnership with the
Government of the Commonwealth of Dominica (GoCD) and the Government of Saint
Lucia (GoSL). This study builds on the initial work done by the World Bank’s
Global Program for Resilient Housing (GPRH). We thank Mike Fedak and
Christopher Williams for their assistance in providing data access.

## References

  * [1] Building Footprint Delineation for Disaster Risk Reduction and Response (Part I).  https://colab.research.google.com/github/GFDRR/caribbean-rooftop-classification/blob/master/tutorials/01_building_delineation.ipynb.
  * [2] Google Open Buildings.  https://sites.research.google/open-buildings/.  Accessed on 15.09.2023.
  * [3] Humanitarian OpenStreetMap Team (HOT).  https://www.hotosm.org/.
  * [4] Microsoft Building Footprints.  https://www.microsoft.com/en-us/maps/building-footprints.  Accessed on 15.09.2023.
  * [5] OpenAerialMap.  https://map.openaerialmap.org/.  Accessed on 15.09.2023.
  * [6] OpenStreetMap Buildings.  https://osmbuildings.org/.  Accessed on 15.09.2023.
  * [7] Rooftop Classification from Drone Imagery for Disaster Risk Reduction and Reponse (Part II).  https://colab.research.google.com/github/GFDRR/caribbean-rooftop-classification/blob/master/tutorials/02_building_classification.ipynb.
  * [8] Global Program for Resilient Housing (GPRH).  https://www.worldbank.org/en/topic/disasterriskmanagement/brief/global-program-for-resilient-housing, January 2022.
  * [9] World Bank.  Capturing Housing Data in Small Island Developing States, 2022.  Washington, DC. Creative Commons Attribution CC BY 4.0.
  * [10] M Buyukdemircioglu, R Can, and S Kocaman.  Deep learning based roof type classification using very high resolution aerial imagery.  The International Archives of Photogrammetry, Remote Sensing and Spatial Information Sciences, 43:55–60, 2021.
  * [11] Jia Deng, Wei Dong, Richard Socher, Li-Jia Li, Kai Li, and Li Fei-Fei.  Imagenet: A large-scale hierarchical image database.  In 2009 IEEE conference on computer vision and pattern recognition, pages 248–255. Ieee, 2009.
  * [12] David H Douglas and Thomas K Peucker.  Algorithms for the reduction of the number of points required to represent a digitized line or its caricature.  Cartographica: the international journal for geographic information and geovisualization, 10(2):112–122, 1973.
  * [13] Kaiming He, Xiangyu Zhang, Shaoqing Ren, and Jian Sun.  Deep residual learning for image recognition.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 770–778, 2016.
  * [14] Xingliang Huang, Libo Ren, Chenglong Liu, Yixuan Wang, Hongfeng Yu, Michael Schmitt, Ronny Hänsch, Xian Sun, Hai Huang, and Helmut Mayer.  Urban Building Classification (UBC)-A Dataset for Individual Building Detection and Classification From Satellite Imagery.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 1413–1421, 2022.
  * [15] Alexander Kirillov, Eric Mintun, Nikhila Ravi, Hanzi Mao, Chloe Rolland, Laura Gustafson, Tete Xiao, Spencer Whitehead, Alexander C Berg, Wan-Yen Lo, et al.  Segment anything.  arXiv preprint arXiv:2304.02643, 2023.
  * [16] Luca Medeiros.  Language Segment-Anything.  https://github.com/luca-medeiros/lang-segment-anything, 2023.
  * [17] Rafael Müller, Simon Kornblith, and Geoffrey E Hinton.  When does label smoothing help?  Advances in neural information processing systems, 32, 2019.
  * [18] Government of the Commonwealth of Dominica.  Post-Disaster Needs Assessment Hurricane Maria September 18, 2017.  https://www.gfdrr.org/sites/default/files/publication/Dominica_mp_012418_web.pdf, September 2017.
  * [19] Government of the Commonwealth of Dominica.  Dominica Climate Resilience and Recovery Plan 2020–2030.  https://odm.gov.dm/wp-content/uploads/2022/02/CRRP-Final-042020.pdf, 2020\.
  * [20] Lucas Osco.  AI-RemoteSensing: A collection of Jupyter and Google Colaboratory notebooks dedicated to leveraging Artificial Intelligence (AI) in Remote Sensing applications, June 2023.
  * [21] Lucas Prado Osco.  AI-RemoteSensing.  https://github.com/LucasOsco/AI-RemoteSensing, 2023.
  * [22] Tahmineh Partovi, Friedrich Fraundorfer, Seyedmajid Azimi, Dimitrios Marmanis, and Peter Reinartz.  Roof Type Selection based on patch-based classification using deep learning for high-resolution Satellite Imagery.  International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences-ISPRS Archives, 42(W1):653–657, 2017.
  * [23] Karen Simonyan and Andrew Zisserman.  Very deep convolutional networks for large-scale image recognition.  arXiv preprint arXiv:1409.1556, 2014.
  * [24] Roman A Solovyev.  Roof material classification from aerial imagery.  Optical Memory and Neural Networks, 29:198–208, 2020.
  * [25] Christian Szegedy, Vincent Vanhoucke, Sergey Ioffe, Jon Shlens, and Zbigniew Wojna.  Rethinking the inception architecture for computer vision.  In Proceedings of the IEEE conference on computer vision and pattern recognition, pages 2818–2826, 2016.
  * [26] Mingxing Tan and Quoc Le.  Efficientnet: Rethinking model scaling for convolutional neural networks.  In International conference on machine learning, pages 6105–6114. PMLR, 2019.
  * [27] Isabelle Tingzon, Nuala Margaret Cowan, and Pierre Chrzanowski.  Fusing VHR Post-disaster Aerial Imagery and LiDAR Data for Roof Classification in the Caribbean using CNNs.  arXiv preprint arXiv:2307.16177, 2023.
  * [28] Luis Triveno, Sarah Antos, Jan Koers, and Victor Endo.  Coupling Imagery from Drones and Street-View with Proper Incentives To Promote Sustainable Urban Cadasters in Developing Countries.  In Proceedings of the World Bank Conference on Land and Poverty, Washington, DC, USA, pages 25–29, 2019.
  * [29] Qiusheng Wu and Lucas Prado Osco.  samgeo: A Python package for segmenting geospatial data with the Segment Anything Model (SAM).  Journal of Open Source Software, 8(89):5663, 2023.

## Appendix

Table 3: Aircraft- and drone-derived aerial images used in this study. Coverage | Resolution | Year | Source | Data Provider |  | Building
---
Count
| Dominica
---
(nationwide)
20.0 cm/px | 2018-2019 | Aircraft | GoCD | 5,936
Colihaut | 2.7 cm/px | 2017 | Drone | GoCD | 373
Coulibistrie | 2.3 cm/px | 2017 | Drone | GoCD | 158
Delices | 4.3 cm/px | 2018 | Drone | OpenAerialMap [5] | 380
Dublanc | 2.9 cm/px | 2017 | Drone | GoCD | 126
Kalinago | 3.3 cm/px | 2018 | Drone | OpenAerialMap [5] | 102
Laplaine | 4.9 cm/px | 2018 | Drone | OpenAerialMap [5] | 456
Marigot | 3.4 cm/px | 2018 | Drone | OpenAerialMap [5] | 387
Pichelin | 3.4 cm/px | 2017 | Drone | GoCD | 149
Roseau | 2.6 cm/px | 2017 | Drone | GoCD | 348
Salisbury | 3-5 cm/px | 2018 | Drone | OpenAerialMap [5] | 280
| Saint Lucia
---
(nationwide)
10.0 cm/px | 2022 | Aircraft | GoSL | 2,485
Castries | 4.5 cm/px | 2019 | Drone | GPRH [8] | 1,084
Dennery | 4.2 cm/px | 2019 | Drone | GPRH [8] | 742
Gros Islet | 3.6 cm/px | 2019 | Drone | GPRH [8] | 864
Table 4: Test set (%) results of different CNN architectures for roof type and
roof material classification trained using (a) only Dominica data, (b) only
Saint Lucia data, and (c) using a combination of Dominica and Saint Lucia.
Models trained on the combined dataset are evaluated using the combined test
sets of Dominica and Saint Lucia.

(a) Dominica |  |
---|---|---
|  | F1 score | Precision | Recall | Accuracy
Roof Type | VGG16 | 86.15 | 86.86 | 85.79 | 85.37
ResNet50 | 86.16 | 85.80 | 86.70 | 85.70
Inceptionv3 | 86.53 | 86.93 | 86.44 | 85.70
EfficientNet-B0 | 87.08 | 87.06 | 87.37 | 86.41
Roof Material | VGG16 | 89.21 | 90.11 | 88.73 | 88.89
| ResNet50 | 89.50 | 90.49 | 88.84 | 89.44
| Inceptionv3 | 88.24 | 89.72 | 87.36 | 88.17
| EfficientNet-B0 | 89.00 | 90.20 | 88.23 | 88.94

(b) Saint Lucia |  |
---|---|---
|  | F1 score | Precision | Recall | Accuracy
Roof Type | VGG16 | 88.02 | 91.29 | 85.36 | 93.10
ResNet50 | 89.49 | 94.11 | 85.97 | 94.28
Inceptionv3 | 87.50 | 94.32 | 83.02 | 93.45
EfficientNet-B0 | 88.10 | 94.77 | 83.51 | 93.29
Roof Material | VGG16 | 88.27 | 92.38 | 85.41 | 91.52
| ResNet50 | 91.55 | 92.99 | 90.28 | 93.20
| Inceptionv3 | 91.42 | 93.28 | 89.91 | 93.15
| EfficientNet-B0 | 91.71 | 93.81 | 90.00 | 93.79

(c) Combined (Dominica + Saint Lucia) |  |
---|---|---
|  | F1 score | Precision | Recall | Accuracy
Roof Type | VGG16 | 88.70 | 88.70 | 88.72 | 89.03
ResNet50 | 88.71 | 88.63 | 88.87 | 89.24
Inceptionv3 | 89.38 | 89.29 | 89.62 | 89.60
EfficientNet-B0 | 90.03 | 90.07 | 90.25 | 90.45
Roof Material | VGG16 | 90.18 | 92.09 | 88.72 | 90.73
| ResNet50 | 89.82 | 91.76 | 88.36 | 90.80
| Inceptionv3 | 90.42 | 92.03 | 89.10 | 90.80
| EfficientNet-B0 | 90.09 | 92.37 | 88.45 | 91.09


Prox}_{\lambda(\mu_{0}\|\cdot\|_{0})}\big{(}x^{k}-\lambda\nabla
f(x^{k}\big{)},\quad\widehat{v}^{k}=\nabla f(\widehat{x}^{k})-\nabla
f(x^{k})+\frac{1}{\lambda}(x^{k}-\widehat{x}^{k}).$ (7.18)

Thus $d^{k}$ can be computed for each $k\in{\rm I\\!N}$ by solving the system
of linear equations $X^{k}d=y^{k}$, where

$(X^{k})_{i}:=\begin{cases}\big{(}\nabla^{2}f(x^{k})\big{)}_{i}&\text{if}\quad\left(\widehat{x}^{k}\right)_{i}\neq
0,\\\
I_{i}&\text{if}\quad\left(\widehat{x}^{k}\right)_{i}=0\end{cases}\quad\text{and
}\;y^{k}:=\begin{cases}(-\widehat{v}^{k})_{i}&\text{if}\quad\left(\widehat{x}^{k}\right)_{i}\neq
0,\\\ 0&\text{if}\quad\left(\widehat{x}^{k}\right)_{i}=0.\end{cases}$ (7.19)

######  Theorem 7.4 (solving the $\ell_{0}$-$\ell_{2}$ regularized least
square problem by GCNM).

Consider the $\ell_{0}$-$\ell_{2}$ model (7.1) with notation (7.3) for $f$ and
$g$. The following assertions hold:

(i) The directions $d^{k}$ satisfying (6.1) in the Newton step of Algorithm 4
always exist and are calculated from the conditions in (7.17)–(7.19).

(ii) Algorithm 4, with the data for $\nabla f$, $\nabla^{2}f$,
$\partial^{2}g$, and $\text{\rm Prox}_{\lambda g}$ that are explicitly
computed in (7.4), (7.5), and (7.7), either stops after finitely many
iterations, or produces a sequence $\\{x^{k}\\}$ such that any accumulation
point $\bar{x}$ of this sequence is a critical point of (7.1).

(iii) If the accumulation point $\bar{x}$ satisfies the property

$\left(A^{*}(A\bar{x}-b)+2\mu_{2}\bar{x}\right)_{i}\neq 0\quad\text{whenever
}\;\bar{x}_{i}=0,$

then the sequence $\\{x^{k}\\}$ globally Q-superlinearly converges to
$\bar{x}$ with the Q-superlinear rate.

[Proof.] To verify (i), we deduce from Corollary 7.2 the implication

$0\in(D\partial\varphi)(x,y)(v),v\in{\rm I\\!R}^{n}\Longrightarrow v=0\;\mbox{
for any }\;(x,y)\in\mbox{\rm gph}\,\partial\varphi,$

which implies that $\partial\varphi$ is strongly metrically subregular at
$(x,y)$ due to the Levy–Rockafellar criterion. Thus (i) holds by Lemma 3.1. To
proceed with (ii), observe that since $g$ is bounded from below by $0$, it
follows from [35, Exercise 1.24] that $g$ is prox-bounded with threshold
$\lambda_{g}=\infty$. By the calculation in (7.4), $\nabla f$ is Lipschitz
continuous on ${\rm I\\!R}^{n}$ with modulus $L_{f}=\lambda_{\text{\rm
max}}(A^{*}A)+2\mu_{2}$ expressed via the largest eigenvalue of $A^{*}A$.
Therefore, Assumption 2 is also guaranteed for (7.1). Applying Theorem 6.1
tells us that any accumulation point of $\\{x^{k}\\}$ is a critical point of
(7.1), which justifies (ii).

Let us finally verify (iii). We know from (ii) that the accumulation point
$\bar{x}$ is critical for (7.1), which implies that
$0\in\partial\varphi(\bar{x})=\nabla f(\bar{x})+\partial g(\bar{x})$, i.e.,
$-\nabla f(\bar{x})\in\partial g(\bar{x})$. It follows from Lemma 7.3 and the
formula $\nabla f(\bar{x})=A^{*}(A\bar{x}-b)+2\mu_{2}\bar{x}$ that $g$ is
subdifferentially continuous at $\bar{x}$ for $-\nabla f(\bar{x})$. The prox-
regularity of $g$ at $\bar{x}$ for $-\nabla f(\bar{x})$ is a consequence of
the variational convexity of $g$ at $\bar{x}$ for $-\nabla f(\bar{x})$. It
further follows from Corollary 7.2 that $\partial^{2}\varphi(\bar{x},0)$ is
positive-definite, which ensures that $\bar{x}$ is a tilt-stable local
minimizer of $\varphi$ by [32, Theorem 1.3]. The graph of $\partial g$ is the
union of finitely many closed convex sets, and hence $\partial g$ is
semismooth∗ on its graph by [9]. All of this allow us to deduce from Theorem
6.1 the $Q$-superlinear convergence of the iterative sequence $\\{x^{k}\\}$
generated by Algorithm 4 to the solution $\bar{x}$ and thus complete the proof
of the theorem.

Having the obtained results and calculations for solving the
$\ell_{0}$-$\ell_{2}$ model (7.1), we conduct numerical experiments for our
GCNM Algorithm 4 and compare them with those obtained by implementing the
proximal gradient method (PGM) (also known as the forward-backward splitting
method) developed in [2, Section 5.1]. In the conducted numerical experiments,
the matrix $A$ and vector $b$ are randomly generated with the entries chosen
from a uniform distribution of random numbers. The initial points in all the
experiments are set to be $x^{0}:=1_{{\rm I\\!R}^{n}}:=(1,\ldots,1)$. Each
method is terminated at the iterate $x^{k}$ if $\|x^{k}-\widehat{x}^{k}\|\leq
10^{-6}$. We also stop the algorithms if they reach the maximum number of
iterations of $50,000$. To simplify the numerical implementations for solving
problem (7.1), the regularization parameters are chosen as either
$\mu_{0}=\mu_{2}=10^{-1}$, or $\mu_{0}=\mu_{2}=10^{-2}$. The achieved
numerical results are presented in Tables 1 and 2. In these tables, “CPU time”
stands for the time needed to achieve the prescribed accuracy of approximate
solutions (the smaller the better). As we can see from the results presented
in these tables, our GCNM algorithm achieves very small residuals and requires
a shorter computational time in comparison with PGM.

Test | Size Problem | Number of iterations | CPU time | Residual
---|---|---|---|---
m | n | GCNM | PGM | GCNM | PGM | GCNM | PGM
1 | 10 | 100 | 3 | 10403 | 0.0174 | 1.3522 | 1.88E-17 | 9.99E-07
2 | 100 | 1000 | 5 | 50000 | 0.3554 | 13.8951 | 7.85E-18 | 7.22E-06
3 | 1000 | 100 | 2 | 2881 | 0.0181 | 1.6866 | 9.34E-18 | 1.00E-06
4 | 100 | 10 | 2 | 374 | 0.002 | 0.0598 | 0 | 9.85E-07
5 | 100 | 100 | 2 | 9091 | 0.0087 | 1.7437 | 2.94E-17 | 1.00E-06
6 | 1000 | 1000 | 3 | 50000 | 0.3789 | 85.0849 | 2.25E-17 | 2.03E-06
Table 1: $\mu_{0}=\mu_{2}=10^{-1}$ Test | Size Problem | Number of iterations | CPU time | Residual
---|---|---|---|---
m | n | GCNM | PGM | GCNM | PGM | GCNM | PGM
1 | 10 | 100 | 3 | 50000 | 0.0038 | 2.8955 | 8.14E-18 | 2.88E-06
2 | 100 | 1000 | 2 | 7852 | 0.2281 | 4.9574 | 5.90E-18 | 1.00E-06
3 | 1000 | 100 | 2 | 3096 | 0.0177 | 1.3063 | 7.76E-18 | 1.00E-06
4 | 100 | 10 | 2 | 436 | 0.007 | 0.0563 | 1.96E-17 | 9.93E-07
5 | 100 | 100 | 1 | 11978 | 0.0069 | 2.7534 | 6.07E-16 | 1.00E-06
6 | 1000 | 1000 | 2 | 50000 | 0.2371 | 81.676 | 1.68E-17 | 2.38E-06
Table 2: $\mu_{0}=\mu_{2}=10^{-2}$

##  8 Conclusion

In this paper, we propose and develop fast converging local and global
generalized Newton methods to solve problems of nonconvex and nonsmooth
optimization. The developed algorithms are far-going extensions of the
classical Newton method, where the standard Hessian is replaced by its
generalized version specifically designed for nonsmooth (of first-order and
second-order) functions. This construction is coderivative-generated, giving
rise to the names of our generalized Newton methods. The obtained results
demonstrate a high performance of the algorithms and their applications to
solving the $\ell_{0}$-$\ell_{2}$ regularized least square regression model,
with providing numerical experiments.

## Funding and Conflicts of Interests

Research of Boris S. Mordukhovich was partly supported by the US National
Science Foundation under grant DMS-2204519, by the Australian Research Council
under Discovery Project DP190100555, and by the Project 111 of China under
grant D21024. Research of Vo Thanh Phat was partly supported by the US
National Science Foundation under grant DMS-2204519. The authors declare that
the presented results are new, and there is no any conflict of interest.

## References

  * [1] Aragón-Artacho, F. J., Mordukhovich, B.S., Peréz-Aros, P.: Coderivative-based semi-Newton method in nonsmooth difference programming. arXiv:2301.03491 (2023)
  * [2] Attouch, H., Bolte, J., Svaiter, B.F.: Convergence of descent methods for semi-algebraic and tame problems: proximal algorithms, forward–backward splitting, and regularized Gauss–Seidel methods. Math. Program. 137, 91–129 (2013)
  * [3] Bauschke, H.H., Combettes, P.L.: Convex Analysis and Monotone Operator Theory in Hilbert Spaces, 2nd edition. Springer, New York (2017)
  * [4] Chieu, N.H., Lee, G.M., Yen, N.D.: Second-order subdifferentials and optimality conditions for ${\cal C}^{1}$-smooth optimization problems. Appl. Anal. Optim. 1, 461–476 (2017)
  * [5] Dontchev, A.L., Rockafellar, R.T.: Implicit Functions and Solution Mappings: A View from Variational Analysis, 2nd edition. Springer, New York (2014)
  * [6] Drusvyatskiy, D., Lewis, A.S.: Tilt stability, uniform quadratic growth, and strong metric regularity of the subdifferential. SIAM J. Optim. 23, 256–267 (2013)
  * [7] Facchinei, F., Pang, J.-C.: Finite-Dimensional Variational Inequalities and Complementarity Problems, Vols. I and II. Springer, New York (2003)
  * [8] Gfrerer, H.: On directional metric regularity, subregularity and optimality conditions for nonsmooth mathematical programs. Set-Valued Var. Anal. 21, 151–176 (2013)
  * [9] Gfrerer, H., Outrata, J.V.: On a semismooth∗ Newton method for solving generalized equations. SIAM J. Optim. 31, 489–517 (2021)
  * [10] Gfrerer, H., Outrata, J.V.: On (local) analysis of multifunctions via subspaces contained in graphs of generalized derivatives. J. Math. Anal. Appl. 508, 125–895 (2022)
  * [11] Ginchev, I., Mordukhovich, B.S.: On directionally dependent subdifferentials. C. R. Acad. Bulg. Sci. 64, 497–508 (2011)
  * [12] Hazimeh, H., Mazumder, R.: Fast best subset selection: coordinate descent and local combinatorial optimization algorithms. Operations Research, 68, 1517–1537 (2020)
  * [13] Izmailov, A.F., Solodov, M.V.: Newton-Type Methods for Optimization and Variational Problems. Springer, New York (2014)
  * [14] Kha, P.T., Khanh, P.D.: Variational inequalities governed by strongly pseudomonotone operators, Optimization, 71, 1983-2004 (2021)
  * [15] Khanh, P.D., Khoa, V. V. H., Mordukhovich, B.S., Phat, V.T.: Second-order subdifferential optimality conditions in nonsmooth optimization. arXiv:2312.1627 (2023)
  * [16] Khanh, P.D., Khoa, V. V. H., Mordukhovich, B.S., Phat, V.T.: Variational convexity of functions in Banach spaces, in Functional Analysis and Continuous Optimization, J. M. Amigó, M. J. Cánovas et al. (eds.). Springer Proceedings in Mathematics and Statistics 424, pp. 237–260, Springer, Cham, Switzerland (2023)
  * [17] Khanh, P.D., Khoa, V. V. H., Mordukhovich, B.S., Phat, V.T.: Variational and strong variational convexity in infinite-dimensional variational analysis. arXiv:2309.13779 (2023)
  * [18] Khanh, P.D., Mordukhovich, B.S., Phat, V.T.: A generalized Newton method for subgradient systems. Math. Oper. Res. 48, 1811–1845 (2023)
  * [19] Khanh, P.D., Mordukhovich, B.S., Phat, V.T.: Variational convexity of functions and variational sufficiency in optimization. SIAM J. Optim. 33, 1121–1158 (2023)
  * [20] Khanh, P.D., Mordukhovich, B.S., Phat, V.T., Tran, D. B.: Generalized damped Newton algorithms in nonsmooth optimization via second-order subdifferentials. J. Global Optim. 86, 93–122 (2023)
  * [21] Khanh, P.D., Mordukhovich, B.S., Phat, V.T., Tran, D. B.: Globally convergent coderivative-based generalized Newton methods in nonsmooth optimization. Math. Program. DOI: 10.1007/s10107-023-01980-2 (2023)
  * [22] Klatte, D., Kummer, B.: Nonsmooth Equations in Optimization: Regularity, Calculus, and Application., Kluwer, Boston (2002)
  * [23] Mordukhovich, B.S.: Sensitivity analysis in nonsmooth optimization. In: Field, D.A., Komkov, V.(eds) Theoretical Aspects of Industrial Design, pp. 32–46. SIAM Proc. Appl. Math. 58. Philadelphia, PA (1992)
  * [24] Mordukhovich, B.S.: Complete characterizations of openness, metric regularity, and Lipschitzian properties of multifunctions. Trans. Amer. Math. Soc. 340, 1–35 (1993)
  * [25] Mordukhovich, B.S.: Variational Analysis and Generalized Differentiation, I: Basic Theory, II: Applications. Springer, Berlin (2006)
  * [26] Mordukhovich, B.S.: Variational Analysis and Applications. Springer, Cham, Switzerland (2018)
  * [27] Mordukhovich, B.S.: Second-Order Variational Analysis in Optimization, Variational Stability, and Control: Theory, Algorithms, Applications. Springer, Cham, Switzerland (2024)
  * [28] Mordukhovich, B.S., Outrata, J.V.: On second-order subdifferentials and their applications. SIAM J. Optim. 12, 139–169 (2001)
  * [29] Mordukhovich, B.S., Sarabi, M.E.: Generalized Newton algorithms for tilt-stable minimizers in nonsmooth optimization. SIAM J. Optim. 31, 1184–1214 (2021)
  * [30] Patrinos, P., Bemporad, A.: Proximal Newton methods for convex composite optimization. In: IEEE Conference on Decision and Control, 2358–2363 (2013)
  * [31] Poliquin, R.A., Rockafellar, R.T.: Prox-regular functions in variational analysis. Trans. Amer. Math. Soc. 348, 1805–1838 (1996)
  * [32] Poliquin, R.A., Rockafellar, R.T.: Tilt stability of a local minimum. SIAM J. Optim. 8, 287–299 (1998)
  * [33] Stella, L., Themelis, A., Patrinos, P.: Forward–backward quasi-Newton methods for nonsmooth optimization problems. Comput. Optim. Appl. 67, 443–487 (2017)
  * [34] Stella, L., Themelis, A., Patrinos, P.: Forward-backward envelope for the sum of two nonconvex functions: further properties and nonmonotone linesearch algorithms. SIAM J. Optim. 28, 2274–2303 (2018)
  * [35] Rockafellar, R.T., Wets R.J-B.: Variational Analysis. Springer, Berlin (1998)
  * [36] Rockafellar R.T.: Variational convexity and local monotonicity of subgradient mappings. Vietnam Journal of Mathematics, Vietnam J. Math. 47, 547–561 (2019)
  * [37] Rockafellar, R.T.: Augmented Lagrangians and hidden convexity in sufficient conditions for local optimality. Math. Program. 198, 159–194 (2023)