mathaOMLtxmimitmatha118

|  Molecular Dynamics Simulation of the Coalescence of Surfactant-Laden
Droplets†
---|---
|  Soheil Arbabi,a Piotr Deuar,a Mateusz Denys,a Rachid Bennacer,b Zhizhao
Che,c and Panagiotis E. Theodorakis∗a
|  We investigate the coalescence of surfactant-laden water droplets by using
several different surfactant types and a wide range of concentrations by means
of a coarse-grained model obtained by the statistical associating fluid
theory. Our results demonstrate in detail a universal mass transport mechanism
of surfactant across many concentrations and several surfactant types during
the process. Coalescence initiation is seen to occur via a single pinch due to
aggregation of surface surfactant, and its remnants tend to become engulfed in
part inside the forming bridge. Across the board we confirm the existence of
an initial thermal regime with constant bridge width followed by a later
inertial regime with bridge width scaling roughly as the square root of time,
but see no evidence of an intermediate viscous regime. Coalescence becomes
slower as surfactant concentration grows, and we see evidence of the
appearance of a further slowdown of a different nature for several times the
critical concentration. We anticipate that our results provide further
insights in the mechanisms of coalescence of surfactant-laden droplets.

##

††footnotetext: a Institute of Physics, Polish Academy of Sciences, Al.
Lotników 32/46, 02-668 Warsaw, Poland. ††footnotetext: b Université Paris-
Saclay, CentraleSupélec, ENS Paris-Saclay, CNRS, LMPS - Laboratoire de
Mécanique Paris-Saclay, 91190, Gif-sur-Yvette, France††footnotetext: c State
Key Laboratory of Engines, Tianjin University, 300350 Tianjin, China.
††footnotetext: † Electronic Supplementary Information (ESI) available:
[details of any supplementary information available should be included here].
See DOI: 10.1039/cXsm00000x/

##  1 Introduction

While ubiquitous in nature, droplet coalescence is also an important process
in various industrial applications, where the rate of coalescence can
determine their performance. For example, in the context of bio-related
microfluidic technologies, slowing down coalescence in bio-particle
encapsulation on lab-on-chip devices is often desirable, and can be achieved
by using various additives, such as surfactants.1, 2 In contrast, speeding up
the rate of coalescence by adding surfactant could be advantageous in
applications, such as coatings3 and superspreading.4, 5

Experimental, theoretical, and numerical studies of coalescing droplets have
thus far mainly focused on cases without additives, e.g. water or polymer
droplets.6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23,
24, 25, 26, 27, 28, 29, 30 From the point of view of numerical simulations,
these have by and large provided descriptions of the macroscopic and dynamic
properties of coalescence, 31, 32, 7, 33, 8, 20, 21, 22 but they generally
continue to suffer from inadequate resolution at the pinching point between
droplets at the initial stage of coalescence, despite progress in this area.12
Moreover, a detailed molecular-level description of the mass transport
mechanism of surfactant is beyond the reach of any continuum model.10, 34

Various attempts have been made to tackle this challenge. For example,
numerical models based on advection–diffusion equations and chemical kinetic
fluxes attempt to incorporate the mass transport of surfactant into the
equations, 35 but these are only as good as the assumptions put into the
model, while a detailed molecular description still remains out of reach for
continuum simulation. When it comes to experimental techniques,36, 37, 38, 34,
39, 40, 41, 42, 43, 44 high-speed imaging and particle-image velocimetry have
been applied in investigating coalescence of surfactant-laden droplets, and
have mainly focused on the macroscopic description of this phenomenon, much as
in the case of numerical simulation. Due to device limitations, capturing with
high resolution the early stages of coalescence during experiments poses a
major challenge. Both experiments and continuum modelling are unable to
provide a detailed description of the mass transport mechanism at the
molecular level. In the case of molecular simulations, all-atom molecular
dynamics (MD) has explained different microscopic aspects of coalescence, such
as the role of the capillary waves on droplets’ surface, but these have only
been in the context of pure water droplets,9 except for one recent study which
considered a few cases.45

On the whole, though, the effect of surfactant in droplet coalescence has
largely remained unexplored, with implications actually reaching beyond this,
for example, in the context of surfactant-laden coalescence between micelle
and bilayer46. The initial work45 suggests that it plays an important role,
e.g. affecting the coalescence rate and dampening internal droplet dynamics.
Reduction of the surface tension at the liquid–gas (LG) interface is usually
expected, but the type of surfactant (e.g., chemical groups and molecular
architecture) and its concentration are also anticipated to differently affect
the coalescence process. Differences are expected even when the concentration
is higher than the critical aggregation concentration (CAC) and surface
tension at the LG interface is expected to remain constant.47, 48, 49, 50, 51
On that account, unveiling the details of this phenomenon from a molecular-
level perspective is key for fundamentally understanding the underlying
mechanisms that can lead to tailor-made surfactant designs for applications.

Here we report on a more detailed and robust study of droplet coalescence
using the methods of Ref. 45, and a broader range of surfactants and their
concentrations. We employed a high-fidelity coarse-grained force-field
enabling us to simulate with MD the coalescence of surfactant-laden droplets.
We discuss the details of coalescence at each stage of the process by
comparing systems with different surfactant type and concentration. Our data
indicate that the underlying mechanisms of coalescence are consistent between
the different cases, which point to universal features of surfactant transport
and droplet dynamics for the coalescence of water droplets with surfactant.
Differences in the behaviour of different surfactants are also highlighted
when relevant. The dynamic and static characteristics are unveiled by
discussing both macroscopic- and molecular-level quantities, such as the
bridge growth dynamics and the velocity of droplet approach, as well as
describing the molecular transport of surfactant within the liquid phase.
Thus, we anticipate that the present study lays a comprehensive account of the
coalescence of surfactant-laden droplets.

In the following section, we provide background information for droplet
coalescence. Then, Section 3 gives details on our model and methodology, while
Section 4 presents and discusses the results of the coalescence simulations.
Section 5, draws some broader conclusions.

Fig. 1:  Coalescence of droplets with C10E4 surfactant at low concentration
(6.25 wt%, top of each panel) below $\text{CAC}=7.5$ wt% and concentrations
above CAC (24.18 wt%, bottom of each panel). Snapshots show the initial
approach of the droplets (a), their pinching and initial formation of the
bridge (b), see also ESI for a movie showing the pinching of the droplets),
configurations with a partially (c) and fully (d) developed bridge, and the
final equilibrium state (e). External or cross-section views are shown to
highlight the bulk, surface, and bridge structure of the droplets. Magnified
views of a micelle, a surfactant molecule, and a water bead are shown in panel
(a). Water vapour surrounding the droplets is not shown for the sake of
clarity.

##  2 Background

A macroscopic description of the droplet coalescence can generally be divided
into three stages (see Fig. 1). The approach of the droplets, which leads to
the initiation of the phenomenon (pinching) as a result of the inter-molecular
interactions; the subsequent growth of the bridge formed between the droplets;
and the final stage towards equilibrium, which manifests itself by the
formation of a single spherical droplet. The system is driven towards the
equilibrium state as it seeks to minimise the free-energy associated with the
surface tension at the droplet boundary. Hence, on the coarsest level, the
formation of larger droplets is energetically favourable. In particular, for a
droplet with $N$ beads, the droplet radius behaves as $R_{0}\propto N^{1/3}$
(Fig. 1a) since the volume scales as $N$, while surface tension force as
$\propto N^{2/3}$.

The rate of bridge growth after the initial contact is key for characterising
the coalescence process and can generally be described from the point of view
of fluid dynamics by two successive regimes, namely the viscous regime (VR)
and the inertial regime (IR).6, 52 Moreover, recent MD simulations have
uncovered the existence of a third, thermal regime (TR) during the droplet
pinching at the very initial stage of the coalescence process,9, 45 which will
be discussed later in more detail. The characteristic velocity in the VR can
be defined as $\varv_{\rm v}=\gamma/\eta$, where $\gamma$ is the surface
tension and $\eta$ the viscosity, which suggests that the capillary number
${\rm Ca}=\eta\varv_{\rm v}/\gamma\sim 1$, with the characteristic time scale
for droplet-sized effects being $t_{\rm v}=R_{0}\eta/\gamma$ 10, with $R_{0}$
the droplet radius as in Fig. 1. As coalescence proceeds to the subsequent IR,
the bridge flow can be characterised by the Weber number, namely ${\rm
We}=\rho\varv_{\rm i}^{2}R_{0}/\gamma\sim 1$ indicating the limit that inertia
effects will take over the surface tension effects, while $\rho$ is the
density of the fluid. The bridge velocity is found to scale as
$\varv_{i}\sim\sqrt{\gamma/\rho R_{0}}$, and thus the characteristic time
scale is $t_{\rm i}=\sqrt{\rho R_{0}^{3}/\gamma}$. For many purposes,
especially in the VR, the more relevant length scale is the radius of the
bridge, $b$ (Fig. 1), and the corresponding Reynolds number can be defined as
Re $=\rho\varv b/\eta$, which in the VR is $\rho\gamma b/\eta^{2}$. On the one
hand, since the bridge radius is very small in the initial stage, the viscous
forces are dominant regardless of the values of $\gamma$ and $\eta$ and Re
$\ll 1$. On the other hand, Re $\geq 1$ reflects the IR.53 Hence, it can be
taken that the crossover between the viscous and the inertial regimes would
take place for ${\rm Ca}\sim 1$ or ${\rm Re}\sim 1$. Taking the IR expression
$\varv_{\rm i}\sim\sqrt{\gamma/b\rho}$ using bridge size, we obtain an
estimate of the crossover bridge-radius, $b_{\rm c}=\eta^{2}/\rho\gamma$, and
a characteristic time $t_{\rm c}=\eta^{3}/\rho\gamma^{2}$, which for water
droplets would roughly correspond to $b_{\rm c}\sim 15$ nm and $t_{\rm c}=0.1$
ns, highlighting the fast time scales of coalescence.11

In the VR, where inter-molecular forces are playing a dominant role in pulling
the droplets together, a linear scaling has been proposed for the growth of
the bridge radius with time, i.e. $b\propto t$, as well as logarithmic
corrections, $b\propto t\ln t$.34, 10 In the case of the IR, a power-law
scaling law has been suggested for the bridge, namely $b\propto\sqrt{t}$.34,
10 Experimental studies on the coalescence of water droplets are consistent
with this and have shown that 0.7 ns after the first contact of the droplets,
drag forces give way to the inertial ones and the bridge radius has been
reported to scale with time as $b\propto(R_{0}\gamma/\rho)^{1/4}t^{1/2}$ where
$R_{0}$ is the initial droplet radius. 10, 11, 54, 53, 12 Other works have
proposed scaling regimes that depend on the ratio of characteristic scales to
the viscous length scale $l_{v}=\eta^{2}/\rho\gamma$.55, 52 Moreover, it is
argued that the inertia of the droplets cannot be neglected at the initial
stage of coalescence. Then, the initial stage would be better described as
inertially limited viscous (ILV) regime and a linear scaling with time for the
bridge radius has been proposed, which, according to numerical simulations, is
only realised when the coalescing drops are initially separated by a finite
distance.20 In the case of miscible and immiscible droplets, a similar viscous
dominated regime has been suggested, but immiscible droplets seem to develop a
bridge slower.56

