# Influence of surfactant concentration on drop production by bubble bursting

Juliette Pierre<EMAIL_ADDRESS>Mathis Poujol  Thomas
Séon  Sorbonne Université, CNRS, UMR 7190, Institut Jean Le Rond
$\partial$’Alembert, F-75005 Paris, France

###### Abstract

Bubble bursting at the surface of the sea water produce drops and is the main
source of sea spay aerosol. The mechanisms underlying the drops production
from a single bubble bursting event have been intensively studied and the
influence of the bubble size and liquid parameters (density, viscosity and
surface tension) has been unified. However, despite the diversity of the
surfactant molecules present in the oceans, their influence has been
overlooked. In this paper we experimentally explore the influence of the
surfactant concentration (Sodium Dodecyl Surfate - SDS) in a water solution on
a single bubble collapse and subsequent drop production. We show that these
surfactant molecules have an astonishing effect. In particular, we
quantitatively show that they modify the bubble collapse, they induce less,
smaller and faster drops and they can even completely prevent the drop
production for a particular concentration. These results are presented in the
existing dimensionless framework and they allow us to affirme that these
effects are mainly a consequence of the surface tension gradients (Marangoni
stresses) more than just the surface tension lowering. Finally, this study
shows that the role of water contamination by surface-active agents is
important and needs to be understood to improve the prediction of the sea
spray aerosol.

###### pacs:

Valid PACS appear here

##  I Introduction

As a wave crashes in the ocean it entrains air below the surface. After a
turbulent break-up cascade [1] a population of bubbles is produced [2, 3, 4]
and while small bubbles may be dissolved into the water, larger bubbles rise
back to the surface and collapse [5]. The bursting of a bubble starts with the
break up of the thin liquid film that separates the air cavity from the
atmosphere and ends up with the fragmentation of a rising jet. Through these
two fragmentation events, bubble bursting produces film drops [6, 7] and jet
drops [8, 9, 10, 11, 12, 13] and constitutes one of the main sources of the
ocean spray [14]. By evaporating, the sea spray transports in the atmosphere
water vapor, important for the thermodynamics of the atmosphere, and salt
crystals that affect the radiative balance of the atmosphere and form cloud
condensation nuclei [15, 16]. And, in a significant way, they carry also heat,
dissolved gases, surfactants, biological materials [5]. Finally, uncertainties
in predicting sea spray aerosols characteristics directly impacts our ability
to perform weather prediction and earth system modeling [17, 18].

Since the pioneer work of D. Blanchard [19] there have been a number of -
experimental, numerical and theoretical - combining studies on a single bubble
bursting, that brought comprehensive data on the size and speed of the jet
drops produced by bubble bursting in water [20, 21, 22, 23]. Applying these
results to the bubble size distribution produced under a breaking wave enabled
a rough estimation of the statistics of jet drop production [24]. However, the
ocean surface is partly covered by a biofilm, which can be modeled with
surfactants [25]. The surface-active contaminations are known to modify the
static and dynamic behaviors of bubbles, including their coalescence,
lifetimes, and bursting [26, 27, 28]. Consequently, the influence of the
physicochemistry of the interface has to be taken into account in the process
of bubble collapse at the interface and in the subsequent drop production.
Néel & Deike (2021) [28] considered a monodisperse assembly of millimetric air
bubbles produced identically in the bulk for a wide range of surface
contamination and showed that, depending on the contamination, the bubble
distribution that bursts can be very distinct from the initial distribution.
There have been various experiments that attempted to described the role of
the physicochemical parameters on the production of droplets by bursting
bubbles [29, 3, 30, 5], but there are large variations in protocols, and the
influence of surfactants on the drop production remains largely unclear. All
these experiments are realized on a large collection of bubbles, with
different size distributions, suggesting the need to carry out a study on a
single bursting bubble. Recently, Constante-Amores et al. (2021) [31] studied
the effect of surfactant on the dynamics of a bubble bursting, using numerical
simulations, accounting for sorption kinetics and diffusive effects. At one
fixed bubble size and one surface contamination, they showed that the presence
of surfactant affects the dynamics of the system through Marangoni-induced
flow and is responsible for delaying the collapse and generating slower and
fewer drops.

In this article, we study experimentally the effect of Sodium Dodecyl Sulfate
(SDS) surfactant on the dynamics of a bubble bursting through an interface.
After describing the experimental setup, we show qualitatively that the
surfactants have an astonishing influence on the jet dynamics subsequent to
the bubble collapse, and on the jet drops production. The following is
dedicated to quantify this effect by varying the surfactant concentration and
the bubble size. We start by studying the influence of the surfactants on the
bubble collapse time, before characterizing the variation of the number, size
and speed of the ejected drops as a function of the control parameters.
Finally, we focus on the influence of the surfactant concentration on the
cavity collapse and the capillary waves dynamics.

Figure 1: Sequences of bursting bubbles (bubbles start to burst at 0 ms) with
comparable radius in three different liquids : (a) water, $\gamma$ = 70
mN.m-1, R = 830 $\mu$m ; (b) water-ethanol solution with respectively 89.5%
and 10.5% of total weight, $\gamma$ = 48 mN.m-1, R = 830 $\mu$m ; (c) water-
SDS solution with 3.4 mM of SDS, $\gamma$ = 40 mN.m-1, R = 840 $\mu$m. The
surfactants have a strong influence on the jet dynamics, independently of
their influence on the surface tension.

##  II Experimental setup

The experiment consists in releasing a single air bubble from a submerged
needle in a liquid and recording the upward jet and released drops after the
bubble bursts at the free surface. Air bubbles are generated in a
parallelipedal glass tank (20 cm length, 14 cm width, 9.5 cm depth) filled
with either tap water or an aqueous solution of SDS (Sodium Dodecyl sulfate -
purchased from Sigma Aldrich) surfactant with a mass concentration ranging
from 0.5 g/L to 10 g/L, i.e. 1.7 mM to 34.7 mM. For SDS at the ambient
temperature the critical micelle concentration (CMC) is found to be around 8mM
[32], which means that the SDS concentration in our solutions varies from C=
0.2 CMC to 4.3 CMC. Bubbles are generated using a syringe pump filled with
air. Three different needles are used, with internal diameter varying from
0.08 mm to 1.5 mm enabling to create bubble with thee different radii : 0.8,
1.1, 1.7 $\pm$ 0.1 mm. The bubbles rise to the surface and briefly float
before bursting. Considering the elliptic shape of the floating bubble, we
defined an equivalent radius as $R=(a^{2}b)^{1/3}$ with $a$ and $b$
respectively the semi-major and semi-minor axes of the ellipse. The surface
tension of each solution is measured using the pendant drop technique [33].

In all experiments a digital high speed camera (Phantom V2511) is used to
image the rising jet and releasing drops from the side, above the free
surface. In a few experiments a second digital high speed camera (Photron
SA-5) is added to image the collapse of the submerged cavity below the free
surface.

##  III Qualitative description

Figure 1 presents three sequences of bubble busting. In each case the bubble
radius is almost constant and the liquid is different : (a) tap water
($\gamma$ = 70 mN.m-1), (b) a water-ethanol mixture with a surface tension
$\gamma$ = 48 mN.m-1, and (c) a water-SDS solution with a surface tension
$\gamma$ = 40 mN.m-1. We observe on (a) and (b) that the decrease of surface
tension does not affect much the drop size and velocity. These observations
have been reported quantitatively in the literature [21, 9, 20, 11, 22]. In
the sequence (c), the bubble bursts in a liquid with surfactants concentrated
at 0.4 CMC and with a surface tension very close to that of sequence (b). The
result is remarkable. The presence of surfactants completely changes the jet
dynamics, the jet velocity is so low that it can barely reach the free surface
and cannot produce any droplet. In the following, our goal is to examine more
quantitatively the influence of the surfactant concentration on the drop
dynamics, and to look at where and how, in the cavity collapse process, the
surfactants can have such a strong influence.

##  IV The bubble collapse

Figure 2:  (a) Bubble collapsing time $\Delta t$ as a function of the
surfactant concentration, adimensionalized using the critical micelle
concentration (CMC). The CMC is taken equal to 8 mM [32]. Inset: for each
solution the surface tension is measured using the pendant drop technique and
reported as a function of the surfactant concentration. (b) Collapsing time,
$\Delta t$, normalized by the capillaro-inertial time scale, $\rho
R^{3}/\gamma$, as a function of the dimensionless surfactant concentration,
for three different bubble radius R : 0.8, 1.1 and 1.7 mm. The 0.3 prefactor
is added in the normalization of $\Delta t$ in order to have the dimensionless
collapsing time equal to one. The color of the markers is associated to the
surfactant concentration in the liquid.

Before the cavity collapses, the bubble is floating at the free surface. As it
is static, its shape is due to an equilibrium between capillarity and gravity
and is obtained by integration of the Young–Laplace equation [34, 7, 35].
Surfactants have no more influence on the static bubble shape than through
their modification of the surface tension. The static bubble shape before
bursting cannot be responsible for the modification of the jet dynamics
between figure 1(b) and (c).

We then focus on the influence of the surfactant concentration on the bubble
collapse duration $\Delta t$. $\Delta t$ is defined as the time elapsed
between the hole nucleation in the cap film and the cavity reversal, when the
depth of the immersed cavity starts to decrease. Figure 2 (a) presents $\Delta
t$ as a function of C(CMC) the surfactant concentration adimensionalized using
the CMC. Note that the color of the markers is associated to the surfactant
concentration in the liquid. The different markers at one concentration show
$\Delta t$ for experiments in similar conditions, they therefore only reflect
the dispersion of the results. We observe that, independently of the
dispersion, the bubble collapsing time increases with the surfactant
concentration, reaches a maximum close to the CMC and decreases. In other
words, up to the CMC, the cavity is slower to collapse as the surfactants are
more concentrated and, above the CMC, it becomes faster again. This non-
monotonic variation of the collapsing time with the surfactant concentration
is surprising. In particular because the variation of the surface tension
$\gamma$ with the dimensionless surfactant concentration $C$ is monotonic, as
verified in the inset that presents the variation of $\gamma$ with $C$. In
this inset $\gamma$ expectedly decreases with $C$, from the water surface
tension, until it reaches a plateau, around 36 mN/m, beyond the CMC.
Consequently, the non-monotonic variation of $\Delta t$ with $C$ indicates
that the surfactant dynamics should play a role in the cavity collapse.

In this capillaro-inertial collapse, $\Delta t$ is expected to scale as the
capillaro-inertial time : $\rho R^{3}/\gamma$, with $\rho$ the liquid density
and $R$ the bubble radius [35]. Figure 2 (b) presents the collapsing time
$\Delta t$ normalized by this capillaro-inertial time scale as a function of
the dimensionless surfactant concentration $C$, for the three different bubble
radii $R=$ 0.8, 1.1 and 1.7 mm. As expected, without surfactant ($C=0$), all
the bubbles with different radius collapse, demonstrating that this
adimensionalized collapsing time is relevant. The prefactor 0.3 is added to
the capillaro-inertial time so that the normalized times collapse around 1.
However, as the surfactants are added in solution, the dimensionless
collapsing time increases and become more dispersed with increasing bubble
size. The dimensionless time and its dispersion both reach a maximum between
half the CMC and the CMC (C= 0.5 - 1) and then decrease. As we approach the
CMC, the experiments with different bubble radii do not collapse anymore and
the small bubbles are relatively slower to collapse compared to the larger
ones.

By normalizing $\Delta t$ by the relevant capillaro-inertial collapsing time,
figure 2 (b) enables to decorrelate the respective influence of the surface
tension and the surfactant dynamics. As the data do not rescale in the
presence of the surfactants, they are expected to be a consequence of the
particular dynamics of the surfactants, independently of their influence on
the measured liquid surface tension $\gamma$ displayed in the inset. Indeed,
in processes that involve surface stretching and/or capillary waves, as it is
the case in our experiment, gradients of surfactants can appear, generating
Marangoni stresses that affect the dynamics [36, 37]. Constante-Amores et al.
[31] showed that in the insoluble surfactant limit, the collapse yields to an
over-concentration of surfactants at the apex of the cavity when the capillary
waves focus, source of a strong Marangoni stress that can delay the cavity
collapse. The presence of insoluble surfactants is also known to influence the
dispersion of surface waves and to enhance capillary wave damping due to the
interfacial rigidification [38, 39].