All-atom molecular dynamics simulations of pure two-dimensional (cylindrical)
water droplets have found that multiple precursor bridges develop via thermal
fluctuations at the droplet’s surface, which initially connect the droplets
and then grow, thus identifying a thermal regime at the onset of coalescence.9
After a certain threshold, when the radius of the growing bridge becomes
larger than a thermal length scale,
$l_{T}\approx{\left(k_{B}T/\gamma\right)}^{1/4}R_{0}^{/1/2}$, the hydrodynamic
regime is recovered, continuum models can describe the process, and one
expects VR, IR, or ILV scaling. $k_{B}$ is Boltzmann’s constant and $T$ the
temperature. Since surface tension appears in the denominator, the addition of
surfactant, which lowers the surface tension, is expected to lead to the
increase of the thermal length, $l_{T}$, towards an upper limit that is
attained above the CAC. Moreover, these all-atom MD simulations have found
that, in the TR, the bridge radius expands linearly in time with a velocity
much faster than the viscous–capillary regime due to the thermal, molecular
‘jumps’ at the droplets’ surface where multiple bridges are formed.9 The
initial thermal regime has also been recently observed in the case of coarse-
grained models for pure and surfactant-laden three-dimensional (spherical)
droplets, without the formation of multiple thermal bridges but a single
pinching point instead.45

The above considerations indicate that coalescence is still under intensive
debate even in the case of pure water droplets. Moreover, only a few studies
have dealt with droplet coalescence in the presence of surfactant.40, 43, 57,
41, 42, 44, 58, 59, 60, 38, 36, 61, 37, 62, 63, 64, 65, 66 While it is
generally expected that surfactant would decrease the surface tension of the
droplets and a delay in the process would be potentially forecast,39, 61
comprehensive understanding is currently lacking, which calls for a systematic
investigation of different surfactants for a range of concentrations. It has
been experimentally shown that the presence of surfactant would locally reduce
the surface tension on the droplet surface and a nonuniform surfactant
concentration would lead to surface tension gradients that would eventually
lead to surface flow (Marangoni flow) and a rearrangement of the surfactant
molecules to counteract the gradient would delay the coalescence process.36
One might also attempt to see whether the diffusion and adsorption of
surfactant at the droplets’ LG interface would further affect the coalescence
process, for example, by influencing the flow field of the water molecules,45
especially when surfactant concentration is above the CAC and surfactant
aggregates are present in large amounts within the liquid bulk. MD simulations
based on a high-fidelity coarse-grained force-field are capable of addressing
these and other questions.

##  3 Model and Methods

An important motivation for choosing the force-field for the problem at hand
is the ability to simulate relatively large droplets that could allow for the
investigation of the surfactant mass transport mechanism with MD. Based on the
fact that coarse-grained force-fields would be a natural choice in this case,
and our success in previous investigations on the mass transport mechanism of
surfactant in the context of the superspreading phenomenon, 67, 5, 68, 4, 69
we have embarked on carrying out our studies here by using MD simulations
based on the SAFT (statistical associating fluid theory) force-field.70, 71,
72, 73, 74, 75 More specifically, a force-field based on the SAFT-$\gamma$ Mie
theory76 is used, which can accurately reproduce relevant key properties of
water–surfactant systems, such as their phase behaviour and surface
tension.67, 5, 77, 78, 79

In the case of the SAFT force-field, interactions between different coarse-
grained (CG) beads within a distance smaller than $r_{\rm c}$ are described
via the Mie potential, which is mathematically expressed as

$U(r_{\rm ij})=C\epsilon_{\rm ij}\left[\left({\frac{\sigma_{\rm ij}}{r_{\rm
ij}}}\right)^{\lambda_{\rm ij}^{\rm r}}-\left({\frac{\sigma_{\rm ij}}{r_{\rm
ij}}}\right)^{\lambda_{\rm ij}^{\rm a}}\right],r_{\rm ij}\leq r_{\rm c},$ (1)

where

$C=\left(\frac{\lambda_{\rm ij}^{\rm r}}{\lambda_{\rm ij}^{\rm r}-\lambda_{\rm
ij}^{\rm a}}\right){\left(\frac{\lambda_{\rm ij}^{\rm r}}{\lambda_{\rm
ij}^{\rm a}}\right)}^{\frac{\lambda_{\rm ij}^{\rm a}}{\lambda_{\rm ij}^{\rm
r}-\lambda_{\rm ij}^{\rm a}}}.$ (2)

i and j are the bead types, $\sigma_{\rm ij}$ indicates the effective bead
size, and $\epsilon_{\rm ij}$ is the interaction strength between any beads of
type i and j. $\lambda_{\rm ij}^{a}=6$ and $\lambda_{\rm ij}^{r}$ are Mie
potential parameters, while $r_{\rm ij}$ is the distance between two CG beads.
Units are chosen for the length, $\sigma$, energy, $\epsilon$, mass, $m$, and
time $\tau$, which in real units would correspond to: $\sigma=0.43635$ nm,
$\epsilon/k_{B}=492$ K, $m=44.0521$ amu and
$\tau=\sigma{(m/{\epsilon})}^{0.5}=1.4062$ ps. All simulations are carried out
in the NVT ensemble by using the Nosé–Hoover thermostat as implemented in the
LAMMPS package 80 with an integration time-step $\delta
t=0.005\leavevmode\nobreak\ \tau$. Our simulations took place at room
temperature ($T=25^{\leavevmode\nobreak\ \circ}$C), which in the simulation
units corresponds to $T=0.6057\leavevmode\nobreak\ \epsilon/k_{B}$. Finally, a
universal cutoff for all nonbonded (Mie) interactions is set to
$r_{c}=4.583\leavevmode\nobreak\ \sigma$.

Fig. 2:  Coarse-grained representation of water and surfactant molecules.
Surfactant’s hydrophobic beads are shown in red, while the hydrophilic parts
of the surfactant are in yellow. A cyan bead corresponds to two water
molecules (a). C10E8 (b), C10E4 (c), and Silwet-L77 (d) CG surfactant models.

We have considered surfactants of type CiEj, such as C10E8 and C10E4 (Figs 2b,
c) and a trisiloxane surfactant.4, 68, 69, 67, 48, 49, 50, 51 In the case of
CiEj surfactants, a hydrophobic alkane CG ‘C’ bead represents a $\rm-
CH_{2}-CH_{2}-CH_{2}-$ group of atoms, while a hydrophilic CG ‘EO’ bead
represents an oxyethylene group $\rm-CH_{2}-O-CH_{2}$. Silwet-L77 (Fig. 2d) is
also considered as a trisiloxane surfactant with the same number of beads as
in the case of C10E8, but different hydrophobic chemical units and
architecture, where M type beads represent a chemical group
$\rm(-CH_{3})_{3}-Si-O_{1/2}$ and D type $\rm O_{1/2}-(CH_{3})_{2}-Si-
O_{1/2}$. Finally, a water CG ‘W’ bead corresponds to two water molecules
(Fig. 2a). The nonbonded interaction parameters between the above chemical
groups, which can be used in practice to simulate a wide range of surfactants
with different molecular architecture are reported in Table 1, while the mass
of each CG bead is documented in Table 2.

Table 1:  Summary of Mie interaction parameters (Eq. 1). $\lambda_{\rm ij}^{\rm a}=6$.  i–j | $\sigma_{\rm ij}\leavevmode\nobreak\ [\sigma]$ | $\epsilon_{\rm ij}\leavevmode\nobreak\ [\epsilon/k_{B}]$ | $\lambda_{\rm ij}^{\rm r}$
---|---|---|---
W–W | 0.8584 | 0.8129 | 8.00
W–C | 0.9292 | 0.5081 | 10.75
W–EO | 0.8946 | 0.9756 | 11.94
W–M | 1.0491 | 0.8132 | 13.72
W–D | 0.9643 | 0.6311 | 10.38
C–C | 1.0000 | 0.7000 | 15.00
C–EO | 0.9653 | 0.7154 | 16.86
M–M | 1.2398 | 0.8998 | 26.00
M–D | 1.1550 | 0.7114 | 18.83
M–EO | 1.0853 | 0.8262 | 22.18
D–D | 1.0702 | 0.5081 | 13.90
D–EO | 1.0004 | 0.6355 | 16.21
EO–EO | 0.9307 | 0.8067 | 19.00
Table 2: Mass of the CG beads. Bead Type | Mass [m]
---|---
W | 0.8179
C | 0.9552
EO | 1.0000
M | 1.8588
D | 1.6833

To tether beads together in the case of surfactant chains, a bond potential is
required, which in the case of this model is harmonic, i.e.,

$V(r_{\rm ij})=0.5k(r_{\rm ij}-\sigma_{\rm ij})^{2}$ (3)

where the harmonic constant $k=295.33$ $\epsilon/\sigma^{2}$. Moreover, EO
beads experience a harmonic angle potential,

$V_{\theta}(\theta_{\rm ijk})=0.5k_{\theta}(\theta_{\rm ijk}-\theta_{0})^{2},$
(4)

where $\theta_{\rm ijk}$ is the angle between consecutive beads i, j and k
(here, i, j, k indicate the order of EO beads instead of bead type),
$k_{\theta}=4.32$ $\epsilon/$rad2, and $\theta_{0}=2.75$ rad is the
equilibrium angle. Further discussion on the model can be found in previous
studies. 68, 67, 69

To prepare the initial configuration of each system, individual droplets were
first equilibrated in the NVT ensemble. The total number of beads in the
simulations was $10^{5}$ per initial droplet, with approximately 5%
evaporation into the gas. Droplet diameters were $\sim 53\sigma$, which is
about 23 nm, similar to that of several previous studies.9, 45 Careful
consideration was given not only to observing the energy of the system, but,
also, making sure that the distribution of surfactant clusters has reached a
dynamic equilibrium and that each of them was able to diffuse a distance many
times its size. After equilibration of the individual droplets, the volume of
the simulation box was doubled and the two droplets (and the surrounding gas)
were placed next to each other as shown in Fig. 1a. In this way, roughly the
same thermodynamic conditions are maintained and further evaporation that
would reduce the number of water molecules of the droplets is avoided. The
final size of the simulation box is also chosen large enough to avoid the
interaction of mirror images of the droplets over the periodic boundary
conditions. Figure 1 illustrates typical snapshots at different stages during
coalescence for cases below and above CAC. For our study, we have considered a
range of different surfactant concentrations up to about 6$\times$CAC, which
cover the relevant span of phenomena. A summary of the mean values of various
properties for our systems is given in Table 3. We can see that increasing
surfactant concentration slightly increases the size of the droplet. Also,
note that the CAC in the case of Silwet-L77 in terms of wt% is almost double
that of C10E4 and C10E8 surfactants.

Table 3:  Properties of individual droplets (equilibrium)

Concentration ($wt\%$) |  Diameter ($\sigma$) | Water Beads* | # Molecules
---|---|---|---
C10E4 |  |  |
6.25 | 53.08 | 90466.09 | 714
12.37 | 53.02 | 85496.45 | 1429
24.18 | 53.62 | 75572.27 | 2857
35.48 | 54.14 | 65746.73 | 4286
46.02 | 54.63 | 55966.55 | 5714
CAC $\approx$ 7.5 wt% |  |  |
C10E8 |  |  |
6.25 | 52.74 | 90519.09 | 455
12.37 | 52.98 | 85500.09 | 909
24.18 | 53.40 | 75722.36 | 1818
35.48 | 53.91 | 65862.36 | 2727
46.02 | 54.33 | 56488.91 | 3636
CAC $\approx$ 7.5 wt% |  |  |
Silwet-L77 |  |  |
7.6 | 52.62 | 90438.45 | 455
14.8 | 53.01 | 85574.64 | 909
28.2 | 53.89 | 75786.36 | 1818
40.3 | 54.78 | 66312.82 | 2727
51.2 | 55.71 | 56787.91 | 3636
CAC $\approx$ 16.23 wt% |  |  |

* Indicates the average number of water beads.

Fig. 3:  Bridge analysis and domains for analysing the mass transport
mechanism. Six different regions for the position of the molecules within the
droplets are considered, namely, bridge bulk and surface, left and right
surface and bulk. (a) A snapshot obtained by the MD simulation with a clearly
developed bridge length. (b) Analysis of an $X-Z$ projection for identifying
the bridge based on the method described in Sec. 3. (c) Different regions of
the droplets considered for analysing the mass transport mechanism of
surfactants between these regions (see ESI for further details).

To obtain reliable estimates of the bridge growth dynamics and the mass
transport mechanism, snapshots of the system are made frequently (every 250 MD
time steps) for the initial $4\times 10^{5}$ MD time steps, and cluster
analysis is performed to identify the beads belonging to the liquid phase
(droplets), which is used for our properties analysis. The bridge region is
chosen to be a slab with a width in the $X$ (approach) direction that is
calculated for each configuration. In particular, the left and right limits of
the slab are determined by analysing the grid points on the $X-Z$ plane at
position $Y=0$ according to the procedure shown in Fig. 3b. We fit a circle
around each droplet and note the surface grid positions at the central $X=0$
position, shown by the red points in Fig. 3b. Horizontal lines are drawn in
the $X$ direction passing through these red points to touch the fitted
circles, thus defining the rectangle in green. The vertical sides of the
rectangle give the limits of the bridge slab in the $X$ direction, and its
width. All molecules with centres having $X$ coordinates inside these limits
are labelled as belonging to the bridge in a given snapshot. On the other
hand, the bridge radius $b$ shown in Figs 1c and 3 is calculated using the
distances between extrema of the positions of the beads belonging to the grids
located at $X=0$. That is this distance is first calculated separately for the
$Z$ coordinate to give a distance $2b_{Z}$, and then for the $Y$ coordinate to
give $2b_{Y}$. The final bridge radius estimate is then given by
$b=(b_{Z}+b_{Y})/2$.

We also describe the mass transport mechanism by tracking each surfactant
molecule during the coalescence and then identifying the probability of
adsorption of each molecule to different regions in the liquid phase, namely
the left and right LG surfaces and the bulk, and the surface and interior of
the bridge (Fig. 3b, c). Further discussion and details on the calculation of
probabilities related to the mass transport mechanism can be found in the
Electronic Supplementary Information (ESI). Finally, we have calculated the
density profiles of the water and surfactant molecules during coalescence, the
flow field at different times, as well as the approach distance and velocity
of approach of the droplets and their asphericity, for which further details
are discussed in the ESI.

Fig. 4:  Pinching of droplets with C10E4 surfactant in the TR for
concentration above (a, 35.4 wt%) and below (b, 6.25 wt%) CAC. For each case,
a magnified view of the initial bridge formation is shown, as well as the
distribution of surfactant (hydrophobic beads in red and hydrophilic in yellow
colour) and water molecules (cyan colour). (a) $t=t_{\rm
c}+8.75\leavevmode\nobreak\ \tau$, (b) $t=t_{\rm c}+8.75\leavevmode\nobreak\
\tau$,  Fig. 5:  Density profiles of surfactant beads, C10E4 (46.2 wt%) (upper
row) and water beads (lower row) on the mid plane (width of
$6\leavevmode\nobreak\ \sigma$) at different stages after the time $t_{\rm c}$
at which a permanent contact between the droplets is established (a) $t=t_{\rm
c}\leavevmode\nobreak\ \tau$, (b) $t=t_{\rm c}+71.25\leavevmode\nobreak\
\tau$, (c) $t=t_{\rm c}+98.75\leavevmode\nobreak\ \tau$.

##  4 Results and Discussion

Fig. 6:  Mass transport mechanism of surfactant (C10E8) during the coalescence
process, for concentrations above (upper panel) and below (lower panel) CAC.
The size of the arrows reflects the probabilities associated with surfactant
transport to the different droplet areas (see Table S4 of the ESI for further
details). Above the CAC (a–c, 46.2 wt%) snapshots were obtained at times (a)
$t_{\rm c}+22.5\leavevmode\nobreak\ \tau$, (b) $t_{\rm
c}+172.5\leavevmode\nobreak\ \tau$, (c) $t_{\rm c}+480.0\leavevmode\nobreak\
\tau$ while below the CAC (d–e, 6.2 wt%) times shown are (d) $t_{\rm
c}+25.0\leavevmode\nobreak\ \tau$, (e) $t_{\rm c}+66.3\leavevmode\nobreak\
\tau$, (f) $t_{\rm c}+116.3\leavevmode\nobreak\ \tau$, soon after the end of
the thermal regime (a,d), the development of the bridge and the formation of
new aggregates (b) or surfactant monomers remaining in the bridge region (e),
and the full development of the bridge (c, f). Magnified views of the bridge
region and its cross-section (showing only surfactant hydrophobic beads in the
bridge region, red) are attached above and to the right of the snapshots,
respectively. Figure S4 in the ESI illustrates results for the C10E4
surfactant.

###  4.1 Coalescence mechanism

The initial stage of droplets’ coalescence manifests itself via their pinching
as illustrated by the snapshots of Fig. 4 and movies in the ESI. In this early
stage of the phenomenon, a previous study on coalescence with pure water
droplets has shown via all-atom MD that multiple bridges form on the surface
of the droplets and the overall radius of the affected region grows linearly
over time.9 Here, we rather see that the droplet pinching initially involves
only the surfactant molecules at the droplets’ surfaces, particularly when the
surfactant concentration is higher. Here and in those previous results the
size of the droplets are about the same. Fig. 4 illustrates the single
pinching of the droplets based on our model, whereby hydrophobic parts of the
surfactants come together forming aggregates, while water molecules remain
further from the pinching point (Fig. 4). By inspecting the density profile of
the droplets during coalescence and considering a cross-section on the $X-Y$
plane passing through the centre of mass of the droplets for high-
concentration droplets, we can better highlight the aggregates at the cross-
section (Fig. 5). One can observe that the water density between the droplets
is negligible (the grid used for the calculation of the density profile does
not resolve the bead size), in contrast to the density of surfactant.
Therefore, the latter is solely responsible for the droplet pinching. At later
times and when the bridge has developed past the very initial pinching stage,
water molecules become part of the bridge (Fig. 5b). However, the bridge
region is still dominated by the presence of surfactant, with surfactant
density values similar to those inside aggregates. It is notable that
surfactant continues to be present in the bridge in significant amounts even
at later stages of coalescence, when the bridge has been almost fully
developed (Fig. 5c). The reasons for this will become more apparent when
further details on the surfactant mass transport are unveiled.

After the initial droplet pinching (Fig 4), the contact between the droplet
persists and the growth of the bridge, which manifests by the large change in
curvature at the surface, takes place as shown in Fig. 6. Moreover, we
analysed the surfactant mass transport mechanism during coalescence, and the
main adsorption processes of surfactant were monitored. These processes are
described by the probabilities of surfactant remaining at a particular region
or moving between the different regions shown in Fig. 3c during coalescence.
These are documented in the ESI for all surfactant types and concentrations
considered in our study. On the basis of these probabilities we have indicated
by arrows the main surfactant movements (Fig. 6). In particular, the
concentration of surfactant is high at the initial contact of the droplets due
to the preexisting surfactant at the droplet surfaces (Fig. 6a, d), which
becomes trapped in between the droplets. As the bridge gradually grows, we
observe that most of the surfactant, which was initially on the surfaces of
the droplets, preferably moves towards the surface of the bridge, which is
energetically more favourable. However, the bridge surface has limited place
for accommodating all surfactant from the initial LG surfaces of the droplets,
since in this area between the droplets the surface excess concentration is
doubled at the initial approach. As a result, some surfactant remains in the
bulk and forms new aggregates, which is most clearly seen in the case of
higher surfactant concentrations (Fig. 6b, c and Figs. S1–S3 and movie in the
ESI). In contrast, when the concentration is lower than CAC, the bridge
surface is able to accommodate surfactant molecules that existed on the LG
surfaces before the coalescence (see Tables S1–S3 for detailed surfactant
counts).

Fig. 7:  Density profiles of a $X-Y$ cross-section ($X$ and $Y$ are the
coordinates of the centre point of each grid voxel) passing through the centre
of mass of the droplets with (a) C10E4 (46.2 wt%) $t=t_{\rm
c}+288.75\leavevmode\nobreak\ \tau$ and (b) Silwet-L77 (51.2 wt%), $t=t_{\rm
c}+437.50\leavevmode\nobreak\ \tau$. Upper panels show the number density of
the surfactant beads, lower panels of the water.

We observe that most of the surfactant activity takes place in the bridge area
and that the dominant process is transport towards the bridge surface. We do
not see any clear evidence of Marangoni flow in the studied system, possibly
due to there not being enough space or time for it to develop. By analysing
our data, we have verified that the qualitative features of the mass transport
mechanism are independent of the surfactant type or concentration, while the
largest quantitative differences were seen in the case of Silwet-L77. In
particular, Silwet-L77 tends to form a larger number of aggregates but their
overall density is lower than that of the C10E4 and C10E8 surfactants, which
also implies that higher quantities of water molecules are found among
surfactant in the case of droplets with Silwet-L77. To illustrate these
effects, the density profiles of droplets with either C10E4 or Silwet-L77 with
high surfactant concentration are plotted in Fig. 7. Finally, in terms of the
mass transport (Table S4 in the ESI contains details), we find that Silwet-L77
have a higher tendency to move towards the surface of the bridge in comparison
with C10E4 and C10E8 surfactants, which is in line with previous observations
in systems with surfactant-laden droplets.69 In turn, this results in the
formation of fewer aggregates in the bridge in the case of droplets with
Silwet-L77 surfactant.

Fig. 8:  Flow field of water inside the droplets for the case without
surfactant (a–c), with C10E4 (46.2 wt%, d–f), and with Silwet-L77 (51.2 wt%,
g–i) surfactant. Cross-sections through the centre of mass are shown. Red
reflects the intensity of flow motion (only water) toward the bridge $X$
component of velocity $\varv$, blue away from the bridge. The averaged data of
red and blue grids indicating the total flow towards the bridge and away of it
in the droplets and in the bridge region are reported in the ESI and discussed
further in the main text. Note that the white space between the water areas
(e.g. in the bridge) includes surfactant aggregates and surfactant on the
surface. The snapshots were obtained at times (a, d, g) $t_{\rm c}$ (b)
$t_{\rm c}+90.75\leavevmode\nobreak\ \tau$, (c) $t_{\rm
c}+275.25\leavevmode\nobreak\ \tau$, (e) $t_{\rm c}+183.75\leavevmode\nobreak\
\tau$, (f) $t_{\rm c}+2446\leavevmode\nobreak\ \tau$, (h) $t_{\rm
c}+425\leavevmode\nobreak\ \tau$, and (i) $t_{\rm
c}+2437.5\leavevmode\nobreak\ \tau$.

###  4.2 Water flow

In a previous study,45 we have shown that the coalescence of water droplets is
characterised by intense internal flow variations during the thermal regime,
which attenuate with the increase of surfactant concentration. Figure 8
illustrates the flow of water in the case of water droplets and those with
surfactant after the initial thermal regime and at the initial stages of the
power-law (IR) regime. Firstly, as expected, the data suggest that there is an
increased water flow toward the bridge (red colour) throughout the coalescence
process. This is particularly visible in the bulk in the IR phase, as can be
better seen in data included in Tables S6–S8 in the ESI. These are separately
averaged for the red and the blue grids and concern the pure water and
surfactant-laden droplets at high concentration. Moreover, surfactant
attenuates the free flow of water molecules toward the bridge, because it
reduces the surface tension at the droplets’ LG surface and moreover forms
aggregates in the bulk that hinder the flow directly. The formation of
aggregates at the pinching point is also seen, which manifests by the empty
spots in the flow field of Fig. 8. These observations are valid for all of the
different surfactants studied here and for the whole range of concentrations
above CAC.

Fig. 9:  Bridge length $b$ as a function of time from the first contact of the
droplets at $t_{c}$ for different concentrations of C10E4 surfactant, as
indicated. Power-law fits $\sim t^{\beta}$ are shown, labeled “IR”. $l_{s}$ is
the thermal length in the case of surfactant-laden droplets above CAC, while
$l_{w}$ is the thermal length in the case of pure water droplets.