In the aim of contextualizing our results within the existing literature, we
need to discuss the surfactant dynamics in our experiment. Two kinetics must
be considered: the surface diffusivity and the surfactant rearrangement
between the surface and the bulk. First, the surface diffusivity $D_{s}$ of
SDS is around 10-9 m2.s-1 [40]. Thus, the Peclet number
$\text{Pe}=\sqrt{\gamma R/\rho}/D_{s}$ that measures the relative importance
of surface convection of surfactant to its diffusion is $\mathcal{O}(10^{4})$
in our experiment. The surfactant surface diffusion is therefore not strong
enough to mitigate its advection.

Second, to estimate the adsorption-desorption dynamics of surfactants one can
compare the characteristics time of the sorption rates to the diffusion time
of the surfactants from the liquid-gas interface to the bulk. The Langmuir
model gives the characteristics time for the sorption kinetics :
$\tau_{b}=\Gamma_{m}(k_{\text{des}}\Gamma_{m}+k_{\text{ads}}C)$ with
$\Gamma_{m}=4.10^{-6}$ mol.m-2 [41] the maximum surface packing concentration
of SDS at the air-water interface and $k_{\text{des}}$ and $k_{\text{ads}}$
respectively the desorption and adsorption rate [42]. The typical diffusive
time-scale can be express as $t_{\text{diff}}=\Gamma^{2}/(D_{v}C^{2})$ with
$\Gamma$ the surface concentration at equilibrium [43] and $D_{v}=5.10^{-10}$
m2.s-1 [44] the diffusion coefficient of the surfactants in the bulk liquid.
Even if the characteristic time $\tau_{b}$ is not easy to estimate, in
particular due to the lack of precision in the sorption rates $k_{\text{des}}$
and $k_{\text{ads}}$, the diffusion time, of the order of the millisecond,
seems to remain larger than the time of the adsorption-desorption kinetics.
This indicates that the dynamics of surfactants is limited by the diffusion,
which is of the same order of magnitude than the time of collapse $\Delta t$.
Thus, we expect to have a comparable time between surfactant dynamics and
bubble collapse and, subsequently, to have solubility effects during the
collapse. However, as the surprising variation of the collapsing time $\Delta
t$ with the surfactant concentration $C$ is interpreted as a consequence of
the variation of the local concentration of the surfactants, we expect the
solubility not to be dominant. Consequently, the characteristics times of the
surfactant dynamics might be underestimated, or the dynamics might be
influenced by the presence of dodecanol, an insoluble impurity which may be
either present from the original preparation or produced by hydrolysis of the
aqueous SDS on standing [41]. The use of numerical simulation would now be an
asset in the aim of interpreting this complex situation [36, 31].

##  V Jet drops characteristics

Figure 3: (a) Number of ejected drops as a function of the SDS concentration.
(b) Capillary number of the first ejected drop as a function of La(1+2.2Bo)
(c) Laplace number of the first ejected drop as function of the Laplace number
of initial bubble. In all graphs color codes for the SDS concentration and
symbol codes for the bubble size. In (b) and (c) the dashed lines represent
the unified relations from [45]. The gray shadows traduce the dispersion of
the experimental and numerical values all reported in [22].

As shown in figure 1, after a bubble has collapsed, an upward jet usually
rises and produces the so-called jet drops. Figure 3 (a) presents the number
of these drops produced when a bubble bursts, as a function of the surfactant
concentration, for 3 different bubble radii. For bubbles bursting in water
(C=0 CMC, empty markers) the 3 different bubbles produce droplets and the
smaller the bubble, the more drops are produced [22]. When surfactants are
added to the water, our largest bubble size (triangle) cannot produce drops
anymore, regardless of the surfactant concentration tested. This is a strong
result, for this bubble size, surface contamination completely kills the drop
production. Then, for smaller bubble, we observe that there is dispersion in
the number of produced droplets. Indeed, for C $\simeq$ 0.2 CMC (yellow
markers), R = 0.8 and 1.1 mm are superimposed and produce either 1, 2 ou 3
droplets. There is the same kind of dispersion for C $\geq$ 0.9 CMC.
Nevertheless, despite this dispersion, the trend is clear: (i) surface
contamination can prevent the drop production and (ii) when droplets are
produced there are less numerous and their number seems to decrease down to a
minimum around half the CMC, before increasing again. These results are
crucial, they signify that the size distribution of ejected jet drops produced
in pure water [24] might be very different than the one produced in water with
surfactant. They precise and experimentally validate the recent numerical
results of Constante-Amores et al. [31] that show that a reduction in the
number of ejected droplets arises with surfactant-laden flow due to Marangoni
flow.

To go further in influence of the surfactant contamination in the jet drop
production, the speed and size of the first ejected droplet are quantified as
a function of the surfactant concentration. Based on a large amount of
numerical and experimental results, previous studies have demonstrated that
the problem has two control parameters : the main one, the Laplace number
(La), which compares the capillaro-inertial forces with the viscous forces,
and the Bond number (Bo), which compares the gravitational forces with the
capillary ones [21, 12, 45, 22]. They are defined as:

La $\displaystyle=\frac{\rho\gamma R}{\mu^{2}}$ (1) Bo
$\displaystyle=\frac{\rho gR^{2}}{\gamma}$ (2)

where $R$ is the bubble radius, $\mu$ the water viscosity, $\rho$ the water
density, $\gamma$ the surface tension and $g$ the acceleration of gravity. The
first drop speed $V_{d}$ and size $R_{d}$ are also adimensionalized using,
respectively, the visco-capillary velocity $V_{\mu}=\mu/\gamma$ and length
$l_{\mu}=\mu^{2}/(\rho\gamma)$, yielding the dimensionless drop speed and
size:

$\displaystyle\text{Ca}_{d}$ $\displaystyle=\frac{V_{d}\mu}{\gamma}$ (3)
$\displaystyle\text{La}_{d}$ $\displaystyle=\frac{\rho\gamma R_{d}}{\mu^{2}}$
(4)

Within this dimensionless framework, previous studies [45, 22] have proposed
universal rescalings able to fully describe the first drop velocity and size.
These scalings are respectively represented with dashed line on figures 3 (b)
and (c). They gather a large range of bubble size and liquid parameters
($\rho$, $\mu$, $\gamma$). The grey zone around the dashed line represents the
error bar including all the experiences. The Bond number appears in the x-axis
as a correction term for the drop velocity (Cad) and plays no role for the
drop size (Lad).

On these plots we add here the values measured with bubbles bursting in our
solutions of SDS mixed to water. The surfactant concentration is represented
using the same colors as in figures 2 and 3(a). First, expectedly, the drop
velocity and size from bubble busting in water, with empty markers, fall onto
the universal scalings, and many markers are on the x-axis because no drop are
produced for the largest bubble and for the concentration close to half of the
CMC, as shown in figure 3(a). Secondly, in figure 3(b), we observe that the
more concentrated is the solution, the higher the drop velocity is above the
scaling. In figure 3(c), even with a small amount of surfactant, the drop size
falls quite far below the universal scaling. It also seems that the drop size
is less affected by the Laplace number as the surfactant concentration is
high. These variations indicate that the influence of the surfactants is
highly non trivial, undoubtedly dependent on the local gradient concentration
along the jet. Finally, we observe again that the dispersion is quite larger
than without surfactant, may be because the dynamics is very sensitive to the
balance between the coupled dynamics of surfactants and the jet. Consequently,
the next step of this study will probably need a statistical characterization
to properly capture the influence of the surfactant concentration in the drop
production [46].

##  VI Capillary waves focusing

Figure 4: Sequences showing the collapse of the submerged cavity for bubbles
with almost the same radius in three different liquids: (a) water,
$\gamma=70\,$mN.m-1, $R=1.1\,$mm, (b) water-SDS solution at 0.4 CMC,
$R=1.1\,$mm, (c) water-SDS solution at 4.3 CMC, $R=1.0\,$mm.

The jet dynamics strongly depends on the capillary waves focusing at the
bottom of the cavity [9, 47]. Figure 4 presents 3 sequences of cavity collapse
with almost the same bubble radius and 3 different surfactant concentrations:
(a) no surfactant, (b) 0.4 CMC and (c) 4.3 CMC. The capillary waves
propagation is different in these 3 sequences. The clearest difference appears
between (a) and (b), and lies in the wave shape, in particular in the second
half of the collapse sequences. The shape of the lowest collapsing cavity
shown in the second to last images is completely different, and undoubtedly
explains the strong difference in the drop production dynamics (see $C=0.4$
CMC in figure 3). As shown by Constante-Amores et al. [31] for precise values
of the control parameters, the interfacial surfactant concentration reaches
its maximum value as the surfactant-laden capillary waves converge on the
cavity apex. The Marangoni stresses that drive motion from high to low surface
concentration regions can explain the shape of the cavity.

For the highest concentration (c) the shape of the capillary waves, and of the
lowest cavity, looks quite similar to the one in water. This seems to indicate
that for concentration higher than the CMC, the Marangoni stresses are lower,
due to lower concentration gradient, and this can be explained by a smaller
diffusive time of the surfactants for high concentration, giving a surfactant
soluble behavior to the collapse. This similarity in the cavity collapse
between no contamination and high contamination probably explains why the drop
characteristics are closer between water and the highest concentration than
water and 0.4 CMC, for which there is no drop produced.

##  VII Conclusion

Surfactants are most often present in the liquids where bubbles are bursting
(ocean, sparkle wine, soda…) and have barely been taken into account in the
experiments and the models. In particular, to our knowledge, experiments of a
single bubble bursting in a surfactant laden liquid have never been carried
out. Here, we have shown how the SDS strongly influences the cavity collapse
and the drops production, for different values of the bubble size and of the
surfactant concentration. In particular, we highlight that the contamination
induces: (i) a maximum in the bubble collapse duration around the CMC, (ii)
smaller and faster drops and (iii) less drops, with no drop at all for a
particular concentration of half the CMC. We also show that these effects are
a consequence of the surface tension gradients (Marangoni stresses) and not
just the surface tension lowering. The exact role of the Marangoni flows is
not known and needs to be clarified by quantifying the surface tension
gradients appearing during the bubble collapse.

In the following, motivated by this study, more experiments should be done, in
particular with insoluble surfactants to examine the influence of the
solubility. As it is complex to only change one parameter at a time, it would
undoubtedly be interesting to carry out a large campaign of numerical
simulations. Indeed, the cavity collapse, cavity reversal, jet dynamics and
end pinching are very complex phenomenon and their dynamics involve a high
dependency on the local gradient concentration. In simulations Marangoni
stresses can be turned off while surfactant-induced lowering of surface
tension can be retained, thereby determining which of the two effects is the
dominant mechanism by which surfactants affect the flow [36]. On the other
hand, simulations would also enable a statistical characterization of the drop
production [46] as it seems to be very sensitive to the experimental
conditions and probably to the initial conditions. Finally, in any way, the
surfactant concentration needs to be taken into account in the experiments and
simulations to improve the prediction of a real bubble bursting spray as the
sea spray.