###  4.3 Bridge Growth

The growth rate of the bridge radius, $b$, is a key parameter that
characterises the dynamics of the coalescence process. From a fluid dynamics
perspective, the VR and IR have mainly been proposed, where the former
manifests itself through a linear dependence on time, $b\sim t$, while the
latter is described by a slower power-law dependence, namely $b\sim
t^{\beta}$55, 52, 81, 10 with expected $\beta\approx\tfrac{1}{2}$. In
addition, all-atom molecular dynamics have identified a thermal regime at the
very early stage of coalescence, which persists over the length scale $l_{T}$
as mentioned in Sec. 1.9 Well beyond this length, the hydrodynamic description
is considered valid. In view of the importance of $b$ in describing the
dynamics of coalescence and the various scenarios discussed thus far in the
literature, we have embarked here on investigating the time evolution of $b$
for a range of different surfactants and concentrations. Figure 9 presents
results for droplets with C10E4 surfactant at different concentrations, both
below and above the CAC. In Fig. S6 in the ESI, results for other surfactants
and concentrations are shown. Overall, our model captures two different
coalescence regimes, namely an initial TR and a subsequent power-law regime
characterised by power law exponents $\beta$ in the range $0.46-0.71$, mean
$\beta=0.57$ close to those predicted by IR scaling. We observe different
exponents for different surfactant concentration, with more or less random
variation from case to case, which may be due to finite size effects and
different internal aggregate configurations. Prefactors generally decrease
with growing surfactant concentration, leading to a slowdown of growth. Hence,
as in the case of all-atom water simulation,9 our CG model is also confirming
the existence of a TR and an inertial-like power-law regime. The TR persists
over a time $\mathcal{O}(10\leavevmode\nobreak\ \tau)$ during which the length
of the bridge is well described by the length $l_{T}$ both for the pure and
the surfactant-laden droplets. In Fig. 10, we compare the bridge growth for
three different surfactants below and well above CAC to each other and pure
water. Differences in the surface tension above CAC for the different
surfactants are generally expected to be small, and our results for $b(t)$
also differ a little, which agrees with the expectation that bridge size in
this regime is determined by the thermal length $l_{T}\propto\gamma^{-1/4}$.

Fig. 10:  Comparison of bridge radius growth $b(t)$ for different surfactants,
below (a) and significantly above (b) CAC. Power-law fits $\sim t^{\beta}$ are
shown, for times counted from the first contact of the droplets at $t_{c}$.

Following the TR, we observe the transition to the power-law regime. The time
for the transition to occur is longer when droplets have a higher surfactant
concentration. The power-law regime is characterised by exponents close and
generally above $0.5$. The mean over all cases is $\beta=0.57$ and standard
deviation $0.06$, while the reference pure water case gives $\beta=0.541$. It
is also clearly seen that adding surfactant will lead to slower dynamics
overall, regardless of exponent variability, as well as a delay in the start
of the bridge growth — see ESI Fig. S5. We can also see a trend that C10E8
clearly has higher $\beta$ exponent than pure water and than C10E4 (Compare
Fig. S6 in the ESI with Fig. 9). Otherwise, the main variation lies in the
reduction of the growth prefactor as surfactant concentration grows, as can be
seen by the progressive shifting of the evolution to the right in the
logarithmic plots of Figs. 9, 10, and ESI Fig. S6. This trend is broadly
similar for all surfactants that we studied, with the addition that Silwet-L77
displayed a significantly slower growth in all respects, and in particular
lower maximum velocity (see ESI, Figs. S7–S8).

Fig. 11:  System length, $X$, in the coalescence direction for different
concentrations of C10E4. The inset shows the instantaneous velocity of
approach $u=dX/dt$. The maximum speed occurs at droplet-to-bridge width ratios
of $L_{y}/l_{br}$=1.22 ($0.0$ wt%), 1.22 ($6.2$ wt%), 1.30 ($12.3$ wt%)), 1.40
($24.2$ wt%), 1.80 ($35.4$ wt%), 1.80 ($46.2$ wt%).

###  4.4 Velocity of approach and asphericity

A complementary measure of the coalescence dynamics is obtained by monitoring
the total length of the droplets in the $X$ direction, whence the velocity of
approach $u=\dot{X}$ can be obtained.82, 25 Data for droplets with C10E4
surfactant are presented in Fig. 11, and show that the fastest approach takes
place in the middle of the process at a time about 150 $\tau$, which is
similar for different surfactants. Although the exponents $\beta$ of the
bridge radius are roughly similar for all cases, we observe that the droplets
with surfactant consistently exhibit a smaller maximum speed. Moreover, the
process is shorter in time in the case of pure water droplets. A different
behaviour appears well above CAC (e.g. four times above), when the process
becomes significantly slower, sometimes rather suddenly. A conjecture is that
this is due to increased internal rigidity of the aggregates present in the
bulk. Monitoring the ratio of droplet radius to bridge radius
$L_{y}/l_{br}\approx R_{0}/b$ one finds that maximum coalescence velocity
occurs for very broad bridges $L_{y}/l_{br}\approx 1.2$ when the surfactant
concentration is low but earlier and for much smaller bridge size
$L_{y}/l_{br}\approx 1.8$ when well above CAC. This would be consistent with
the conjecture that a major slowdown occurs once significant internal
rearrangements, as necessitated by a broad bridge in high concentrations,
become necessary. Further results for different surfactants and comparative
plots are shown in the ESI (Figs S7 and S8).

Monitoring the asphericity of the two coalescing droplets (Figs S9 and S10 in
the ESI), as a complimentary measure to monitor the lateral changes in the
dimensions of the system, we can observe that this generally follows the
behaviour of the velocity of approach. Pure water droplets complete the
coalescing process obtaining the final spherical shape much faster than in the
case of droplets with surfactant. As with approach velocity and system size,
we see a change in behaviour when surfactant concentration increases, but the
change is gradual and does not occur at the CAC but well above.

##  5 Conclusions

In this study, we have analysed various macroscopic and microscopic properties
of droplets laden with surfactant and unveiled the surfactant mass transport
mechanism during the process for different surfactants and a range of
concentrations. We have demonstrated that the underlying mechanisms and flow
patterns of coalescence are universal, qualitatively and often quantitatively
independent of the type of surfactant. Particular differences between
Silwet-L77 and CiEj surfactants have been noted in the bridge dynamics for
concentrations below the CAC. We have also observed that water molecules are
not part of the initial pinching process of coalescence at larger (above CAC)
surfactant concentration, which is driven by surfactant aggregation at the
droplet’s surface.

Other features include the engulfment of part of the initial contact surface
inside the forming bridge, the existence of an initial thermal regime,
followed by a later power-law growth of the bridge radius with power exponents
close to $\frac{1}{2}$, as expected in the inertial regime. We saw no evidence
of an intermediate viscous regime. One possibility is that the droplets,
despite being large computationally, were still too small for the VR to emerge
out of the thermal regime before the IR is activated. Coalescence also
universally becomes slower as surfactant concentration grows, and we see
evidence of the appearance of a further slowdown with different qualities for
several times the critical concentration, using several different indicators.
We conjecture that this is due to the appearance of greater internal stiffness
caused by closely-packed surfactant aggregates in the bulk. The range of
concentrations studied here reflects the need in practical applications where
surfactant concentration is usually above the CAC.83 We anticipate that our
results demonstrate the mechanisms of a fundamental process in nature and
technological applications, which remain universal.

## Conflicts of interest

There are no conflicts to declare

## Acknowledgements

This research has been supported by the National Science Centre, Poland, under
grant No. 2019/34/E/ST3/00232. We gratefully acknowledge Polish high-
performance computing infrastructure PLGrid (HPC Centers: ACK Cyfronet AGH)
for providing computer facilities and support within computational grant no.
PLG/2022/015747.

## Notes and references

  * Feng _et al._ 2015 S. Feng, L. Yi, L. Zhao-Miao, C. Ren-Tuo and W. Gui-Ren, _Chinese J. Anal. Chem._ , 2015, 43, 1942–1954.
  * Baret 2012 J.-C. Baret, _Lab Chip_ , 2012, 12, 422–433.
  * Ristenpart _et al._ 2006 W. Ristenpart, P. McCalla, R. Roy and H. Stone, _Phys. Rev. Lett._ , 2006, 97, 064501.
  * Theodorakis _et al._ 2014 P. E. Theodorakis, E. A. Müller, R. V. Craster and O. K. Matar, _Curr. Opin. Colloid Interface Sci._ , 2014, 19, 283–289.
  * Theodorakis _et al._ 2019 P. E. Theodorakis, E. R. Smith, R. V. Craster, E. A. Müller and O. K. Matar, _Fluids_ , 2019, 4, 176.
  * Paulsen _et al._ 2014 J. D. Paulsen, R. Carmigniani, A. Kannan, J. C. Burton and S. R. Nagel, _Nat. Commun._ , 2014, 5, 3182.
  * Yoon _et al._ 2007 Y. Yoon, F. Baldessari, H. D. Ceniceros and L. G. Leal, _Phys. Fluids_ , 2007, 19, 102102.
  * Khodabocus _et al._ 2018 M. I. Khodabocus, M. Sellier and V. Nock, _Adv. Math. Phys._ , 2018, 2018, .
  * Perumanath _et al._ 2019 S. Perumanath, M. K. Borg, M. V. Chubynsky, J. E. Sprittles and J. M. Reese, _Phys. Rev. Lett._ , 2019, 122, 104501.
  * Eggers _et al._ 1999 J. Eggers, J. R. Lister and H. A. Stone, _J. Fluid Mech._ , 1999, 401, 293–310.
  * Aarts _et al._ 2005 D. G. Aarts, H. N. Lekkerkerker, H. Guo, G. H. Wegdam and D. Bonn, _Phys. Rev. Lett._ , 2005, 95, 164503.
  * Sprittles and Shikhmurzaev 2012 J. Sprittles and Y. Shikhmurzaev, _Phys. Fluids_ , 2012, 24, 122105\.
  * Dudek _et al._ 2020 M. Dudek, D. Fernandes, E. Helno Herø and G. Øye, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2020, 586, 124265.
  * Rahman _et al._ 2019 M. M. Rahman, W. Lee, A. Iyer and S. J. Williams, _Phys. Fluids_ , 2019, 31, 012104.
  * Berry and Dagastine 2017 J. D. Berry and R. R. Dagastine, _J. Colloid Interface Sci._ , 2017, 487, 513–522.
  * Somwanshi _et al._ 2018 P. M. Somwanshi, K. Muralidhar and S. Khandekar, _Phys. Fluids_ , 2018, 30, 092103.
  * Kirar _et al._ 2020 P. K. Kirar, K. Alvarenga, P. Kolhe, G. Biswas and K. Chandra Sahu, _Phys. Fluids_ , 2020, 32, 052103.
  * Bayani _et al._ 2018 S. Bayani, Y. Tabe, Y. T. Kang, S. H. Lee and C. K. Choi, _J. Flow Vis. Image Process._ , 2018, 25, 191–205.
  * Brik _et al._ 2021 M. Brik, S. Harmand and I. Zaaroura, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2021, 629, 127429.
  * Anthony _et al._ 2020 C. R. Anthony, M. T. Harris and O. A. Basaran, _Phys. Rev. Fluids_ , 2020, 5, 033608.
  * Kern _et al._ 2022 V. R. Kern, T. Sæter and A. Carlson, _Phys. Rev. Fluids_ , 2022, 7, L081601.
  * Heinen _et al._ 2022 M. Heinen, M. Hoffmann, F. Diewald, S. Seckler, K. Langenbach and J. Vrabec, _Phys. Fluids_ , 2022, 34, 042006.
  * Geri _et al._ 2017 M. Geri, B. Keshavarz, G. H. McKinley and J. W. M. Bush, _J. Fluid Mech._ , 2017, 833, R3.
  * Abouelsoud and Bai 2021 M. Abouelsoud and B. Bai, _Phys. Fluids_ , 2021, 33, 063309.
  * Arbabi and Theodorakis 2023 S. Arbabi and P. E. Theodorakis, _Macromol. Theory Simul._ , 2023, 32, 2300017.
  * Dekker _et al._ 2022 P. J. Dekker, M. A. Hack, W. Tewes, C. Datt, A. Bouillant and J. H. Snoeijer, _Phys. Rev. Lett._ , 2022, 128, 028004.
  * Calvo _et al._ 2019 E. Calvo, E. de Malmazet, F. Risso and O. Masbernat, _Ind. Eng. Chem. Res._ , 2019, 58, 15573–15587.
  * Sivasankar _et al._ 2022 V. S. Sivasankar, D. R. Hines and S. Das, _Langmuir_ , 2022, 38, 14084–14096.
  * R. Otazo _et al._ 2019 M. R. Otazo, R. Ward, G. Gillies, R. S. Osborne, M. Golding and M. A. K. Williams, _Soft Matter_ , 2019, 15, 6383–6391.
  * Vannozzi 2019 C. Vannozzi, _Phys. Fluids_ , 2019, 31, 082112.
  * Yeo _et al._ 2003 L. Y. Yeo, O. K. Matar, E. S. P. de Ortiz and G. F. Hewitt, _J. Colloid Interface Sci._ , 2003, 257, 93–107.
  * Hu _et al._ 2000 Y. Hu, D. Pine and L. G. Leal, _Phys. Fluids_ , 2000, 12, 484–489.
  * Mansouri _et al._ 2014 A. Mansouri, H. Arabnejad and R. Mohan, Fluids Engineering Division Summer Meeting, 2014, p. V01AT05A006.
  * Duchemin _et al._ 2003 L. Duchemin, J. Eggers and C. Josserand, _J. Fluid Mech._ , 2003, 487, 167–178.
  * Karapetsas _et al._ 2011 G. Karapetsas, R. V. Craster and O. K. Matar, _J. Fluid Mech._ , 2011, 670, 5–37.
  * Nowak _et al._ 2016 E. Nowak, N. M. Kovalchuk, Z. Che and M. J. Simmons, _Colloids Surf. A Physicochem. Eng. Asp._ , 2016, 505, 124–131.
  * Narayan _et al._ 2020 S. Narayan, A. E. Metaxas, R. Bachnak, T. Neumiller and C. S. Dutcher, _Curr. Opin. Colloid Interface Sci._ , 2020, 50, 101385.
  * Nowak _et al._ 2017 E. Nowak, Z. Xie, N. M. Kovalchuk, O. K. Matar and M. J. Simmons, _Soft Matter_ , 2017, 13, 4616–4628.
  * Leal 2004 L. Leal, _Phys. Fluids_ , 2004, 16, 1833–1851.
  * Politova _et al._ 2017 N. I. Politova, S. Tcholakova, S. Tsibranska, N. D. Denkov and K. Muelheims, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2017, 531, 32–39.
  * Dong _et al._ 2017 T. Dong, W. H. Weheliye, P. Chausset and P. Angeli, _Phys. Fluids_ , 2017, 29, 102101.
  * Dong _et al._ 2019 T. Dong, W. H. Weheliye and P. Angeli, _Phys. Fluids_ , 2019, 31, 012106\.
  * Weheliye _et al._ 2017 W. H. Weheliye, T. Dong and P. Angeli, _Chem. Eng. Sci._ , 2017, 161, 215–227.
  * Kovalchuk _et al._ 2019 N. M. Kovalchuk, M. Reichow, T. Frommweiler, D. Vigolo and M. J. H. Simmons, _Langmuir_ , 2019, 35, 9184–9193.
  * Arbabi _et al._ 2023 S. Arbabi, P. Deuar, M. Denys, R. Bennacer, Z. Che and P. E. Theodorakis, _Phys. Fluids_ , 2023, 35, 063329.
  * Li _et al._ 2008 S. Li, X. Zhang, W. Dong and W. Wang, _Langmuir_ , 2008, 24, 9344–9353.
  * Ritacco _et al._ 2010 H. A. Ritacco, F. Ortega, R. G. Rubio, N. Ivanova and V. M. Starov, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2010, 365, 199–203.
  * Ritacco _et al._ 2010 H. A. Ritacco, V. B. Fainerman, F. Ortega, R. G. Rubio, N. Ivanova and V. M. Starov, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2010, 365, 204–209.
  * Ivanova _et al._ 2010 N. Ivanova, V. Starov, R. Rubio, H. Ritacco, N. Hilal and D. Johnson, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2010, 354, 143–148.
  * Ivanova _et al._ 2012 N. Ivanova, Z. Zhantenova and V. Starov, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2012, 413, 307–313.
  * Ivanova _et al._ 2009 N. Ivanova, V. Starov, D. Johnson, N. Hilal and R. Rubio, _Langmuir_ , 2009, 25, 3564–3570.
  * Paulsen 2013 J. D. Paulsen, _Phys. Rev. E_ , 2013, 88, 063010.
  * Gross _et al._ 2013 M. Gross, I. Steinbach, D. Raabe and F. Varnik, _Phys. Fluids_ , 2013, 25, 052101.
  * Thoroddsen _et al._ 2005 S. Thoroddsen, K. Takehara and T. Etoh, _J. Fluid Mech._ , 2005, 527, 85–114.
  * Paulsen _et al._ 2012 J. D. Paulsen, J. C. Burton, S. R. Nagel, S. Appathurai, M. T. Harris and O. A. Basaran, _Proc. Natl. Acad. Sci. U.S.A._ , 2012, 109, 6857–6861.
  * Xu _et al._ 2022 H. Xu, T. Wang and Z. Che, _J. Colloid Interface Sci._ , 2022, 628, 869–877.
  * Soligo _et al._ 2019 G. Soligo, A. Roccon and A. Soldati, _J. Comput. Phys._ , 2019, 376, 1292–1311.
  * Botti _et al._ 2022 T. C. Botti, A. Hutin, E. Quintella and M. S. Carvalho, _Soft Matter_ , 2022, 18, 1423–1434.
  * Constante-Amores _et al._ 2021 C. Constante-Amores, A. Batchvarov, L. Kahouadji, S. Shin, J. Chergui, D. Juric and O. Matar, _J. Fluid Mech._ , 2021, 925, A15.
  * Kasmaee _et al._ 2018 M. Kasmaee, F. Varaminian, P. Khadiv-Parsi and J. Saien, _J. Mol. Liq._ , 2018, 263, 31–39.
  * Jaensson and Vermant 2018 N. Jaensson and J. Vermant, _Curr. Opin. Colloid Interface Sci._ , 2018, 37, 136–150.
  * Ivanov _et al._ 1999 I. B. Ivanov, K. D. Danov and P. A. Kralchevsky, _Colloids Surf. A Physicochem. Eng._ , 1999, 152, 161–182.
  * Tcholakova _et al._ 2004 S. Tcholakova, N. D. Denkov and T. Danner, _Langmuir_ , 2004, 20, 7444–7458.
  * Langevin 2019 D. Langevin, _Curr. Opin. Colloid Interface Sci._ , 2019, 44, 23–31.
  * Velev _et al._ 1993 O. D. Velev, T. D. Gurkov and R. P. Borwankar, _J. Colloid Interface Sci._ , 1993, 159, 497–501.
  * Suja _et al._ 2018 V. C. Suja, A. Kar, W. Cates, S. M. Remmert, P. D. Savage and G. G. Fuller, _Proc. Natl. Acad. Sci. U.S.A._ , 2018, 115, 7919–7924.
  * Theodorakis _et al._ 2015 P. E. Theodorakis, E. A. Müller, R. V. Craster and O. K. Matar, _Soft Matter_ , 2015, 11, 9254–9261.
  * Theodorakis _et al._ 2015 P. E. Theodorakis, E. A. Müller, R. V. Craster and O. K. Matar, _Langmuir_ , 2015, 31, 2304–2309.
  * Theodorakis _et al._ 2019 P. E. Theodorakis, E. R. Smith and Müller, _Colloids Surf. A: Physicochem. Eng. Asp._ , 2019, 581, 123810.
  * Chapman _et al._ 1989 W. G. Chapman, K. E. Gubbins, G. Jackson and M. Radosz, _Fluid Phase Equil._ , 1989, 52, 31–38.
  * Müller and Gubbins 2001 E. A. Müller and K. E. Gubbins, _Industrial & engineering chemistry research_, 2001, 40, 2193–2211.
  * Avendaño _et al._ 2011 C. Avendaño, T. Lafitte, A. Galindo, C. S. Adjiman, G. Jackson and E. A. Müller, _J. Phys. Chem B_ , 2011, 115, 11154–11169.
  * Avendaño _et al._ 2013 C. Avendaño, T. Lafitte, A. Galindo, C. S. Adjiman, E. A. Müller and G. Jackson, _J. Phys. Chem B_ , 2013, 117, 2717–2733.
  * Sergi _et al._ 2012 D. Sergi, G. Scocchi and A. Ortona, _J. Chem. Phys._ , 2012, 137, 094904\.
  * Müller and Jackson 2014 E. A. Müller and G. Jackson, _Annu. Rev. Chem. Biomol. Eng._ , 2014, 5, 405–427.
  * Lafitte _et al._ 2013 T. Lafitte, A. Apostolakou, C. Avendaño, A. Galindo, C. S. Adjiman, E. A. Müller and G. Jackson, _J. Chem. Phys._ , 2013, 139, 154504\.
  * Lobanova 2014 O. Lobanova, _PhD thesis_ , Imperial College London, 2014.
  * Lobanova _et al._ 2016 O. Lobanova, A. Mejia, G. Jackson and E. A. Müller, _J. Chem. Thermodyn._ , 2016, 93, 320–336.
  * Morgado _et al._ 2016 P. Morgado, O. Lobanova, E. A. Müller, G. Jackson, M. Almeida and E. J. Filipe, _Mol. Phys._ , 2016, 114, 2597–2614.
  * Thompson _et al._ 2022 A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton, _Comp. Phys. Comm._ , 2022, 271, 108171.
  * Frenkel 1945 J. Frenkel, _J. phys._ , 1945, 9, 385.
  * Verdier and Brizard 2002 C. Verdier and M. Brizard, _Rheol. Acta_ , 2002, 41, 514–523.
  * Gonçalves _et al._ 2022 I. M. Gonçalves, I. Castro, F. Barbosa, V. Faustino, S. O. Catarino, A. Moita, J. M. Miranda, G. Minas, P. C. Sousa and R. Lima, _Processes_ , 2022, 10, 2698.


# MC-GEN:Multi-level Clustering for Private Synthetic Data Generation

Mingchen Li111This work was done by Mingchen Li and Di Zhuang when they were
at the University of South Florida. , Di Zhuang222Di Zhuang is with Snap Inc.,
2772 Donald Douglas Loop N, Santa Monica, CA 90405., and J. Morris Chang
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
Department of Electrical Engineering, University of South Florida, Tampa, FL
33620

###### Abstract

Nowadays, machine learning is one of the most common technology to turn raw
data into useful information in scientific and industrial processes. The
performance of the machine learning model often depends on the size of
dataset. Companies and research institutes usually share or exchange their
data to avoid data scarcity. However, sharing original datasets that contain
privacy information can cause privacy leakage. Utilizing synthetic datasets
which have similar characteristics as a substitute is one of the solutions to
avoid the privacy issue. Differential privacy provides a strong privacy
guarantee to protect the individual data records which contain sensitive
information. We propose MC-GEN, a privacy-preserving synthetic data generation
method under differential privacy guarantee for multiple classification tasks.
MC-GEN builds differentially private generative models on the multi-level
clustered data to generate synthetic datasets. Our method also reduced the
noise introduced from differential privacy to improve the utility. In
experimental evaluation, we evaluated the parameter effect of MC-GEN and
compared MC-GEN with three existing methods. Our results showed that MC-GEN
can achieve significant effectiveness under certain privacy guarantees on
multiple classification tasks.