## References

  * [1] Chris Garrett, Ming Li, and David Farmer.  The connection between bubble size spectra and energy dissipation rates in the upper ocean.  Journal of physical oceanography, 30(9):2163–2171, 2000.
  * [2] Grant B Deane and M Dale Stokes.  Scale dependence of bubble creation mechanisms in breaking waves.  Nature, 418(6900):839–844, 2002.
  * [3] Kimberly A Prather, Timothy H Bertram, Vicki H Grassian, Grant B Deane, M Dale Stokes, Paul J DeMott, Lihini I Aluwihare, Brian P Palenik, Farooq Azam, John H Seinfeld, et al.  Bringing the ocean into the laboratory to probe the chemical complexity of sea spray aerosol.  Proceedings of the National Academy of Sciences, 110(19):7550–7555, 2013.
  * [4] Luc Deike, W Kendall Melville, and Stéphane Popinet.  Air entrainment and bubble statistics in breaking waves.  Journal of Fluid Mechanics, 801:91–129, 2016.
  * [5] Luc Deike.  Mass transfer at the ocean–atmosphere interface: The role of wave breaking, droplets, and bubbles.  Annual Review of Fluid Mechanics, 54, 2021.
  * [6] D. C. Blanchard and L. D. Syzdek.  Film drop production as a function of bubble size.  J. Geophys. Res., 93(C4):3649–3654, 1988.
  * [7] H. Lhuissier and E. Villermaux.  Bursting bubble aerosols.  Journal of Fluid Mechanics, 696:5–44, 2011.
  * [8] D. E. Spiel.  More on the births of jet drops from bubbles bursting on seawater surfaces.  J. Geophys. Res., 102(C3):5815–5821, 1997.
  * [9] Elisabeth Ghabache, Arnaud Antkowiak, Christophe Josserand, and Thomas Séon.  On the physics of fizziness: How bubble bursting controls droplets ejection.  Physics of Fluids (1994-present), 26(12):–, 2014.
  * [10] E. Ghabache and T. Séon.  Size of the top jet drop produced by bubble bursting.  Phys. Rev. Fluids, 1(051901), 2016.
  * [11] C Frederik Brasz, Casey T Bartlett, Peter LL Walls, Elena G Flynn, Yingxian Estella Yu, and James C Bird.  Minimum size for the top jet drop from a bursting bubble.  Physical Review Fluids, 3(7):074001, 2018.
  * [12] Alfonso M Gañán-Calvo.  Revision of bubble bursting: universal scaling laws of top jet drop size and speed.  Physical review letters, 119(20):204502, 2017.
  * [13] Francisco J Blanco-Rodríguez and JM Gordillo.  On the sea spray aerosol originated from bubble bursting jets.  Journal of Fluid Mechanics, 886, 2020.
  * [14] F Veron, C Hopkins, EL Harrison, and JA Mueller.  Sea spray spume droplet production in high wind speeds.  Geophysical Research Letters, 39(16), 2012.
  * [15] E. R. Lewis and S. E. Schwartz.  Sea Salt Aerosol Production. Mechanisms, Methods, Measurements, and Models.  American Geophysical Union, Washington, DC,, geophysical monograph 152\. edition, 2004.
  * [16] Fabrice Veron.  Ocean spray.  Annual Review of Fluid Mechanics, 47(1):507–538, 2015.
  * [17] Gerrit de Leeuw, Edgar L Andreas, Magdalena D. Anguelova, C. W. Fairall, Ernie R. Lewis, Colin O’Dowd, Michael Schulz, and Stephen E. Schwartz.  Production flux of sea spray aerosol.  Rev. Geophys., 49(2), 05 2011.
  * [18] Luc Deike and W Kendall Melville.  Gas transfer by breaking waves.  Geophysical Research Letters, 45(19):10–482, 2018.
  * [19] D. C Blanchard.  The electrification of the atmosphere by particles from bubbles in the sea.  Progress In Oceanography, 1:73 – 112, IN7, 113–202, 1963.
  * [20] Thomas Séon and Gérard Liger-Belair.  Effervescence in champagne and sparkling wines: From bubble bursting to droplet evaporation.  The European Physical Journal Special Topics, 226(1):117–156, 2017\.
  * [21] L. Duchemin, S. Popinet, C. Josserand, and S. Zaleski.  Jet formation in bubbles bursting at a free surface.  Phys. Fluids, 14(9):3000–3008, 2002.
  * [22] Alexis Berny, Luc Deike, Thomas Séon, and Stéphane Popinet.  Role of all jet drops in mass transfer from bursting bubbles.  Physical Review Fluids, 5(3):033605, 2020.
  * [23] Alfonso M. Ganan-Calvo.  Scaling laws of top jet drop size and speed from bubble bursting including gravity and inviscid limit.  Phys. Rev. Fluids, 3:091601, Sep 2018.
  * [24] Alexis Berny, Stéphane Popinet, Thomas Séon, and Luc Deike.  Statistics of jet drop production.  Geophysical Research Letters, 48(10), 2021.
  * [25] Oliver Wurl, E Wurl, L Miller, K Johnson, and Svein Vagle.  Formation and global distribution of sea-surface microlayers.  Biogeosciences, 8(1):121–135, 2011.
  * [26] S Poulain, E Villermaux, and L Bourouiba.  Ageing and burst of surface bubbles.  Journal of fluid mechanics, 851:636–671, 2018.
  * [27] Daniel B Shaw and Luc Deike.  Surface bubble coalescence.  Journal of Fluid Mechanics, 915, 2021.
  * [28] B Néel and L Deike.  Collective bursting of free-surface bubbles, and the role of surface contamination.  Journal of Fluid Mechanics, 917, 2021.
  * [29] RL Modini, LM Russell, GB Deane, and MD Stokes.  Effect of soluble surfactant on bubble persistence and bubble-produced aerosol particles.  Journal of Geophysical Research: Atmospheres, 118(3):1388–1400, 2013.
  * [30] Patricia K Quinn, Douglas B Collins, Vicki H Grassian, Kimberly A Prather, and Timothy S Bates.  Chemistry and related properties of freshly emitted sea spray aerosol.  Chemical reviews, 115(10):4383–4399, 2015.
  * [31] Cristian R Constante-Amores, Lyes Kahouadji, Assen Batchvarov, Seungwon Shin, Jalel Chergui, Damir Juric, and Omar K Matar.  Dynamics of a surfactant-laden bubble bursting through an interface.  Journal of Fluid Mechanics, 911, 2021.
  * [32] V. Thominet, C. Stenvot, and D. Langevin.  Light scattering study of the viscoelasticity of soluble monolayers.  Journal of colloid and Interface Science, 126(1):54, 1987.
  * [33] Joseph D. Berry, Michael J. Neeson, Raymond R. Dagastine, Derek Y.C. Chan, and Rico F. Tabor.  Measurement of surface and interfacial tension using pendant drop tensiometry.  Journal of colloid and interface science, 454:226–237, 2015.
  * [34] Yoshiaki Toba.  Drop production by bursting of air bubbles on the sea surface (ii) theoretical study on the shape of floating bubbles.  Journal of the Oceanographical Society of Japan, 15(3):121–130, 1959.
  * [35] Mathis Poujol, Régis Wunenburger, François Ollivier, Arnaud Antkowiak, and Juliette Pierre.  Sound of effervescence.  Physical Review Fluids, 6(1):013604, 2021.
  * [36] Pritish M Kamat, Brayden W Wagoner, Sumeet S Thete, and Osman A Basaran.  Role of marangoni stress during breakup of surfactant-covered liquid threads: reduced rates of thinning and microthread cascades.  Physical Review Fluids, 3(4):043602, 2018.
  * [37] Harishankar Manikantan and Todd M. Squires.  Surfactant dynamics: hidden variables controlling fluid flows.  Journal of Fluid Mechanics, 892, apr 2020.
  * [38] J Lucassen and Robert S Hansen.  Damping of waves on monolayer-covered surfaces: I. systems with negligible surface dilational viscosity.  Journal of Colloid and Interface Science, 22(1):32–44, 1966.
  * [39] TJ Asaki, DB Thiessen, and PL Marston.  Effect of an insoluble surfactant on capillary oscillations of bubbles in water: observation of a maximum in the damping.  Physical review letters, 75(14):2686, 1995.
  * [40] Patrick T McGough and Osman A Basaran.  Repeated formation of fluid threads in breakup of a surfactant-covered jet.  Physical review letters, 96(5):054502, 2006.
  * [41] JR Lu, IP Purcell, EM Lee, EA Simister, RK Thomas, AR Rennie, and J Penfold.  The composition and structure of sodium dodecyl sulfate-dodecanol mixtures adsorbed at the air-water interface: a neutron reflection study.  Journal of Colloid and Interface science, 174(2):441–455, 1995\.
  * [42] Chien-Hsiang Chang and Elias I Franses.  Adsorption dynamics of surfactants at the air/water interface: a critical review of mathematical models, data, and mechanisms.  Colloids and Surfaces A: Physicochemical and Engineering Aspects, 100:1–45, 1995.
  * [43] Isabelle Cantat, Sylvie Cohen-Addad, Florence Elias, François Graner, Reinhard Höhler, Olivier Pitois, Florence Rouyer, and Arnaud Saint-Jalmes.  Foams: structure and dynamics.  OUP Oxford, 2013.
  * [44] Koji Kinoshita, Elisa Parra, and David Needham.  Adsorption of ionic surfactants at microscopic air-water interfaces using the micropipette interfacial area-expansion method: Measurement of the diffusion coefficient and renormalization of the mean ionic activity for sds.  Journal of colloid and interface science, 504:765–779, 2017.
  * [45] Luc Deike, Elisabeth Ghabache, Gérard Liger-Belair, Arup K Das, Stéphane Zaleski, Stéphane Popinet, and Thomas Séon.  Dynamics of jets produced by bursting bubbles.  Physical Review Fluids, 3(1):013603, 2018.
  * [46] Alexis Berny, Luc Deike, Stéphane Popinet, and Thomas Séon.  How size and speed of jet drops are robust to initial conditions.  Subm. to Physical Review Fluids, 2021.
  * [47] JM Gordillo and J Rodríguez-Rodríguez.  Capillary waves control the ejection of bubble bursting jets.  Journal of Fluid Mechanics, 867:556–571, 2019.


$U_{2}(t)=1-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,H_{2}(s)U_{2}(s),\qquad
U_{1}(t)^{\dagger}=1+\mathrm{i}\int_{0}^{t}\mathrm{d}s\,U_{1}(s)^{\dagger}H_{1}(s)^{\dagger},$
(A.14)

one gets (A.13) by integration by parts (A.9), for $A(t)=U_{1}(t)^{\dagger}$
with $a(t)=\mathrm{i}U_{1}(t)^{\dagger}H_{1}(t)^{\dagger}$, and
$B(t)=U_{2}(t)$ with $b(t)=-\mathrm{i}H_{2}(t)U_{2}(t)$. ∎

######  Corollary 5.

Let $H=H^{\dagger}\in L^{1}_{\mathrm{loc}}(\mathbb{R})$ be the bounded self-
adjoint generator of a continuous propagator $U$. Then, $U(t)$ is unitary for
all $t$,

$U(t)^{\dagger}U(t)=U(t)U(t)^{\dagger}=1.$ (A.15)

######  Corollary 6.

Let $H_{j}=H_{j}^{\dagger}\in L^{1}_{\mathrm{loc}}(\mathbb{R})$ be the bounded
self-adjoint generators of continuous propagators $U_{j}$, for $j=1,2$. Then,

$U_{2}(t)-U_{1}(t)=-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,U_{1}(t)U_{1}(s)^{\dagger}[H_{2}(s)-H_{1}(s)]U_{2}(s).$
(A.16)

Now, the next theorem provides an explicit derivation of the integration-by-
part lemma (Lemma 1) for locally integrable generators.

######  Theorem 9 (Integration-by-part lemma for locally integrable
generators).

Let $U_{j}$ be the continuous propagators generated by bounded self-adjoint
$H_{j}\in L^{1}_{\mathrm{loc}}(\mathbb{R})$, for $j=1,2$. Let

$S_{21}(t)=\int_{0}^{t}\mathrm{d}s\,[H_{2}(s)-H_{1}(s)].$ (A.17)

Then,

$U_{2}(t)-U_{1}(t)=-\mathrm{i}S_{21}(t)U_{2}(t)-\int_{0}^{t}\mathrm{d}s\,U_{1}(t)U_{1}(s)^{\dagger}[H_{1}(s)S_{21}(s)-S_{21}(s)H_{2}(s)]U_{2}(s).$
(A.18)

###### Proof.

By Corollary 6, $U_{2}(t)-U_{1}(t)=U_{1}(t)D(t)$, with

$D(t)=-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,U_{1}(s)^{\dagger}[H_{2}(s)-H_{1}(s)]U_{2}(s).$
(A.19)