###### keywords:

Synthetic data generation, differential privacy, machine learning

††journal: Knowledge-based Systems

##  1 Introduction

With the advance of data science, the demand for data collection has increased
dramatically over the past years. A large number of data has been collected by
organizations, such as governments, corporations, and research institutes.
Machine learning plays a critical role to learn useful information from these
collected data. The performance of the machine learning model often depends on
the amount of data. When the organization does not have sufficient data, data
sharing happens among organizations that have similar research interests. It
not only provides scalability of researches but also makes up the disadvantage
of insufficient data on undergoing researches. For instance, the research
institutes do not have sufficient data for an illness diagnostic system,
hospitals can share their data (patient records) to help them. Such data
usually carries privacy information that can cause privacy leakage of
individuals. Data sharing happens very often in machine learning works. As
such, sharing data in a privacy-preserving way is of vital importance.

Various methods have been proposed for privacy-preserving data sharing.
Synthetic data release allows the data owner to publish datasets to data users
without disclosing original data. The data users can make use of synthetic
datasets for different purposes, such as machine learning, data mining, etc.
Differential privacy (DP)[1] has emerged as a provable privacy guarantee for
data at the individual sample level. Releasing differentially private
synthetic datasets prevents individual records identification and information
inference with some background knowledge. Therefore, introducing differential
privacy into synthetic data release turns to be a good solution for privacy-
preserving data sharing in machine learning to protect individual records in
the original datasets.

Designing a powerful privacy-preserving synthetic data generation method for
machine learning purposes is of great challenge. First, a privacy-preserving
method usually introduces perturbations on data samples that hurt much on the
utility of data. Mitigating the perturbation to reach a certain level of
utility is not easy. Second, machine learning has multiple tasks, like support
vector machine, logistic regression, random forest, and k-nearest neighbor. An
effective synthetic data generation method should be applicable to different
tasks. Third, some data has a complex distribution. It is hard to form an
accurate generator base on the whole data distribution.

In recent years, several works related to private synthetic data release have
been proposed in the literature [2], [3], [4], [5], [6], [7], [8], [9]. One
kind of method is to produce the synthetic database that preserved privacy and
utility regarding query sets. These algorithm [2], [3] usually designed to
answer different query classes, data samples in the dataset is not actually
release. In [4], [5], they proposed algorithms for synthetic data release
based on the noisy histograms, which is more focusing on the categorical
feature variables. Algorithm [6], [7], [9] generate synthetic dataset under
statistical model with some preprocessing on the original datasets. However,
they only consider to generate the synthetic data based on the whole data
distribution.

Our primary motivation is to share a private synthetic dataset that has a
similar statistical characteristic as the original dataset while data owners
have no privacy concerns. We proposed MC-GEN, an approach that using multi-
level clustering (sample level and feature level) and differentially private
multivariate gaussian generative model to release synthetic datasets which
satisfy $\epsilon$-differential privacy. By applying multi-level clustering to
the generative model, MC-GEN reduced the total amount of noise on synthetic
datasets. Extensive experiments on different original datasets have been
conducted to evaluate MC-GEN, which demonstrates that synthetic datasets
generated by MC-GEN maintain the utility of classification tasks while
preserving privacy. We also compared MC-GEN with other existing methods, and
our results show that MC-GEN outperforms other existing methods in terms of
effectiveness.

The main contributions of MC-GEN are summarized as following:

  * 1.

We proposed an innovative, effective synthetic data generating method, MC-GEN,
which is applicable to multiple classification tasks.

  * 2.

We demonstrated MC-GEN engaging less differentially private noise to achieve
same level of privacy while preserving the utility of synthetic datasets.

  * 3.

We investigated the effect of the parameter in MC-GEN. Extensive comparison
experiments were conducted to show our method is more effective than other
existing synthetic data generating methods on classfication tasks.

The rest of paper is organized as follows: Section 2 introduces some related
works. Section 3 presents some preliminary of our approach. Section 4
describes our methodology and generative model about MC-GEN. Section 5
presents the experimental evaluation of MC-GEN. Section 6 makes the
conclusions.

##  2 Related Work

Moritz et al. [2] proposed an algorithm combined the multiplicative weighs
approach and exponential mechanism for differentially private database
release. It finds and improves an approximation dataset to better reflect the
original data distribution by using multiplicative weighs update rule. The
weight for each data record answer to desired queries in the approximation
dataset will scale up or down depends on its contribution to the query result.
The queries are sampled and measured by using the exponential mechanism and
Laplace mechanism, which guarantee differential privacy. This work only
consider the privacy solution to a set of linear queries.

Some of the research works applied the conditional probabilistic model on
synthetic data generation. A new privacy notion, “ Plausible Deniability” [6],
has been proposed and achieved by applying a privacy test after generating the
synthetic data. The generative model proposed in this paper is a probabilistic
model which captures the joint distribution of features based on correlation-
based feature selection (CFS)[10]. The original data transforms into synthetic
data by using conditional probabilities in a probabilistic model. This work
also proved that by adding some randomizing function to “ Plausible
Deniability” can guarantee differential privacy. The synthetic data generated
by this approach contains a portion of the original data, which may lead to
the potential privacy issue. PrivBayes [7] generated the synthetic data by
releasing a private Bayesian network. Starting from a randomly selected
feature node, it extend the network iteratively by selecting a new feature
node from the parent set using the exponential mechanism. They also applied
the Laplace mechanism on the conditional probability to achieve the private
Bayesian network. Bayesian network is a good approach for discrete data.
However, using the encoding method to represent the dataset contains many
numerical data introduced more noise to synthetic datasets.

Some works form a generative model on preprocessing data. A non-interactive
private data release method has been proposed in [8]. They projected the data
into a lower dimension and proved it is nearly Gaussian distribution. After
projected the original data, the Gaussian mixture model (GMM) is used to model
the original data based on the estimated statistical information. The
synthetic data is generated by adding differential privacy noise on GMM
parameters. However, there is some information loss while projection the data
into a lower dimension. Josep Domingo-Ferrer et al. [9] proposed the release
of a differentially private dataset based on microaggregation[11], which
reduced the noise required by differential privacy based on k$-$anonymity.
They clustered original data into clusters, and records in each cluster will
be substitute by a representative record (mean) computed by each cluster. The
synthetic data is generated by applying the Laplace mechanism on
representative records. They also mentioned the idea to consider the feature
relationship while clustering, but the paper does not come up with a detailed
methodology. Using the representative records as a seed data to generate
synthetic data may lose some variance in the original datasets, which makes
the synthetic data not accurate.

##  3 Preliminary

This section introduces background knowledge of this paper.

###  3.1 Differential Privacy [3]

Differential privacy proposed by D.work et al. [3] is one of the most popular
privacy definition which guarantee strong privacy protection on individual
samples. The formal definition of $\epsilon$\- differential privacy is defined
as below:

######  Definition 1.

(Differential Privacy)

A randomized mechanisms $M$ provides $\epsilon$ \- differential privacy if,
for all datasets $D_{1}$, $D_{2}$ that differ in a single record (neighboring
datasets) and for all possible outputs S $\subset$ Range(M), it satisfies

$\frac{Pr(M(D_{1})\in S)}{Pr(M(D_{2})\in S)}\leq exp(\epsilon)$ (1)

The privacy parameter (privacy budget) $\epsilon$ used to present the privacy
level achieved by the output of the mechanisms. If $\epsilon$ = 0, it ensure
the perfect privacy guarantee, which means $D_{1}$, $D_{2}$ can not be
distinguished by observing the output of algorithm, if $\epsilon$ = $\infty$,
there is no privacy guarantee at all. In general, differential privacy
mechanism can be achieved by adding noise to the output of a query function
with a privacy parameter between 0 and 1. The mechanism used to achieve
$\epsilon$ \- differential privacy is called differentially private sanitizer
and it always associated with the sensitivity, also known as $L_{1}$ \-
Sensitivity, defined as follows:

######  Definition 2.

($L_{1}$ \- Sensitivity) For a query function $f:D^{n}\to\mathbb{R}^{d}$ which
maps datasets to real number domain, the sensitivity of the query function f
for the neighboring datasets $D$, $D^{\prime}$ is defined as follows:

$\Delta f=\max_{D,D^{\prime}}||f(D)-f(D^{\prime})||$ (2)

Laplace mechanism achieve $\epsilon$-differential privacy by adding noise
drawn from Laplace distribution [12] to the output.

######  Definition 3.

(Laplace Mechanism) For a function $f:D^{n*d}\to\mathbb{R}^{n*d}$ which maps
datasets to real number domain, the mechanism M defined in the following
equation provides $\epsilon$-differential privacy:

$M(D)=f(D)+Laplace(\Delta f/\epsilon)$ (3)

###  3.2 Agglomerative Hierarchical Clustering

An essential part of our methodology is hierarchical clustering [13].
Hierarchical clustering is an algorithm that clusters input samples into
different clusters based on the proximity matrix of samples. The proximity
matrix contains the distance between each cluster. The agglomerative
hierarchical clustering is a bottom-up approach, starting from assigning each
data sample to its own group and merge the pairs of clusters that have the
smallest distance to move up until there is only a single cluster left.

###  3.3 Microaggregation

Microaggregation [11] is a kind of dataset anonymization algorithm which can
achieve k-anonymity. There are serval setting for microaggregation, notice,
the microaggregation mentioned here is a simple heuristic method called
maximum distance to average record (MDAV) proposed by Domingo-Ferrer et al.
[11]. MDAV clusters samples into clusters, in which each cluster contains
exactly k records, except the last one. Records in the same cluster are
supposed to be as similar as possible in terms of distance. Furthermore, each
record in the cluster will be replaced by a representative record of the
cluster to complete the data anonymization.

###  3.4 Multivariate Gaussian Generative Model

The multivariate Gaussian generative model keeps the multivariate Gaussian
distribution, which is parameterized by the mean $\mu$ and covariance matrix
$\Sigma$. Formally, the density function of multivariate gaussian distribution
is given by:

$f(x)=\frac{1}{(2\pi)^{n/2}|\Sigma|^{1/2}}exp(-\frac{1}{2}(x-\mu)^{T}\Sigma^{-1}(x-\mu))$
(4)

Data samples drawn from the Multivariate Gaussian generative model are under
Multivariate gaussian distribution.

##  4 Methodology

In this section, we present the methodology of our approach.

###  4.1 Problem Statement

Given a numerical dataset $D^{n\times d}$ ($n$ samples and $d$ features) which
contains sensitive privacy records. Date owner expects to share dataset $D$ to
an untrust party in a secure way. In this paper, we aim to use synthetic
dataset $D^{\prime}$ to substitute original dataset $D$ in data sharing to
prevent privacy leakage. The synthetic dataset $D^{\prime}$ not only maintain
some certain information in original datasets but also protected by a certain
level of privacy.

Figure 1: Synthetic Data Use Case

###  4.2 Methodology Overview

Fig. 2 illustrates the design of our approach. It includes four processes
worked collectively to generate the synthetic dataset which satisfies
differential privacy:

  * 1.

Data preprocessing: Combining independent feature sets and Microaggregation
[11] to produce data clusters.

  * 2.

Statistic extraction: Extract the representative statistical information from
each data cluster.

  * 3.

Differential privacy sanitizer: Introduce differential private noise on
extracted statistical information.

  * 4.

Data generator: Generate synthetic data sample by sample from perturbed
generative model.

Figure 2: MC-GEN Algorithm

We will explain each component in our design in the following sections.

###  4.3 Data Preprocessing

The convention way to generate data clusters only groups the data in sample
level. For example, microaggregation (MDAV) [11] cluster the data in full
feature dimension and add the differentially private noise on the
representative records. Two kinds of error that may be introduced to
jeopardize the utility of this process:

  * 1.

The false feature correlation introduced by sample level clustering. When
modeling the output clusters from sample level clustering, some clusters may
carry some correlation which not exist while looking into all data samples.
This false feature correlation may apply unnecessary constraints when modeling
the data clusters, which may lead the synthetic data into a different shape.

  * 2.

The noise variance introduced by DP mechanism. Intuitively, the less noise
introduced from differential privacy results in higher utility. Hence,
reducing the noise from DP mechanism also helps us to improve the data
utility.

To smooth these two errors, we design multi-level clustering that combines
independent feature sets (IFS) with microaggregation in our approach
consecutively. Namely, we not only cluster the data in sample level, but also
feature level. Feature level clustering help the generative model to capature
the correct correlation of features. Comparing to sample level clustering,
using multi-level clustering also reduced the total noise variance introduced
to the synthetic datasets. In the next section, we show the proof that MC-GEN
reduces the amount of noise on synthetic datasets.

####  4.3.1 Feature Level Clustering

Given a numerical dataset $D^{n\times d}$, we divide d data features into m
independent feature sets using agglomerative hierarchical clustering. A
distance function that converts Pearson correlation to distance has been
designed to form the proximity matrix in hierarchical clustering. Features
that have higher correlation should have lower distance and lower correlation
corresponding to the higher distance. This approach results in that features
in the same set are more correlated to each other and less correlated to the
features in other feature sets. However, hierarchical clustering need to
specify the number of clusters to be divided. We use Davies-Bouldin [14] index
to choose the best split from all possible number of clusters. The distance
function in the proximity matrix shown as below:

$d=2n[1-Corr(X,Y)]$ (5)

####  4.3.2 Sample Level Clustering

Based on the output of feature level clustering, we applied microaggregation
[11] on each independent feature set (IFS). The purpose of microaggregation is
to assign the homogeny samples into the same cluster, which preserves more
information from the original data. On the other hand, referring to [11],
sensitivity on each sample cluster can be potentially reduced compared to the
global sensitivity. This reduction will result in involving less noise in the
differential privacy mechanism. In other words, it enhances the data utility
under the same level of privacy guarantee.

###  4.4 Statistic Extraction and Differentially Private Sanitizer

Multi-level clustering outputs several clusters of at least k records each.
Assuming each cluster forms a multivariate Gaussian distribution, the mean
($\mu$) and covariance matrix ($\Sigma$) are computed for each cluster $c$. As
shown in Fig. 3, $\mu\\_DP$ and $\Sigma\\_DP$ can be achieved by adding
differentially private noise on $\mu$ and $\Sigma$. The privacy sanitizer adds
noise on $\mu$ and $\Sigma$ to ensure the model satisfied differential
privacy, and the generative model is built based on $\mu\\_DP$ and
$\Sigma\\_DP$. Algorithm 1 illustrates the process of adding DP noise on $\mu$
and $\Sigma$.

Figure 3: Synthetic Data Generator

Input: $\mu_{i},\Sigma_{i}$

Output: $\mu_{i}\\_DP,\Sigma_{i}\\_DP$

1

2for  _$c_{i}\in C_{j}^{m}$_ do

3 $\mu_{i},\Sigma_{i}\leftarrow c_{i}$;

4 $n_{mean}\leftarrow$Draw $d$ samples from $Lap(0,\frac{\Delta}{\lvert
c_{i}\rvert\epsilon_{m}})$;

5 $n_{covarM}\leftarrow$Draw $\frac{d^{2}+d}{2}$ samples from
$Lap(0,\frac{\Delta}{\lvert c_{i}\rvert\epsilon_{m}})$;

6 $\mu_{i}\\_DP=\mu_{i}+n_{mean}$ ;

7 $\Sigma_{i}\\_DP=\Sigma_{i}+n_{covarM}$;

8 return $\mu_{i}\\_DP,\Sigma_{i}\\_DP$

9

Algorithm 1 MC-GEN privacy sanitizer

The differentially private noise added on $\mu$ and $\Sigma$ suppose to be as
little as possible to preserve the utility of synthetic data. It refers to the
generative model captures the statistical information more precisely. Thus, it
is very critical to investigate and control the noise variance introduced by
privacy sanitizer. The following proof demonstrates the noise variance of MC-
GEN.

######  Theorem 1.