Then, an application of (A.10) of Lemma 3 to $A(t)=U_{1}(t)^{\dagger}$ with
$a(t)=\mathrm{i}U_{1}(t)^{\dagger}H_{1}(t)$, $B(t)=S_{21}(t)$ with
$b(t)=H_{2}(t)-H_{1}(t)$, and $C(t)=U_{2}(t)$ with
$c(t)=-\mathrm{i}H_{2}(t)U_{2}(t)$ gives

$\displaystyle D(t)=-\mathrm{i}U_{1}(t)^{\dagger}S_{21}(t)U_{2}(t)$
$\displaystyle{}+\mathrm{i}\int_{0}^{t}\mathrm{d}s\,[\mathrm{i}U_{1}(s)^{\dagger}H_{1}(s)]S_{21}(s)U_{2}(s)$
$\displaystyle{}+\mathrm{i}\int_{0}^{t}\mathrm{d}s\,U_{1}(s)^{\dagger}S_{21}(s)[-\mathrm{i}H_{2}(s)U_{2}(s)],$
(A.20)

and the theorem is proved for $U_{2}(t)-U_{1}(t)=U_{1}(t)D(t)$. ∎

##  Appendix B Ergodic Means

In this appendix, we provide a few basic bounds on ergodic means, which are
used in the main text. We use the spectral norm.

######  Lemma 4.

Let $\\{P_{\ell}\\}$ be a set of $m$ Hermitian projections satisfying
$P_{\ell}=P_{\ell}^{\dagger}$ and $P_{k}P_{\ell}=\delta_{k\ell}P_{\ell}$ for
all $k$ and $\ell$. Then, for any bounded operators $\\{A_{\ell}\\}$, we have

$\biggl{\|}\sum_{\ell}A_{\ell}P_{\ell}\biggr{\|}\leq\sqrt{\sum_{\ell}\|A_{\ell}\|^{2}\vphantom{\biggr{\|}}}.$
(B.1)

Moreover, for any bounded operator $A$ and for any complex numbers
$\\{c_{k,\ell}\\}$, we have

$\biggl{\|}\sum_{k,\ell}c_{k,\ell}P_{k}AP_{\ell}\biggr{\|}\leq\sqrt{m}\max_{k,\ell}|c_{k,\ell}|\|A\|.$
(B.2)

###### Proof.

The spectral norm $\|A\|$ gives the largest singular value of an operator $A$,
and satisfies $\|A\|^{2}=\|A^{\dagger}A\|=\|AA^{\dagger}\|$. Thus,

$\displaystyle\biggl{\|}\sum_{\ell}A_{\ell}P_{\ell}\biggr{\|}^{2}$
$\displaystyle=\biggl{\|}\biggl{(}\sum_{k}A_{k}P_{k}\biggr{)}\biggl{(}\sum_{\ell}A_{\ell}P_{\ell}\biggr{)}^{\dagger}\biggr{\|}=\biggl{\|}\sum_{\ell}A_{\ell}P_{\ell}A_{\ell}^{\dagger}\biggr{\|}$
$\displaystyle\leq\sum_{\ell}\|A_{\ell}P_{\ell}A_{\ell}^{\dagger}\|=\sum_{\ell}\|A_{\ell}P_{\ell}\|^{2}\leq\sum_{\ell}\|A_{\ell}\|^{2},$
(B.3)

that is (B.1).

The second inequality (B.2) derives from (B.1),

$\displaystyle\biggl{\|}\sum_{k,\ell}c_{k,\ell}P_{k}AP_{\ell}\biggr{\|}^{2}$
$\displaystyle=\biggl{\|}\sum_{\ell}\,\biggl{(}\sum_{k}c_{k,\ell}P_{k}A\biggr{)}\,P_{\ell}\biggr{\|}^{2}\leq\sum_{\ell}\,\biggl{\|}\sum_{k}c_{k,\ell}P_{k}A\biggr{\|}^{2}$
$\displaystyle\leq\sum_{\ell}\,\biggl{\|}\sum_{k}c_{k,\ell}P_{k}\biggr{\|}^{2}\|A\|^{2}=\sum_{\ell}\max_{k}|c_{k,\ell}|^{2}\|A\|^{2}\leq
m\max_{k,\ell}|c_{k,\ell}|^{2}\|A\|^{2}.$ (B.4)

∎

In the following two lemmas, $\sum^{\prime}_{k,\ell}$ represents double
summations over $k$ and $\ell$ excluding terms with $k=\ell$.

######  Lemma 5 (Continuous ergodic mean).

Let $H$ be a bounded self-adjoint operator with the finite spectral
representation

$H=\sum_{\ell=1}^{m}E_{\ell}P_{\ell},$ (B.5)

where $\\{E_{\ell}\\}$ is the spectrum of $H$ and $\\{P_{\ell}\\}$ its
spectral projections. Let

$\eta=\mathop{\min_{k,\ell}}_{k\neq\ell}|E_{k}-E_{\ell}|$ (B.6)

be the minimal spectral gap of $H$. Then, for any bounded operator $A$ and for
any $t>0$, we have

$\biggl{\|}\frac{1}{t}\int_{0}^{t}\mathrm{d}s\,\mathrm{e}^{\mathrm{i}sH}A\mathrm{e}^{-\mathrm{i}sH}-\sum_{\ell}P_{\ell}AP_{\ell}\biggr{\|}\leq\frac{2\sqrt{m}}{\eta
t}\|A\|.$ (B.7)

###### Proof.

By using the spectral representation of $H$ in (B.5), we have

$\displaystyle\frac{1}{t}\int_{0}^{t}\mathrm{d}s\,\mathrm{e}^{\mathrm{i}sH}A\mathrm{e}^{-\mathrm{i}sH}-\sum_{\ell}P_{\ell}AP_{\ell}$
$\displaystyle={\sum_{k,\ell}}^{\prime}\frac{1}{t}\int_{0}^{t}\mathrm{d}s\,\mathrm{e}^{\mathrm{i}s(E_{k}-E_{\ell})}P_{k}AP_{\ell}$
$\displaystyle={\sum_{k,\ell}}^{\prime}\frac{\mathrm{e}^{\mathrm{i}t(E_{k}-E_{\ell})}-1}{\mathrm{i}t(E_{k}-E_{\ell})}P_{k}AP_{\ell}.$
(B.8)

By using Lemma 4, this is bounded by

$\displaystyle\biggl{\|}\frac{1}{t}\int_{0}^{t}\mathrm{d}s\,\mathrm{e}^{\mathrm{i}sH}A\mathrm{e}^{-\mathrm{i}sH}-\sum_{\ell}P_{\ell}AP_{\ell}\biggr{\|}$
$\displaystyle=\Biggl{\|}{\sum_{k,\ell}}^{\prime}\frac{\mathrm{e}^{\mathrm{i}t(E_{k}-E_{\ell})}-1}{\mathrm{i}t(E_{k}-E_{\ell})}P_{k}AP_{\ell}\Biggr{\|}$
$\displaystyle\leq\sqrt{m}\mathop{\max_{k,\ell}}_{k\neq\ell}\left|\frac{\sin[t(E_{k}-E_{\ell})/2]}{t(E_{k}-E_{\ell})/2}\right|\|A\|\leq\frac{2\sqrt{m}}{\eta
t}\|A\|.$ (B.9)

∎

######  Lemma 6 (Discrete ergodic mean).

Let $U$ be a unitary operator with the finite spectral representation

$U=\sum_{\ell=1}^{m}\mathrm{e}^{-\mathrm{i}\phi_{\ell}}P_{\ell},$ (B.10)

where $\\{\mathrm{e}^{-\mathrm{i}\phi_{\ell}}\\}$ is the spectrum of $U$ and
$\\{P_{\ell}\\}$ its spectral projections. Let

$\eta=\mathop{\min_{k,\ell}}_{k\neq\ell}|\mathrm{e}^{-\mathrm{i}\phi_{k}}-\mathrm{e}^{-\mathrm{i}\phi_{\ell}}|=2\mathop{\min_{k,\ell}}_{k\neq\ell}\left|\sin\frac{\phi_{k}-\phi_{\ell}}{2}\right|$
(B.11)

be the minimal spectral gap of $U$. Then, for any bounded operator $A$ and for
any integer $n\geq 1$, we have

$\Biggl{\|}\frac{1}{n}\sum_{j=0}^{n-1}{U^{\dagger}}^{j}AU^{j}-\sum_{\ell}P_{\ell}AP_{\ell}\Biggr{\|}\leq\frac{2\sqrt{m}}{\eta
n}\|A\|.$ (B.12)

###### Proof.

By using the spectral representation of $U$ in (B.10), we have

$\frac{1}{n}\sum_{j=0}^{n-1}{U^{\dagger}}^{j}AU^{j}-\sum_{\ell}P_{\ell}AP_{\ell}={\sum_{k,\ell}}^{\prime}\frac{1}{n}\sum_{j=0}^{n-1}\mathrm{e}^{\mathrm{i}j(\phi_{k}-\phi_{\ell})}P_{k}AP_{\ell}=\frac{1}{n}{\sum_{k,\ell}}^{\prime}\frac{1-\mathrm{e}^{\mathrm{i}n(\phi_{k}-\phi_{\ell})}}{1-\mathrm{e}^{\mathrm{i}(\phi_{k}-\phi_{\ell})}}P_{k}AP_{\ell}.$
(B.13)

By using Lemma 4, this is bounded by

$\displaystyle\Biggl{\|}\frac{1}{n}\sum_{j=0}^{n-1}{U^{\dagger}}^{j}AU^{j}-\sum_{\ell}P_{\ell}AP_{\ell}\Biggr{\|}$
$\displaystyle=\frac{1}{n}\,\Biggl{\|}{\sum_{k,\ell}}^{\prime}\frac{1-\mathrm{e}^{\mathrm{i}n(\phi_{k}-\phi_{\ell})}}{1-\mathrm{e}^{\mathrm{i}(\phi_{k}-\phi_{\ell})}}P_{k}AP_{\ell}\Biggr{\|}$
$\displaystyle\leq\frac{\sqrt{m}}{n}\mathop{\max_{k,\ell}}_{k\neq\ell}\left|\frac{\sin[n(\phi_{k}-\phi_{\ell})/2]}{\sin[(\phi_{k}-\phi_{\ell})/2]}\right|\|A\|\leq\frac{2\sqrt{m}}{\eta
n}\|A\|.$ (B.14)

∎

##  Appendix C Rotating Frames

In Remark 3, the freedom in the choice of a rotating frame is mentioned. Given
$U_{1}(t)$ and $U_{2}(t)$ unitary propagators, the “canonical” gauge for a
rotating frame is given by one of the two evolutions, say $U_{1}(t)$. In this
rotating frame, the distance between the two evolutions is given by the
distance between the relative evolution $U(t)=U_{1}(t)^{\dagger}U_{2}(t)$ and
the identity (zero Hamiltonian),

$\|U_{2}(t)-U_{1}(t)\|=\|U(t)-1\|.$ (C.1)

The Hamiltonian of the relative evolution $U(t)$ is given by

$H(t)=U_{1}(t)^{\dagger}[H_{2}(t)-H_{1}(t)]U_{1}(t),$ (C.2)

and the integral action is simply

$S(t)=\int_{0}^{t}\mathrm{d}s\,H(s).$ (C.3)

The Hamiltonian $H(t)$ may have an involute form and such canonical frame may
be unpractical. However, it has a conceptual merit as we are going to show.

One can show that the converse of (1.8) holds in the canonical frame, under
suitable assumptions on $H(t)$. Namely,

$U_{2}\approx U_{1}\quad\Rightarrow\quad S\approx 0.$ (C.4)

######  Proposition 3.

Let $t\mapsto U(t)$ be the unitary propagator generated by a locally
integrable time-dependent Hamiltonian $H(t)$, with $H(t)$ bounded and self-
adjoint for all $t\in\mathbb{R}$. Then, the integral action
$S(t)=\int_{0}^{t}\mathrm{d}s\,H(s)$ is bounded by

$\|S\|_{\infty,t}\leq(1+\|H\|_{1,t})\|U-1\|_{\infty,t}.$ (C.5)

###### Proof.

We have

$\displaystyle U(t)-1$
$\displaystyle=-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,H(s)U(s)$
$\displaystyle=-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,H(s)[U(s)-1]-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,H(s)$
$\displaystyle=-\mathrm{i}\int_{0}^{t}\mathrm{d}s\,H(s)[U(s)-1]-\mathrm{i}S(t),$
(C.6)

that is,

$S(t)=\mathrm{i}[U(t)-1]-\int_{0}^{t}\mathrm{d}s\,H(s)[U(s)-1].$ (C.7)

By taking the norm the result follows. ∎

##  Appendix D Non-isospectral Hamiltonians

Here, we prove that, if the constant generators of two evolutions are not
isospectral, the two evolutions eventually diverge, regardless of how close
the generators may be.

######  Proposition 4.

Let $H$ and $G$ be two (not necessarily bounded) self-adjoint operators with
pure point spectrum. If $H$ and $G$ are not unitarily equivalent, then

$\sup_{t\in\mathbb{R}}\|\mathrm{e}^{-\mathrm{i}tH}-\mathrm{e}^{-\mathrm{i}tG}\|\geq\sqrt{2}.$
(D.1)

###### Proof.

  1. 1.

Let $g$ be an eigenvalue of $G$ and let $\varphi$ be its associated normalized
eigenvector. One has

$\displaystyle\delta_{\varphi}^{2}(t)$
$\displaystyle=\|(\mathrm{e}^{-\mathrm{i}tH}-\mathrm{e}^{-\mathrm{i}tG})\varphi\|^{2}\vphantom{\int}$
$\displaystyle=\|(\mathrm{e}^{-\mathrm{i}t(H-g)}-1)\varphi\|^{2}\vphantom{\int}$
$\displaystyle=\int|\mathrm{e}^{-\mathrm{i}t(\lambda-g)}-1|^{2}\,\mathrm{d}\mu_{\varphi}(\lambda)$
$\displaystyle=4\int\sin^{2}\\!\left(\frac{\lambda-g}{2}t\right)\mathrm{d}\mu_{\varphi}(\lambda)$
$\displaystyle=2-2\int\cos[(\lambda-g)t]\,\mathrm{d}\mu_{\varphi}(\lambda),$
(D.2)

where $\mu_{\varphi}$ is the spectral measure of $H$ at $\varphi$. Thus, for
any $T>0$,

$\sup_{t\in\mathbb{R}}\delta_{\varphi}^{2}(t)\geq\frac{1}{T}\int_{0}^{T}\delta_{\varphi}^{2}(t)\,\mathrm{d}t=2-2\int\operatorname{sinc}[(\lambda-g)T]\,\mathrm{d}\mu_{\varphi}(\lambda)$
(D.3)

by Fubini’s theorem. Notice that $\operatorname{sinc}(\lambda
T)\to\mathbf{1}_{\\{0\\}}(\lambda)$ pointwise as $T\to+\infty$, where
$\mathbf{1}_{\Omega}(\lambda)$ is the indicator function of the set $\Omega$.
Therefore, by taking the limit $T\to+\infty$, by the dominated convergence
theorem one gets

$\sup_{t\in\mathbb{R}}\delta_{\varphi}^{2}(t)\geq 2-2\mu_{\varphi}(\\{g\\}).$
(D.4)

Notice that this result is valid for every self-adjoint $H$, irrespective of
its spectrum. If $H$ has a pure point spectrum, with the spectral resolution

$H=\sum_{j}h_{j}P_{j},$ (D.5)

then the spectral measure at $\varphi$ of a set $\Omega\subset\mathbb{R}$
reads

$\mu_{\varphi}(\Omega)=\sum_{j\,:\,h_{j}\in\Omega}\|P_{j}\varphi\|^{2},$ (D.6)

which yields

$\delta_{\varphi}^{2}(t)=2-2\sum_{j}\|P_{j}\varphi\|^{2}\cos[(h_{j}-g)t],$
(D.7)

and

$\mu_{\varphi}(\\{g\\})=\sum_{j}\delta_{h_{j},g}\|P_{j}\varphi\|^{2},$ (D.8)

with $\delta_{h,g}$ denoting the Kronecker delta.

  2. 2.

Since by assumption $H$ and $G$ are not unitarily equivalent, either they have
different spectra or they have the same spectrum with different
multiplicities. In the first case, let $g$ be an eigenvalue of, say, $G$ with
$g\notin\mathop{\mathrm{spec}}H$. Then, $\mu_{\varphi}(\\{g\\})=0$ and, by
point 1 above,

$\sup_{t\in\mathbb{R}}\|\mathrm{e}^{-\mathrm{i}tH}-\mathrm{e}^{-\mathrm{i}tG}\|\geq\sup_{t}\delta_{\varphi}(t)\geq\sqrt{2}.$
(D.9)

In the second case, there is a common eigenvalue $g$ of $G$ and $H$ with
different multiplicities. Let $P_{j_{1}}$ and $Q_{j_{2}}$, for some $j_{1}$
and $j_{2}$, be the associated eigenprojections of $H$ and $G$, respectively,
i.e. $h_{j_{1}}=g_{j_{2}}=g$. We have, say,
$\dim\mathop{\mathrm{ran}}Q_{j_{2}}>\dim\mathop{\mathrm{ran}}P_{j_{1}}$, and
thus there exists a unit vector $\varphi\in\mathop{\mathrm{ran}}Q_{j_{2}}$,
with $P_{j_{1}}\varphi=0$. Therefore,

$\mu_{\varphi}(\\{g\\})=\|P_{j_{1}}\varphi\|^{2}=0,$ (D.10)

and (D.9) holds again.

∎

## References

  * [1] Dynamical Systems III, Vol. 3 of Encyclopaedia of Mathematical Sciences, edited by V. I. Arnold (Springer-Verlag, Berlin, 1988).
  * [2] J. A. Sanders, F. Verhulst, and J. Murdock, Averaging Methods in Nonlinear Dynamical Systems, 2nd ed. (Springer, New York, 2007).
  * [3] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa, Generalized Adiabatic Theorem and Strong-Coupling Limits, Quantum 3, 152 (2019).
  * [4] L. Allen and J. H. Eberly, Optical Resonance and Two-Level Atoms, revised ed. (Dover, New York, 1987).
  * [5] M. Mehring, Principles of High Resolution NMR in Solids, 2nd ed. (Springer, Berlin, 1983).
  * [6] P. Krantz, M. Kjaergaard, F. Yan, T. P. Orlando, S. Gustavsson, and W. D. Oliver, A Quantum Engineer’s Guide to Superconducting Qubits, Appl. Phys. Rev. 6, 021318 (2019).
  * [7] A. Laucht, S. Simmons, R. Kalra, G. Tosi, J. P. Dehollain, J. T. Muhonen, S. Freer, F. E. Hudson, K. M. Itoh, D. N. Jamieson, J. C. McCallum, A. S. Dzurak, and A. Morello, Breaking the Rotating Wave Approximation for a Strongly Driven Dressed Single-Electron Spin, Phys. Rev. B 94, 161302 (2016).
  * [8] U. Haeberlen and J. S. Waugh, Coherent Averaging Effects in Magnetic Resonance, Phys. Rev. 175, 453 (1968).
  * [9] O. Gamel and D. F. V. James, Time-Averaged Quantum Dynamics and the Validity of the Effective Hamiltonian Model, Phys. Rev. A 82, 052106 (2010).
  * [10] G. W. Series, A Semi-Classical Approach to Radiation Problems, Phys. Rep. 43, 1 (1978).
  * [11] D. Zeuch, F. Hassler, J. J. Slim, and D. P. DiVincenzo, Exact Rotating Wave Approximation, Ann. Phys. 423, 168327 (2020).
  * [12] S. Blanes, F. Casas, J. A. Oteo, and J. Ros, The Magnus Expansion and Some of Its Applications, Phys. Rep. 470, 151 (2009).
  * [13] Q. Xie and W. Hai, Analytical Results for a Monochromatically Driven Two-Level System, Phys. Rev. A 82, 032117 (2010).
  * [14] H.-J. Schmidt, J. Schnack, and M. Holthaus, Floquet Theory of the Analytical Solution of a Periodically Driven Two-Level System, arXiv:1809.00558 [quant-ph] (2018).
  * [15] T. Chambrion, Periodic Excitations of Bilinear Quantum Systems, Automatica 48, 2040 (2012).
  * [16] N. Augier, U. Boscain, and M. Sigalotti, On the Compatibility between the Adiabatic and the Rotating Wave Approximations in Quantum Control, in Proceedings of the 2019 IEEE 58th Conference on Decision and Control (CDC) (IEEE, New York, 2019), pp. 2292–2297.
  * [17] R. Robin, N. Augier, U. Boscain, and M. Sigalotti, Ensemble Qubit Controllability with a Single Control via Adiabatic and Rotating Wave Approximations, arXiv:2003.05831 [math-ph] (2020).
  * [18] N. Augier, U. Boscain, and M. Sigalotti, Effective Adiabatic Control of a Decoupled Hamiltonian Obtained by Rotating Wave Approximation, arXiv:2005.02737 [math.OC] (2020).
  * [19] D. D’Alessandro, Introduction to Quantum Control and Dynamics, 2nd ed. (CRC Press, Boca Raton, FL, 2022).
  * [20] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa, Eternal Adiabaticity in Quantum Evolution, Phys. Rev. A 103, 032214 (2021).
  * [21] A. Messiah, Quantum Mechanics (Dover, New York, 2017).
  * [22] T. Kato, On the Adiabatic Theorem of Quantum Mechanics, J. Phys. Soc. Jpn. 5, 435 (1950).
  * [23] J. E. Avron and A. Elgart, Adiabatic Theorem without a Gap Condition, Commun. Math. Phys. 203, 445 (1999).
  * [24] A. Joye, General Adiabatic Evolution with a Gap Condition, Commun. Math. Phys. 275, 139 (2007).
  * [25] J. E. Avron, M. Fraas, G. M. Graf, and P. Grech, Adiabatic Theorems for Generators of Contracting Evolutions, Commun. Math. Phys. 314, 163 (2012).
  * [26] P. Facchi and S. Pascazio, Quantum Zeno Subspaces, Phys. Rev. Lett. 89, 080401 (2002).
  * [27] P. Facchi and S. Pascazio, Quantum Zeno Dynamics: Mathematical and Physical Aspects, J. Phys. A: Math. Theor. 41, 493001 (2008).
  * [28] P. Facchi and M. Ligabò, Quantum Zeno Effect and Dynamics, J. Math. Phys. 51, 022103 (2010).
  * [29] L. S. Schulman, Continuous and Pulsed Observations in the Quantum Zeno Effect, Phys. Rev. A 57, 1509 (1998).
  * [30] P. Facchi, S. Tasaki, S. Pascazio, H. Nakazato, A. Tokuse, and D. A. Lidar, Control of Decoherence: Analysis and Comparison of Three Different Strategies, Phys. Rev. A 71, 022302 (2005).
  * [31] E. W. Streed, J. Mun, M. Boyd, G. K. Campbell, P. Medley, W. Ketterle, and D. E. Pritchard, Continuous and Pulsed Quantum Zeno Effect, Phys. Rev. Lett. 97, 260402 (2006).
  * [32] F. Schäfer, I. Herrera, S. Cherukattil, C. Lovecchio, F. S. Cataliotti, F. Caruso, and A. Smerzi, Experimental Realization of Quantum Zeno Dynamics, Nat. Commun. 5, 3194 (2014).
  * [33] Z. Gong, N. Yoshioka, N. Shibata, and R. Hamazaki, Universal Error Bound for Constrained Quantum Dynamics, Phys. Rev. Lett. 124, 210606 (2020).
  * [34] Z. Gong, N. Yoshioka, N. Shibata, and R. Hamazaki, Error Bounds for Constrained Dynamics in Gapped Quantum Systems: Rigorous Results and Generalizations, Phys. Rev. A 101, 052122 (2020).
  * [35] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa, Kolmogorov-Arnold-Moser Stability for Conserved Quantities in Finite-Dimensional Quantum Systems, Phys. Rev. Lett. 126, 150401 (2021).
  * [36] R. P. Feynman, Space-Time Approach to Non-Relativistic Quantum Mechanics, Rev. Mod. Phys. 20, 367 (1948).
  * [37] B. Simon, Functional Integration and Quantum Physics (Academic Press, New York, 1979), Vol. 86.
  * [38] M. Suzuki, Decomposition Formulas of Exponential Operators and Lie Exponentials with Some Applications to Quantum Mechanics and Statistical Physics, J. Math. Phys. 26, 601 (1985).
  * [39] M. Suzuki, General Theory of Fractal Path Integrals with Applications to Many-Body Theories and Statistical Physics, J. Math. Phys. 32, 400 (1991).
  * [40] L. M. Sieberer, T. Olsacher, A. Elben, M. Heyl, P. Hauke, F. Haake, and P. Zoller, Digital Quantum Simulation, Trotter Errors, and Quantum Chaos of the Kicked Top, npj Quant. Inf. 5, 78 (2019).
  * [41] M. C. Tran, Y. Su, D. Carney, and J. M. Taylor, Faster Digital Quantum Simulation by Symmetry Protection, PRX Quantum 2, 010323 (2021).
  * [42] D. Burgarth, P. Facchi, G. Gramegna, and S. Pascazio, Generalized Product Formulas and Quantum Control, J. Phys. A: Math. Theor. 52, 435301 (2019).
  * [43] D. Burgarth, P. Facchi, H. Nakazato, S. Pascazio, and K. Yuasa, Quantum Zeno Dynamics from General Quantum Operations, Quantum 4, 289 (2020).
  * [44] G. Teschl, Ordinary Differential Equations and Dynamical Systems (American Mathematical Society, Rhode Island, 2012).
  * [45] T. Albash and D. A. Lidar, Adiabatic Quantum Computation, Rev. Mod. Phys. 90, 015002 (2018).
  * [46] R. S. Strichartz, A Guide to Distribution Theory and Fourier Transforms (World Scientific, Singapore, 2003).
  * [47] S. Jansen, M.-B. Ruskai, and R. Seiler, Bounds for the Adiabatic Approximation with Applications to Quantum Computation, J. Math. Phys. 48, 102111 (2007).
  * [48] D. Burgarth, P. Facchi, V. Giovannetti, H. Nakazato, S. Pascazio, and K. Yuasa, Non-Abelian Phases from Quantum Zeno Dynamics, Phys. Rev. A 88, 042107 (2013).
  * [49] J. Z. Bernád, Dynamical Control of Quantum Systems in the Context of Mean Ergodic Theorems, J. Phys. A: Math. Theor. 50, 065303 (2017).
  * [50] A. M. Childs, Y. Su, M. C. Tran, N. Wiebe, and S. Zhu, Theory of Trotter Error with Commutator Scaling, Phys. Rev. X 11, 011020 (2021).
  * [51] T. G. Kurtz, A Random Trotter Product Formula, Proc. Amer. Math. Soc. 35, 147 (1972).
  * [52] E. Campbell, Random Compiler for Fast Hamiltonian Simulation, Phys. Rev. Lett. 123, 070503 (2019).
  * [53] A. M. Childs, A. Ostrander, and Y. Su, Faster Quantum Simulation by Randomization, Quantum 3, 182 (2019).
  * [54] Y. Ouyang, D. R. White, and E. T. Campbell, Compilation by Stochastic Hamiltonian Sparsification, Quantum 4, 235 (2020).
  * [55] C.-F. Chen, H.-Y. Huang, R. Kueng, and J. A. Tropp, Concentration for Random Product Formulas, PRX Quantum 2, 040305 (2021).
  * [56] L. Viola and S. Lloyd, Dynamical Suppression of Decoherence in Two-State Quantum Systems, Phys. Rev. A 58, 2733 (1998).
  * [57] L. Viola, E. Knill, and S. Lloyd, Dynamical Decoupling of Open Quantum Systems, Phys. Rev. Lett. 82, 2417 (1999).
  * [58] D. Vitali and P. Tombesi, Using Parity Kicks for Decoherence Control, Phys. Rev. A 59, 4178 (1999).
  * [59] P. Zanardi, Symmetrizing Evolutions, Phys. Lett. A 258, 77 (1999).
  * [60] L.-M. Duan and G.-C. Guo, Suppressing Environmental Noise in Quantum Computation through Pulse Control, Phys. Lett. A 261, 139 (1999).
  * [61] L. Viola, Quantum Control via Encoded Dynamical Decoupling, Phys. Rev. A 66, 012307 (2002).
  * [62] C. Uchiyama and M. Aihara, Multipulse Control of Decoherence, Phys. Rev. A 66, 032313 (2002).
  * [63] L. Viola and E. Knill, Random Decoupling Schemes for Quantum Dynamical Control and Error Suppression, Phys. Rev. Lett. 94, 060502 (2005).
  * [64] L. F. Santos and L. Viola, Enhanced Convergence and Robust Performance of Randomized Dynamical Decoupling, Phys. Rev. Lett. 97, 150501 (2006).
  * [65] L. Viola and L. F. Santos, Randomized Dynamical Decoupling Techniques for Coherent Quantum Control, J. Mod. Opt. 53, 2559 (2006).
  * [66] L. F. Santos and L. Viola, Advantages of Randomization in Coherent Quantum Dynamical Control, New J. Phys. 10, 083009 (2008).
  * [67] R. Hillier, C. Arenz, and D. Burgarth, A Continuous-Time Diffusion Limit Theorem for Dynamical Decoupling and Intrinsic Decoherence, J. Phys. A: Math. Theor. 48, 155301 (2015).
  * [68] A. Hahn, D. Burgarth, and K. Yuasa, Unification of Random Dynamical Decoupling and the Quantum Zeno Effect, New J. Phys. (in press).
  * [69] P. Facchi, D. A. Lidar, and S. Pascazio, Unification of Dynamical Decoupling and the Quantum Zeno Effect, Phys. Rev. A 69, 032314 (2004).