The noise variance introduced by MC-GEN on mean vector $\mu$ is
$\sum\limits_{j=1}^{j}\frac{\Delta IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$

###### Proof.

If applied independent feature sets (IFS) on data features, the dataset D has
been divided to $m$ independent feature sets and data in each feature set has
been cluster into $j$ clusters. Let $\Delta IFS_{m}$ denote the $L_{1}$ \-
sensitivity of $m$th independent feature set, $C_{j}^{m}$ denote the $j$th
cluster in $IFS_{m}$, and $d_{m}$ denote the size of $IFS_{m}$.

$N_{IFS_{m}}=\sum\limits_{j=1}^{j}\frac{\Delta IFS_{m}}{\lvert
C_{j}^{m}\rvert\epsilon_{m}}d_{m}$ (6)
$N_{IFS}=\sum\limits_{m=1}^{m}N_{IFS_{m}}=N_{IFS_{1}}+N_{IFS_{2}}+...+N_{IFS_{m}}\\\
$ (7)

where the $\epsilon$ has been be divided in portion based on the size of each
IFS:

$\left\\{\begin{aligned} \epsilon_{m}&=\frac{d_{i}}{d}\epsilon\\\
\epsilon&=\epsilon_{1}+\epsilon_{2}+...+\epsilon_{m}\end{aligned}\right.$ (8)

where $d$ is number of total feature in dataset $D$.

Based on equation (6) and equation (8), use addition commutative to rewrite
equation (7) as following:

$N_{IFS}=\sum\limits_{j=1}^{j}\sum\limits_{m=1}^{m}\frac{\Delta
IFS_{m}}{\lvert C_{j}^{m}\rvert\epsilon}d\\\ $ (9)

Since data in each IFS has been cluster by microaggregation with the same
cluster size, which means:

$C_{j}^{m}=C_{j}^{m+1}$ (10)

Substituting equation 10 into equation 9

$\displaystyle N_{IFS}$ $\displaystyle=\sum\limits_{j=1}^{j}\frac{\Delta
IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$ (11)

∎

######  Theorem 2.

The noise variance introduced by MC-GEN on covariance matrix $\Sigma$ is
$\sum\limits_{j=1}^{j}\frac{d_{m}\Delta IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$

###### Proof.

Use the same notation to describe the noise variance on covariance matrix:

$N_{IFS_{m}}=\sum\limits_{j=1}^{j}\frac{\Delta IFS_{m}}{\lvert
C_{j}^{m}\rvert\epsilon_{m}}d_{m}^{2}$ (12)

$N_{IFS}=\sum\limits_{m=1}^{m}N_{IFS_{m}}=N_{IFS_{1}}+N_{IFS_{2}}+...+N_{IFS_{m}}\\\
$ (13)

Based on equation (8) and (12), rewrite equation (13):

$N_{IFS}=\sum\limits_{j=1}^{j}\sum\limits_{m=1}^{m}\frac{\Delta
IFS_{m}d_{m}}{\lvert C_{j}^{m}\rvert\epsilon}d\\\ $ (14)

substituting equation (10) into (14), since $d_{m}$ always smaller than $d$:

$\displaystyle N_{IFS}$ $\displaystyle=\sum\limits_{j=1}^{j}\frac{d_{m}\Delta
IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$ (15)

∎

The critical point of MC-GEN is to cluster the original data on feature level
on the top of sample level. The change of noise variance introduced by IFS is
also worth to be investigated.

######  Theorem 3.

The noise variance introduced by sample level clustering on mean vector $\mu$
is $\sum\limits_{j=1}^{j}\frac{\Delta D}{\lvert C_{j}\rvert\epsilon}d$

###### Proof.

If we do not applied independent feature sets on data features, the dataset D
has been cluster into $j$ clusters. Let $\Delta D$ denote the $L_{1}$ \-
sensitivity of dataset D, and $C_{j}$ denote the $j$th cluster. Refer to [11],
the noise added on cluster mean can be written as following:

$N_{noIFS}=\sum\limits_{j=1}^{j}\frac{\Delta D}{\lvert
C_{j}\rvert\epsilon}d=\frac{\Delta D}{\lvert
C_{1}\rvert\epsilon}d+\frac{\Delta D}{\lvert
C_{2}\rvert\epsilon}d+...+\frac{\Delta D}{\lvert C_{j}\rvert\epsilon}d$ (16)

∎

######  Theorem 4.

The noise variance introduced by sample level clustering on covariance matrix
$\Sigma$ is $\sum\limits_{j=1}^{j}\frac{\Delta D}{\lvert
C_{j}\rvert\epsilon}d^{2}$

###### Proof.

If not apply IFS

$N_{noIFS}=\sum\limits_{j=1}^{j}\frac{\Delta D}{\lvert
C_{j}\rvert\epsilon}d^{2}=\frac{\Delta D}{\lvert
C_{1}\rvert\epsilon}d^{2}+\frac{\Delta D}{\lvert
C_{2}\rvert\epsilon}d^{2}+...+\frac{\Delta D}{\lvert
C_{j}\rvert\epsilon}d^{2}$ (17)

∎

Comparing $N_{IFS}$ with $N_{noIFS}$ on mean and covariance matrix:

###### Proof.

For mean vector:

$\displaystyle N_{IFS}$ $\displaystyle=\sum\limits_{j=1}^{j}\frac{\Delta
IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$ (18)
$\displaystyle=\sum\limits_{j=1}^{j}\frac{\Delta IFS_{1}+\Delta
IFS_{2}+...+\Delta IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$
$\displaystyle=\sum\limits_{j=1}^{j}\frac{\Delta D}{\lvert
C_{j}\rvert\epsilon}d=N_{noIFS}$

For convariance matrix:

$\displaystyle N_{IFS}$ $\displaystyle=\sum\limits_{j=1}^{j}\frac{d_{m}\Delta
IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$ (19)
$\displaystyle=\sum\limits_{j=1}^{j}\frac{d_{1}\Delta IFS_{1}+d_{2}\Delta
IFS_{2}+...+d_{m}\Delta IFS_{m}}{\lvert C_{j}\rvert\epsilon}d$
$\displaystyle<\sum\limits_{j=1}^{j}\frac{\Delta IFS_{1}+\Delta
IFS_{2}+...+\Delta IFS_{m}}{\lvert C_{j}\rvert\epsilon}d^{2}$
$\displaystyle=\sum\limits_{j=1}^{j}\frac{\Delta D}{\lvert
C_{j}\rvert\epsilon}d^{2}=N_{noIFS}$

∎

While applying IFS, the total differential privacy noise variance on the
synthetic dataset is reduced (remains the same on the mean vector, reduced on
the covariance matrix). Using feature level clustering not only mitigates the
false feature correlation error but also helps to reduce the noise variance.
Generally, the synthetic datasets generated with less noise will have better
utility.

###  4.5 Synthetic Data Generator

The original multivariate gaussian model is parameterized by mean ($\mu$) and
covariance matrix ($\Sigma$). The algorithm 1 outputs two parameters
$\mu\\_DP$ and $\Sigma\\_DP$ that are protected by differential privacy. Hence
the multivariate Gaussian model that is parameterized by $\mu\\_DP$ and
$\Sigma\\_DP$ is also protected by differential privacy. Depending on the
post-processing invariance of different privacy, all the synthetic data
derived from differential private multivariate Gaussian models are protected
by differential privacy as well. The synthetic data is synthesized sample by
sample from the multivariate gaussian generative model with perturbed
$\mu$($\mu\\_DP$) and perturbed $\Sigma$ ($\Sigma\\_DP$).

##  5 Experimental Evaluation

In this section, we show the experimental design and discuss the results of
our approach.

###  5.1 Experiment Setting

To evaluate the performance of the proposed method, we have implemented MC-GEN
based on JAVA 8. We generated synthetic datasets and performed experiments
under different cluster sizes ($k$) and different privacy parameters
($\epsilon$). The setting of $\epsilon$ varies from 0.1 to 1 and cluster size
k is picked proportionally (20%, 40%, 60%, 80%, 100%) based on corresponding
dataset. For each synthetic dataset, we evaluate the performance of three
classification tasks: support vector machine (SVM), logistic regression,
gradient boosting. We also compared our method with other private synthetic
data release methods. All the experiments are performed in two scenarios:

  * 1.

Original data training, synthetic data testing (Scenario 1):
The classification algorithm trains on original data and tests on the
synthetic data. For each experimental dataset, 20% of the samples are used as
seed dataset to generate the synthetic dataset, and 80% is used as original
data to train the model.

  * 2.

Synthetic data training, original data testing (Scenario 2):
The classification algorithm trains on synthetic data and tests on the
original data. For each experimental dataset, 80% of the samples are used as
seed dataset to generate the synthetic dataset, and 20% is used as original
data to test the model.

Dataset Task |  SVM |  Logistic Regression |  Gradient Boosting
---|---|---|---
Scenario 1
Diabetes |  0.80 |  0.78 |  0.75
Adult |  0.84 |  0.89 |  0.91
Phishing |  0.96 |  0.98 |  0.96
Scenario 2
Diabetes |  0.82 |  0.80 |  0.85
Adult |  0.85 |  0.86 |  0.83
Phishing |  0.96 |  0.94 |  0.94
Table 1: Baseline of Experiments

###  5.2 Experiment Dataset

We conducted three datasets in our experiments from UCI Machine Learning
Repository [15] to examine the performance for different algorithms. For
dataset that contains categorical variable, we use one hot encoding to convert
it to numerical variable. All the features in the dataset are scaled to
[-1,1]:

####  5.2.1 Diabetes

This dataset contains the diagnostic measurements of patient records. It has
768 samples and 8 features. The features include some patient information,
such as blood pressure, BMI, insulin level, and age. The objective is to
identify whether a patient has diabetes.

####  5.2.2 Adult

This dataset is extract from census bureau database. It has 48,842 samples and
14 features. The features include some information, such as age, work class,
education, sex, and capital-gain/loss. The objective is to predict the annual
salary (over 50k or not) of each person.

####  5.2.3 Phishing

This dataset is used to predict the phishing websites. It has 11,055 samples
and 30 features. The features include information related address, HTML and
JavaSrcipt, such as website forwarding, request URL, and domain registration
length.

###  5.3 Experiment Setup and Metric

####  5.3.1 Experiment Setup

In each round of experiments, we assume the data owners randomly select the
original data as seed data based on the scenario. The seed data is input to
MC-GEN to generate the synthetic dataset under different parameter
combinations ($\epsilon$ and $k$). Then, the synthetic dataset is tested on
three classification tasks to evaluate the performance while publishing the
synthetic dataset. For each scenario on each dataset, we ran 100 rounds to get
the average performance.

####  5.3.2 Experiment Metric

F1 score in equation (22) is utilized as criteria to evaluate all models.

$precision=\frac{true\ postive}{true\ positive+false\ positive}$ (20)

$recall=\frac{true\ positive}{true\ positive+false\ negative}$ (21)

$F1=2\times\frac{1}{\frac{1}{recall}+\frac{1}{precision}}=2\times\frac{precision\times
recall}{precision+recall}$ (22)

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 4: Effect of cluster size (k) on diabetes dataset in scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 5: Effect of cluster size (k) on adult dataset in scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 6: Effect of cluster size (k) on phishing dataset in scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 7: Effect of cluster size (k) on diabetes dataset in scenario 2

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 8: Effect of cluster size (k) on adult dataset in scenario 2

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 9: Effect of cluster size (k) on phishing dataset in scenario 2

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 10: Comparison with other generation methods on diabetes dataset in
scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 11: Comparison with other generation methods on adult dataset in
scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 12: Comparison with other generation methods on phishing dataset in
scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 13: Comparison with other generation methods on diabetes dataset in
scenario 1

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 14: Comparison with other generation methods on adult dataset in
scenario 2

(a) SVM results

(b) Logistic regression results

(c) Gradient boosting results

Figure 15: Comparison with other generation methods on phishing dataset in
scenario 2

###  5.4 Effectiveness Analysis

In this section, we evaluate the effectiveness of our approach on three
classification datasets with three classification algorithms. The best
performance of MC-GEN on each dataset is always very close to the baseline
(Table. 1). Baselines were trained and tested on the corresponding model using
original datasets. Each baseline takes the average performance over 50
experiments.

####  5.4.1 Effectiveness Analysis of Different Parameters

The performance of MC-GEN is affected by two parameters: privacy budget
($\epsilon$) and cluster size ($k$). To evaluate the effectiveness of single
parameter, the other parameter remains the same. As shown in Fig.4, 5, 6, 7, 8
and 9, we observe that:

  1. 1.

Privacy budget ($\epsilon$): This parameter controls the noise variance on the
synthetic datasets. We investigate it from 0.1 to 1. In DP mechanism, greater
$\epsilon$ results in smaller-scaled noise, and vice versa. The performance in
all classification tasks increased as the privacy budget $\epsilon$
increasing, since less noise would be added to the synthetic datasets. While
DP sanitizer adding noise on extracted statistical information (Fig. 3),
smaller-scaled noise make $\mu\\_DP$ and $\Sigma\\_DP$ more close to $\mu$ and
$\Sigma$. In other word, the generative model captures more accurate
statistical information. Thus, a greater $\epsilon$ would result in a better
performance.

  2. 2.

Cluster size ($k$): This parameter controls the number of data points in each
cluster. Namely, the number of data points used to build each generative
model. We investigate it from 20% to 100% based on the total number of seed
data. Intuitively, the smaller $k$ would make synthetic data generation more
precise. However, based on the results shown in Fig. 5, 6, 7, 8 and 9, small
$k$ does not always result in the best performance. For instance, Fig. 5 (b),
Fig. 7 (b) and 9 (a), (b) show that a moderate k achieve the best performance.
Fig. 6 and Fig. 9 show that using a large k achieve the best performance. As
the classification task changes, the local optimal k varies in line with the
dataset. There is a unique local optimal k for each dataset under certain
settings. This k helps MC-GEN capture the statistical representation of
dataset (clusters) in a good manner. Therefore, when the classification task
is assigned, it is worth to investigate the local optimal k value for the
dataset to improve the performance.

####  5.4.2 Comparing with Existing Methods

In this section, we compared MC-GEN with three existing private synthetic data
release methods:

  * 1.

PrivBayess [7]: PrivBayes [7] is designed on the Bayesian network with
differential privacy. Starting from a randomly selected feature node, it
extends the network iteratively by selecting a new feature node from the
parent set using the exponential mechanism. Then it applied the Laplace
mechanism on the conditional probability to achieve the private Bayesian
network. The synthetic datasets are generated by using the perturbated
Bayesian network. PrivBayes [7] are evaluated on the code provided by the
paper authors.

  * 2.

RonGauss [8]: RonGauss [8] releases the synthetic data based on Gaussian
mixture model. It builds the Gaussian mixture model on the projected data in a
lower dimension, and generates the synthetic dataset based on the Gaussian
mixture model with differential privacy noise. RonGauss [8] are evaluated on
the code provided by the paper authors.

  * 3.

NoIFS: NOIFS follows the same procedure (Fig. 2) as MC-GEN to generate the
synthetic data. The difference between NOIFS and MC-GEN is that NOIFS only
applies the sample level clustering (MDAV) in the data prepossessing.
Comparing it with MC-GEN helps us to evaluate the effectiveness of feature
level clustering.

The performance of MC-GEN is evaluated under the local optimal parameter
combination. For other methods, we follow the exact settings mentioned in the
corresponding paper. Each method is evaluated on the average of 100
experiments. Fig. 10, 11, 12, 13, 14 and 15 show the results of different
methods under the same privacy budgets. We can observe that MC-GEN outperforms
the other methods in most cases. As shown in the results, The performance of
PrivBayes and NOIFS are similar to MC-GEN on diabetes and adults datasets. MC-
GEN always outperforms NOIFS due to the reduced noise variance (proof in
section 4.4). The performance of PrivBayes varies as the experiment setting
changes. On adults dataset, PrivBayes outperforms MC-GEN (Fig. 11 (a)) and
NoIFS (Fig. 11 (b), (c) and Fig. 14 (c)). It can be derived from the features
of the datasets. Most of the features in adult datasets are categorical
features that would benefit the Bayesian network. For diabetes and phishing
datasets that have mixed feature types, MC-GEN outperforms the PrivBayes
consistently.

##  6 Conclusion

In this paper, we proposed a novel synthetic date releasing method, MC-GEN,
which targets on generating a private synthetic dataset for classification
tasks. We also demonstrate that our method involves less noise on the
synthetic data by introducing feature level clustering. In experiment
evaluation, we investigated the effect of the parameters in MC-GEN. We also
compared MC-GEN with other existing private synthetic data release methods.
The experimental results show that MC-GEN outperforms the other existing
methods in terms of utility.

## References

  * [1] C. Dwork, Differential privacy: A survey of results, in: International Conference on Theory and Applications of Models of Computation, Springer, 2008, pp. 1–19.
  * [2] M. Hardt, K. Ligett, F. McSherry, A simple and practical algorithm for differentially private data release, in: Advances in Neural Information Processing Systems, 2012, pp. 2339–2347.
  * [3] C. Dwork, Differential privacy, Encyclopedia of Cryptography and Security (2011) 338–340.
  * [4] K. Taneja, T. Xie, Diffgen: Automated regression unit-test generation, in: Proceedings of the 2008 23rd IEEE/ACM International Conference on Automated Software Engineering, IEEE Computer Society, 2008, pp. 407–410.
  * [5] D. Su, J. Cao, N. Li, M. Lyu, Privpfc: differentially private data publication for classification, The VLDB JournalThe International Journal on Very Large Data Bases 27 (2) (2018) 201–223.
  * [6] V. Bindschaedler, R. Shokri, C. A. Gunter, Plausible deniability for privacy-preserving data synthesis, Proceedings of the VLDB Endowment 10 (5) (2017) 481–492.
  * [7] J. Zhang, G. Cormode, C. M. Procopiuc, D. Srivastava, X. Xiao, Privbayes: Private data release via bayesian networks, ACM Transactions on Database Systems (TODS) 42 (4) (2017) 25.
  * [8] T. Chanyaswad, C. Liu, P. Mittal, Coupling dimensionality reduction with generative model for non-interactive private data release, arXiv preprint arXiv:1709.00054 (2017).
  * [9] J. Soria-Comas, J. Domingo-Ferrer, Differentially private data sets based on microaggregation and record perturbation, in: Modeling Decisions for Artificial Intelligence, Springer, 2017, pp. 119–131.
  * [10] M. A. Hall, Correlation-based feature selection for machine learning (1999).
  * [11] J. Domingo-Ferrer, V. Torra, Ordinal, continuous and heterogeneous k-anonymity through microaggregation, Data Mining and Knowledge Discovery 11 (2) (2005) 195–212.
  * [12] C. Dwork, F. McSherry, K. Nissim, A. Smith, Calibrating noise to sensitivity in private data analysis, in: Theory of cryptography conference, Springer, 2006, pp. 265–284.
  * [13] S. C. Johnson, Hierarchical clustering schemes, Psychometrika 32 (3) (1967) 241–254.
  * [14] D. L. Davies, D. W. Bouldin, A cluster separation measure, IEEE transactions on pattern analysis and machine intelligence (2) (1979) 224–227.
  * [15] D. Dua, C. Graff, UCI machine learning repository (2017).
URL http://archive.ics.uci.edu/ml