# Fourier Disentangled Multimodal Prior Knowledge Fusion for Red Nucleus
Segmentation in Brain MRI

###### Abstract

Early and accurate diagnosis of parkinsonian syndromes is critical to provide
appropriate care to patients and for inclusion in therapeutic trials. The red
nucleus is a structure of the midbrain that plays an important role in these
disorders. It can be visualized using iron-sensitive magnetic resonance
imaging (MRI) sequences. Different iron-sensitive contrasts can be produced
with MRI. Combining such multimodal data has the potential to improve
segmentation of the red nucleus. Current multimodal segmentation algorithms
are computationally consuming, cannot deal with missing modalities and need
annotations for all modalities. In this paper, we propose a new model that
integrates prior knowledge from different contrasts for red nucleus
segmentation. The method consists of three main stages. First, it disentangles
the image into high-level information representing the brain structure, and
low-frequency information representing the contrast. The high-frequency
information is then fed into a network to learn anatomical features, while the
list of multimodal low-frequency information is processed by another module.
Finally, feature fusion is performed to complete the segmentation task. The
proposed method was used with several iron-sensitive contrasts (iMag, QSM,
R2*, SWI). Experiments demonstrate that our proposed model substantially
outperforms a baseline UNet model when the training set size is very small.

Index Terms—  Multimodal Fusion, Segmentation, Prior Knowledge, Deep Learning,
Red Nucleus, Parkinson Disease

##  1 Introduction

Parkinsonism is clinically defined by the association of bradykinesia, plastic
rigidity and/or asymmetrical resting tremor. Parkinson’s disease is the most
common neurodegenerative cause of parkinsonism. Atypical parkinsonism includes
tauopathies (progressive supranuclear palsy [PSP] and corticobasal
degeneration) and synucleinopathies (multiple system atrophy [MSA] and
dementia with Lewy body [DLB]) [1]. Early and accurate diagnosis of
parkinsonian syndromes is critical to provide appropriate care to patients and
for inclusion in therapeutic trials. The red nucleus is a small nucleus
located in the midbrain. It is a region of particular interest considering its
known pathophysiological changes in PSP [2]. Iron-sensitive MRI can be used to
assess iron deposition in the brain, particularly in the red nucleus. It
provides promising biomarkers for differentiating parkinsonism and for
longitudinal monitoring of disease progression and treatment effects [3].
Different iron-sensitive contrasts can be derived from multi-echo T2*
acquisitions by combining the phase and magnitude information. These contrasts
include quantitative susceptibility mapping (QSM), R2* relaxation rate mapping
(R2*), magnitude images (iMag) and susceptibility-weighted images (SWI).
Combining these different contrasts has the potential to improve the
segmentation of the red nucleus. In the following, we will refer to these
different contrasts as modalities, keeping in mind that some are actually
contrasts derived from a single sequence (iMag, R2* and QSM) while others are
a separate sequence.

Multimodal medical image fusion has attracted great interest in recent years
(please see [4, 5] for reviews). Chartsias et al. [6] proposed a method that
can improve segmentation accuracy of the modality of interest by learning to
leverage the information from different modalities. The key idea is to learn a
disentangled decomposition into anatomical and imaging factors. This research
demonstrates that using multimodal data through disentangled representations
is an effective way to increase the robustness of the model and improve
performance. Guo et al. [7] proposed a model architecture to fuse cross-
modality representation at feature learning, classifier and decision-making
level, respectively. They conducted experiments using 4 imaging modalities and
showed improvement compared to the baseline network using a single modality.
Chen et al. [8] proposed a multimodal learning model to decompose the
disentanglement feature into modality-specific appearance and modality-
invariant content, and share representation between different disentangled
contents. It can achieve good performance on brain tumor segmentation tasks.
All of the above studies demonstrated the benefits of multimodal learning, but
have not effectively exploited the similarities and differences in information
between modalities. Moreover, these methods require to have all modalities as
input. However, it is frequent that one or more modalities are missing, noisy,
or lack annotated data in real-world tasks [9].

In this paper, we propose a multimodal fusion approach to segment the red
nucleus based on various iron-sensitive modalities (iMag, QSM, R2*, SWI). The
underlying hypothesis of our work is that the anatomical structure of the
target will be approximately constant, and that the cross-modal variation will
mainly be in contrast. Such hypothesis is of course not universally true
(there may be sequences that do not show or show only part of some anatomical
structures), however we believe it holds in our application as well as in
several others. If the data can be decoupled into domain-invariant and domain-
shift parts, then a large amount of computational resources can be saved by
learning the domain-invariant part only once. Moreover, another advantage is
that the approach does not require all modalities to be present for all
subjects or to be annotated. Only the target modality needs to be
systematically present and annotated.

##  2 Methods

An overview of our approach can be seen in Figure 1. To ease understanding, we
show an example where the target is QSM data but the same concept can be used
with another modality as a target. The inputs of the approach are the target
modality (here QSM) and a list of other modalities that will bring low-
frequency information (here iMag, R2* and SWI). The proposed model first
disentangles each image into high- and low-frequency parts. The high-frequency
(here $H_{qsm}$) of the target image is passed through a segmentation backbone
(here a U-Net), and the low-frequency list (here
$L_{qsm},L_{imag},L_{R2*},L_{swi}$) is input to a shared weighted
convolutional layer, and is then fused with the output of the segmentation
backbone.

Fig. 1: Processing flow of the proposed method. The figure displays an example
where the target modality is QSM for easy understanding but the principle
would be similar with a different target.

###  2.1 Fourier domain disentangling

Given $p$ modalities ${M_{1},M_{2},\ldots,M_{p}}$ with samples $i_{M_{j}}\in
I_{M_{j}}$, where $I_{M_{j}}\subset\mathbb{R}^{H\times W}$ is a set of images,
and $H$ and $W$ are the height and width. The disentangling operation is
achieved by first transferring to Fourier space, and separating the high- and
low-frequency components as follows:

$\displaystyle{\scriptstyle\mathcal{F}({i_{M_{j}}})=H^{\theta}_{F}({i_{M_{j}}})+L^{\theta}_{F}({i_{M_{j}}})}$
(1) $\displaystyle{\scriptstyle
L^{\theta}_{F}({i_{M_{j}}})=\mathcal{F}({i_{M_{j}}})\left[\frac{H\times(1-\theta)}{2}:\frac{H\times(1+\theta)}{2},\frac{W\times(1-\theta)}{2}:\frac{W\times(1+\theta)}{2}\right]}$
$\displaystyle{\scriptstyle
H^{\theta}_{F}({i_{M_{j}}})=\mathcal{F}({i_{M_{j}}})-L_{F}({i_{M_{j}}})}$

where $\mathcal{F}$ represents the Fourier transform, $L^{\theta}_{F}$ is the
extraction of the low frequency part (in Fourier space, it is the center of
the signal), $H^{\theta}_{F}$ is the extraction of the high frequency part
(defined by subtracting the low-frequency part) and $\theta\in(0,1)$ is a
parameter that controls the high-/low-frequency separation.

Our approach can get benefit from the low-frequency part of other modalities,
even if not all modalities are present for all subjects. Actually, as will be
shown in the experiments, even when only one subject has a given modality it
can still be beneficial. The low-frequency component, taken from a handful of
subjects, can be seen as prior-knowledge.

The visualization of the high- and low-frequency parts of the QSM can be seen
in Figure 2. The high-frequency part basically represents the brain structure
and the low-frequency part the contrast of the modality. The disentangled
representation can help the model to learn brain structure and contrast
information separately.

Fig. 2: Disentangled representation into high- and low-frequency parts of a
QSM MRI. Note that we pad the low-frequency information to the original image
size with 0 and then perform an inverse Fourier transform for visualization.
Parameter $\theta$ (controlling the high-/low- frequency separation) is 0.1 in
this figure.

###  2.2 Model architecture

As previously mentioned, an image is first disentangled into high- and low-
frequency parts. For the target modality $M_{1}$, we further use both high-
and low-frequency parts $H^{1,\theta}_{F}$ and $L^{1,\theta}_{F}$ while, for
other modalities $M_{2},\ldots,M_{p}$, we only use the low-frequency part:
$L^{2,\theta}_{F},\ldots,L^{p,\theta}_{F}$. We then apply the inverse Fourier
transform ($\mathcal{F}^{-1}$) to obtain high- and low-frequency parts in
image space. The results are respectively denoted as $H^{1,\theta}_{F^{-1}}$
and $L^{j,\theta}_{{F}^{-1}}$ where $j\in\\{1,\ldots,p\\}$. Then,
$H^{1,\theta}_{{F}^{-1}}$ is fed to the backbone (a 3D-UNet [10] in our
experiments), the output being denoted as $O^{1,\theta}_{{F}^{-1}}$. The low-
frequency list $L^{1,\theta}_{{F}^{-1}},\ldots,L^{p,\theta}_{{F}^{-1}}$ is
processed by a shared convolutional layer, the outputs being denoted as
$O^{1,\theta}_{L_{F^{-1}}},\ldots,O^{p,\theta}_{L_{F^{-1}}}$.
$O^{1,\theta}_{L_{F^{-1}}},\ldots,O^{p,\theta}_{L_{F^{-1}}}$ are then fused
with $O^{1,\theta}_{H_{F^{-1}}}$ one by one. Each fusion operation sets the
center values of $O^{1,\theta}_{H_{F^{-1}}}$ to $O^{j,\theta}_{L_{F^{-1}}}$:
$O_{H_{F^{-1}}}^{1,\theta}[\frac{H\times(1-\theta)}{2}:\frac{H\times(1+\theta)}{2},\frac{W\times(1-\theta)}{2}:\frac{W\times(1+\theta)}{2}]=O^{j,\theta}_{L_{F}^{-1}}$.
The high and low-frequency fusion operation thus produces $p$ outputs. They
are further processed with convolution, activation function and dropout,
producing $p$ predictions. The final loss function is the sum of differences
between the $p$ predictions and the ground truth.

##  3 Experiments and Results

###  3.1 Dataset

We used a dataset comprising 18 healthy subjects, 46 patients with early
Parkinson’s disease (i.e. with a disease duration below 4 years), and 16
patients with prodromal parkinsonism (idiopathic rapid eye movement sleep
behavior disorder-iRBD), recruited between May 2015 and January 2019 as part
of the ICEBERG cohort.

Each patient had a multi echo 3D GRE (12 echo times ranging from 4 ms to 37
ms) with a full brain coverage at an isotropic voxel resolution of 1 mm3. The
iMag contrast is the average over the echos of the magnitude volumes. The R2*
is the exponential fit of the echo time decay, and the QSM is the fit of
quantitative susceptibility. We performed different experiments in which the
target image to segment was either QSM, R2* or iMag. Additionally, we
extracted low-frequency information from an SWI sample (without annotated
ground truth) to enrich the low-frequency list. Figure 3 displays examples of
these fours modalities. The dataset was split at the participant-level into
training, validation and test sets that contains 51, 13 and 16 participants
separately. The split was stratified for diagnostic class.

Fig. 3: Visualization of four modalities (QSM, R2*, iMag and SWI). The
segmentation of the red nucleus is superimposed on the first three images.

###  3.2 Implementation details

Our code is developed based on the PyTorch framework [11]. We used the open-
source Python library TorchIO [12] for reshaping images to the same size for
each task and for min-max normalization. We used Adam [13] as optimizer, Dice
as loss function and a batch size of 1.

###  3.3 Experiments

In the following experiments, we compare the model performance for different
modality combinations and with different amounts of training data. In all
experiments, the hyperparameters were identical.

We report performance at the voxel level for a total of 7 metrics: Dice score,
precision, recall, mean volume error rate (MVER), mean absolute volume error
rate (MAVER), Pearson’s $r$ and 95 percentile Hausdorff distance.

####  3.3.1 Extracting prior-knowledge from different modalities

We first experimented with varying the number of modalities bringing prior
knowledge (low-frequency part only) in addition to the target modality (both
high- and low-frequency parts). The prior knowledge may come from:

  * •

0 modality. We simply use the target modality (decoupling the data into high
and low-frequency). This is to evaluate the effect of disentanglement.

  * •

1 to 3 modality. For the modalities in our experiment, we mainly have two
situations, similar contrast (QSM and R2*) or inverse contrast (iMag and SWI).
These experiments can assess the effect of adding other modalities, including
the contrast relationship and the number of modalities.

Since we are mainly interested in the performance in the low sample size
regime, we only used 7.5% of the training set for the experiment. Results can
be seen in Table 1.

The improvements obtained when adding prior knowledge modalities are not
consistent: they depend on the metrics and the combination. However, there are
cases where adding more modalities substantially improves a given metric as
can for instance be seen for the 95 percentile Hausdorff.

####  3.3.2 Varying the size of the training set

We then evaluated the performance of our approach for different sizes of the
training set and compared it to a baseline approach (UNet). Specifically, the
models were trained with 7.5% (4 subjects), 15% (8 subjects), 30% (16
subjects), 50% (26 subjects) and 100% (51 subjects) of the training set, while
the validation and test sets were left unchanged. The low-frequency part is
only used as a prior knowledge and we actually used only one subject (from the
training set) for the modalities used for low-frequency. We perform this
experiment based on the best modalities combinations obtained for the Dice
score in the previous results shown in Table 1 (iMag: iMag+R2*+SWI; QSM:
QSM+SWI and R2*: iMag+QSM+R2*). The results obtained with variable training
set sizes are shown in Table 2. Our approach model outperforms the baseline in
the very small training size regime (7.5% data). This demonstrates the
potential of our algorithm with small training data sets, thus avoiding the
expensive labeling.

Table 1: Results of different modality combination when using only 7.5% of the
training set (4 subjects). Results are reported as $\text{mean}\pm\text{SEM}$
where SEM is the standard-error of the mean.

Task | Modality | Voxel level Performance
---|---|---
iMag | QSM | R2* | SWI | Dice | 95 Hausdorff | Precision | Recall | MVER | MAVER | Pearson’s r
iMag | ✓ |  |  |  | 0.69±0.02 | 38.45±6.01 | 0.6±0.03 | 0.82±0.02 | 0.41±0.07 | 0.41±0.07 | 0.7±0.02
✓ | ✓ |  |  | 0.58±0.03 | 54.74±3.91 | 0.45±0.03 | 0.84±0.02 | 0.98±0.14 | 0.98±0.14 | 0.61±0.02
✓ |  | ✓ |  | 0.69±0.02 | 27.19±6.05 | 0.6±0.02 | 0.8±0.03 | 0.35±0.05 | 0.35±0.05 | 0.69±0.02
✓ |  |  | ✓ | 0.7±0.02 | 20.32±5.57 | 0.63±0.02 | 0.8±0.03 | 0.28±0.05 | 0.29±0.05 | 0.7±0.02
✓ | ✓ | ✓ |  | 0.65±0.03 | 46.29±5.53 | 0.56±0.03 | 0.79±0.03 | 0.47±0.08 | 0.47±0.08 | 0.66±0.02
✓ | ✓ |  | ✓ | 0.67±0.02 | 30.77±6.14 | 0.59±0.03 | 0.79±0.02 | 0.38±0.07 | 0.39±0.07 | 0.68±0.02
✓ |  | ✓ | ✓ | 0.72±0.02 | 25.17±6.40 | 0.65±0.02 | 0.81±0.03 | 0.27±0.05 | 0.27±0.05 | 0.72±0.02
✓ | ✓ | ✓ | ✓ | 0.69±0.03 | 24.52±6.49 | 0.61±0.03 | 0.81±0.03 | 0.39±0.09 | 0.39±0.08 | 0.7±0.03
QSM |  | ✓ |  |  | 0.79±0.02 | 29.54±7.27 | 0.72±0.02 | 0.9±0.02 | 0.26±0.03 | 0.26±0.03 | 0.8±0.01
✓ | ✓ |  |  | 0.8±0.02 | 31.88±8.06 | 0.73±0.02 | 0.89±0.02 | 0.23±0.04 | 0.24±0.03 | 0.8±0.02
| ✓ | ✓ |  | 0.79±0.02 | 11.45±5.76 | 0.73±0.02 | 0.86±0.03 | 0.19±0.05 | 0.23±0.04 | 0.79±0.02
| ✓ |  | ✓ | 0.8±0.01 | 12.48±6.06 | 0.71±0.02 | 0.92±0.01 | 0.31±0.03 | 0.31±0.03 | 0.8±0.01
✓ | ✓ | ✓ |  | 0.77±0.02 | 31.62±7.77 | 0.7±0.02 | 0.87±0.02 | 0.26±0.04 | 0.26±0.04 | 0.78±0.01
✓ | ✓ |  | ✓ | 0.78±0.02 | 38.2±7.38 | 0.71±0.02 | 0.88±0.02 | 0.26±0.05 | 0.27±0.05 | 0.79±0.02
| ✓ | ✓ | ✓ | 0.75±0.02 | 44.69±6.86 | 0.65±0.03 | 0.9±0.01 | 0.41±0.06 | 0.41±0.06 | 0.76±0.02
✓ | ✓ | ✓ | ✓ | 0.76±0.02 | 37.06±6.84 | 0.68±0.02 | 0.86±0.02 | 0.28±0.04 | 0.28±0.04 | 0.76±0.02
R2* |  |  | ✓ |  | 0.64±0.03 | 44.11±6.94 | 0.56±0.03 | 0.77±0.04 | 0.39±0.06 | 0.39±0.06 | 0.65±0.03
✓ |  | ✓ |  | 0.58±0.03 | 53.64±5.83 | 0.47±0.03 | 0.8±0.04 | 0.76±0.08 | 0.76±0.08 | 0.61±0.03
| ✓ | ✓ |  | 0.62±0.03 | 53.5±4.86 | 0.5±0.02 | 0.81±0.03 | 0.65±0.07 | 0.65±0.07 | 0.64±0.02
|  | ✓ | ✓ | 0.59±0.02 | 56.33±4.95 | 0.45±0.02 | 0.86±0.02 | 0.96±0.08 | 0.96±0.08 | 0.62±0.02
✓ | ✓ | ✓ |  | 0.68±0.02 | 37.82±6.71 | 0.58±0.03 | 0.82±0.02 | 0.44±0.06 | 0.44±0.06 | 0.69±0.02
✓ |  | ✓ | ✓ | 0.58±0.04 | 59.51±3.15 | 0.47±0.03 | 0.76±0.05 | 0.65±0.07 | 0.65±0.07 | 0.59±0.04
| ✓ | ✓ | ✓ | 0.64±0.03 | 46.64±4.48 | 0.55±0.03 | 0.79±0.03 | 0.49±0.07 | 0.49±0.07 | 0.65±0.03
✓ | ✓ | ✓ | ✓ | 0.64±0.02 | 45.5±4.81 | 0.51±0.03 | 0.86±0.02 | 0.72±0.07 | 0.72±0.07 | 0.66±0.02

Table 2: The baseline model (U-Net) and the proposed approach are compared
when using different subsets of size n=4, 8, 16, 26 or 51 subjects,
respectively corresponding to 7.5%, 15%, 30%, 50% and 100% of the training
set. Results are reported as $\text{mean}\pm\text{SEM}$ where SEM is the
standard-error of the mean.

Model | Task | n | Dice | 95 Hausdorff | Precision | Recall | MVER | MAVER | Pearson’s r
---|---|---|---|---|---|---|---|---|---
UNet | IMag | 4 | 0.37±0.02 | 66.09±3.10 | 0.23±0.01 | 0.93±0.02 | 3.16±0.17 | 3.16±0.17 | 0.46±0.02
8 | 0.33±0.02 | 59.02±1.68 | 0.2±0.01 | 0.98±0.01 | 4.19±0.24 | 4.19±0.24 | 0.44±0.01
16 | 0.76±0.02 | 37.59±5.33 | 0.66±0.03 | 0.92±0.01 | 0.42±0.05 | 0.42±0.05 | 0.78±0.02
26 | 0.83±0.01 | 3.88±1.77 | 0.76±0.02 | 0.93±0.01 | 0.24±0.03 | 0.24±0.03 | 0.84±0.01
51 | 0.8±0.01 | 9.42±3.73 | 0.69±0.02 | 0.96±0.01 | 0.41±0.04 | 0.41±0.04 | 0.82±0.01
QSM | 4 | 0.46±0.03 | 70.04±1.52 | 0.31±0.03 | 0.93±0.03 | 2.27±0.24 | 2.27±0.24 | 0.53±0.03
8 | 0.84±0.01 | 5.66±3.91 | 0.78±0.02 | 0.92±0.02 | 0.18±0.04 | 0.2±0.03 | 0.85±0.01
16 | 0.77±0.01 | 19.97±6.63 | 0.64±0.02 | 0.99±0.01 | 0.56±0.05 | 0.56±0.05 | 0.79±0.01
26 | 0.88±0.01 | 1.0±0.0 | 0.8±0.01 | 0.97±0.01 | 0.21±0.02 | 0.21±0.02 | 0.88±0.01
51 | 0.88±0.01 | 1.0±0.0 | 0.81±0.02 | 0.97±0.01 | 0.2±0.02 | 0.2±0.02 | 0.89±0.01
R2* | 4 | 0.34±0.03 | 69.58±1.68 | 0.21±0.02 | 0.83±0.03 | 3.21±0.27 | 3.21±0.27 | 0.42±0.02
8 | 0.85±0.01 | 1.12±0.09 | 0.78±0.02 | 0.92±0.01 | 0.19±0.03 | 0.19±0.03 | 0.85±0.01
16 | 0.86±0.01 | 5.54±4.55 | 0.82±0.02 | 0.91±0.01 | 0.13±0.03 | 0.14±0.03 | 0.86±0.01
26 | 0.85±0.01 | 1.0±0.0 | 0.77±0.01 | 0.95±0.01 | 0.23±0.02 | 0.23±0.02 | 0.86±0.01
51 | 0.86±0.02 | 1.0±0.0 | 0.79±0.04 | 0.95±0.02 | 0.21±0.07 | 0.21±0.07 | 0.87±0.02
Ours | IMag | 4 | 0.72±0.02 | 25.17±6.40 | 0.65±0.02 | 0.81±0.03 | 0.27±0.05 | 0.27±0.05 | 0.72±0.02
8 | 0.79±0.01 | 20.0±5.44 | 0.69±0.02 | 0.91±0.01 | 0.33±0.03 | 0.33±0.03 | 0.79±0.01
16 | 0.81±0.01 | 18.71±6.24 | 0.74±0.02 | 0.92±0.01 | 0.26±0.04 | 0.26±0.04 | 0.82±0.01
26 | 0.83±0.01 | 7.45±3.75 | 0.75±0.02 | 0.93±0.01 | 0.25±0.04 | 0.25±0.04 | 0.84±0.01
51 | 0.85±0.01 | 7.46±4.41 | 0.76±0.02 | 0.96±0.01 | 0.26±0.02 | 0.26±0.02 | 0.85±0.01
QSM | 4 | 0.8±0.01 | 12.48±6.06 | 0.71±0.02 | 0.92±0.01 | 0.31±0.03 | 0.31±0.03 | 0.8±0.01
8 | 0.85±0.01 | 7.47±4.83 | 0.78±0.02 | 0.95±0.01 | 0.22±0.04 | 0.22±0.04 | 0.86±0.01
16 | 0.86±0.01 | 1.0±0.0 | 0.77±0.02 | 0.97±0.01 | 0.27±0.03 | 0.27±0.03 | 0.87±0.01
26 | 0.87±0.01 | 1.0±0.0 | 0.8±0.01 | 0.96±0.01 | 0.2±0.02 | 0.2±0.02 | 0.88±0.01
51 | 0.88±0.01 | 1.0±0.0 | 0.82±0.01 | 0.96±0.01 | 0.18±0.02 | 0.18±0.02 | 0.89±0.01
R2* | 4 | 0.68±0.02 | 37.82±6.71 | 0.58±0.03 | 0.82±0.02 | 0.44±0.06 | 0.44±0.06 | 0.69±0.02
8 | 0.78±0.01 | 3.12±1.04 | 0.69±0.02 | 0.89±0.02 | 0.3±0.05 | 0.32±0.04 | 0.78±0.01
16 | 0.84±0.01 | 1.06±0.035 | 0.76±0.01 | 0.96±0.01 | 0.28±0.03 | 0.28±0.03 | 0.85±0.01
26 | 0.84±0.01 | 1.03±0.03 | 0.77±0.02 | 0.95±0.01 | 0.24±0.03 | 0.24±0.03 | 0.85±0.01
51 | 0.86±0.01 | 1.0±0.0 | 0.79±0.02 | 0.95±0.01 | 0.2±0.03 | 0.2±0.03 | 0.86±0.01

##  4 Conclusion

In this paper, we proposed a new approach for segmentation of the red nucleus
from brain MRI data. The approach can leverage the low-frequency information
from different contrasts (or modalities). Learning high-frequency and low-
frequency information separately allows the model to focus on brain structure
throughout the learning process while contrast information is fused in the
final stage.

Our experiments show that the proposed approach leads to substantial
improvements when dealing with very small training sets. On the other hand,
the addition of other modalities led to mixed results, providing improvements
only part of the time. Future work shall focus on understanding how and why
our Fourier disentanglement process can lead to such massive improvements in
performance over the baseline in the very small training set regime.

Other future improvements would be to use convolution theory to do
multiplication in the Fourier domain to improve computational efficiency,
exploring other methods for representation disentanglement and applying the
framework to other datasets.

##  5 Compliance with ethical standards

The institutional ethical standard committee approved the study (CPP Paris
VI/RCB: 2014-A00725-42). All participants gave written informed consent.

##  6 Acknowledgments

The research leading to these results has received funding from the French
government under management of Agence Nationale de la Recherche as part of the
”Investissements d’avenir” program, reference ANR-19-P3IA-0001 (PRAIRIE 3IA
Institute) and reference ANR-10-IAIHU-06 (Agence Nationale de la
Recherche-10-IA Institut Hospitalo-Universitaire-6). The ICEBERG study is
supported by the European Research Council (ERC) under grant agreement No.
678304, the European Union’s Horizon 2020 research and innovation program
under grant agreement No. 826421 (TVB-Cloud), Agence Nationale de la Recherche
(ANR) under grant agreements ANR-10-IAIHU-06 (IHU ICM), ANR-11-INBS-0006, and
ANR-19-JPW2-000 (JPND E-DADS), association France Parkinson (PRECISE-PD
project), the Fondation d’Entreprise EDF, Biogen Inc., Fondation Thérèse and
René Planiol, Fondation Saint Michel. It received unrestricted support for
Research on Parkinson’s disease from Energipole (M. Mallart), M. Villain and
the Société Française de Médecine Esthétique (M. Legrand). Guanghui Fu is
supported by the Chinese Government Scholarship provided by China Scholarship
Council (CSC). Lydia Chougar is supported by a Poste d’accueil Inria/AP-HP.

## References

  * [1] Dennis W Dickson,  “Parkinson’s disease and parkinsonism: neuropathology,”  Cold Spring Harbor perspectives in medicine, vol. 2, no. 8, pp. a009258, 2012.
  * [2] David R Williams and Andrew J Lees,  “Progressive supranuclear palsy: clinicopathological concepts and diagnostic challenges,”  The Lancet Neurology, vol. 8, no. 3, pp. 270–279, 2009.
  * [3] Henrik Sjöström, Tobias Granberg, Eric Westman, and Per Svenningsson,  “Quantitative susceptibility mapping differentiates between parkinsonian disorders,”  Parkinsonism & related disorders, vol. 44, pp. 51–57, 2017.
  * [4] Tongxue Zhou, Su Ruan, and Stéphane Canu,  “A review: Deep learning for medical image segmentation using multi-modality fusion,”  Array, vol. 3, pp. 100004, 2019.
  * [5] Haithem Hermessi, Olfa Mourali, and Ezzeddine Zagrouba,  “Multimodal medical image fusion review: Theoretical background and recent advances,”  Signal Processing, vol. 183, pp. 108036, 2021.
  * [6] Agisilaos Chartsias, Giorgos Papanastasiou, Chengjia Wang, Scott Semple, David E Newby, Rohan Dharmakumar, and Sotirios A Tsaftaris,  “Disentangle, align and fuse for multimodal and semi-supervised image segmentation,”  IEEE transactions on medical imaging, vol. 40, no. 3, pp. 781–792, 2020.
  * [7] Zhe Guo, Xiang Li, Heng Huang, Ning Guo, and Quanzheng Li,  “Deep learning-based image segmentation on multimodal medical imaging,”  IEEE Transactions on Radiation and Plasma Medical Sciences, vol. 3, no. 2, pp. 162–169, 2019.
  * [8] Cheng Chen, Qi Dou, Yueming Jin, Hao Chen, Jing Qin, and Pheng-Ann Heng,  “Robust multimodal brain tumor segmentation via feature disentanglement and gated fusion,”  in International Conference on Medical Image Computing and Computer-Assisted Intervention. Springer, 2019, pp. 447–456.
  * [9] Anil Rahate, Rahee Walambe, Sheela Ramanna, and Ketan Kotecha,  “Multimodal co-learning: challenges, applications with datasets, recent advances and future directions,”  Information Fusion, vol. 81, pp. 203–239, 2022.
  * [10] Özgün Çiçek, Ahmed Abdulkadir, Soeren S Lienkamp, Thomas Brox, and Olaf Ronneberger,  “3d u-net: learning dense volumetric segmentation from sparse annotation,”  in International conference on medical image computing and computer-assisted intervention. Springer, 2016, pp. 424–432.
  * [11] Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, et al.,  “Pytorch: An imperative style, high-performance deep learning library,”  Advances in neural information processing systems, vol. 32, 2019\.
  * [12] Fernando Pérez-García, Rachel Sparks, and Sebastien Ourselin,  “Torchio: a python library for efficient loading, preprocessing, augmentation and patch-based sampling of medical images in deep learning,”  Computer Methods and Programs in Biomedicine, vol. 208, pp. 106236, 2021.
  * [13] Diederik P Kingma and Jimmy Ba,  “Adam: A method for stochastic optimization,”  arXiv preprint, vol. 1412.6980, 2014.