#  Direct numerical simulations of turbulent jets: vortex-interface-surfactant
interactions

C. R. Constante-Amores1,2,3 T. Abadie3 L. Kahouadji3 S. Shin4 J. Chergui5 D.
Juric5,6 A. A. Castrejon-Pita1 and O. K. Matar3<EMAIL_ADDRESS>1Department of Engineering Science, University of Oxford, Oxford OX1 3PJ, UK
2Mathematical Institute, University of Oxford, Oxford OX2 6GG, UK
3Department of Chemical Engineering, Imperial College London, London SW7 2AZ,
UK
4Department of Mechanical and System Design Engineering, Hongik University, KR
5 Université Paris Saclay, Centre National de la Recherche Scientifique
(CNRS), Laboratoire Interdisciplinaire des Sciences du Numérique (LISN), 91400
Orsay, France
6Department of Applied Mathematics and Theoretical Physics, University of
Cambridge, Cambridge CB3 0WA, UK

###### Abstract

We study the effect of insoluble surfactants on the spatio-temporal evolution
of turbulent jets. We use three-dimensional numerical simulations and employ
an interface-tracking/level-set method that accounts for surfactant-induced
Marangoni stresses. The present study builds on our previous work (Constante
et al., 2021, J. Fluid Mech., 922, A6) in which we examined in detail the
vortex-surface interaction in the absence of surfactants. Numerical solutions
are obtained for a wide range of Weber and elasticity numbers in which
vorticity production is generated by surface deformation and surfactant-
induced Marangoni stresses. The present work demonstrates, for the first time,
the crucial role of Marangoni stresses, brought about by surfactant
concentration gradients, in the formation of coherent, hairpin-like vortex
structures. These structures have a profound influence on the development of
the three-dimensional interfacial dynamics. We also present theoretical
expressions for the mechanisms that influence the rate of production of
circulation in the presence of surfactants for a general, three-dimensional,
two-phase flow and highlight the dominant contribution surfactant-induced
Marangoni stresses.

##  1 Introduction

The atomisation of a liquid jet has driven interest in the fluid mechanics
community because of its occurrence in both natural and industrial
applications (e.g., propellant combustion, pharmaceutical sprays, etc.). The
process results in a ‘cascade mechanism’ for fluid fragmentation (Plateau,
1873; Eggers, 1997; Marmottant & Villermaux, 2004; Constante-Amores et al.,
2020a): from the growth of linear modes through a Kelvin-Helmholtz instability
to the development of nonlinearities leading to capillary breakup events via
long filament pinch-off that can be modulated by a Rayleigh-Plateau
instability or controlled by an ‘end-pinching’ mechanism. The understanding of
the interfacial dynamics relies on the characterisation of the vortex-
interface interactions. For instance, Jarrahbashi et al. (2016), Zandian et
al. (2018, 2019) and Constante-Amores et al. (2021a) reported that their
interplay determines the interfacial dynamics for turbulent jets; Hoepffner &
Paré (2013) showed that vorticity production results in a change in the
capillary retraction of a liquid thread. Theoretically, Longuet-Higgins
(1992), Wu (1995), Lundgren & Koumoutsakos (1999) demonstrated that vorticity
production depends on the velocity field and the interfacial curvature for the
condition of zero shear stress at a free surface. Additionally, Brøns et al.
(2014) and Terrington et al. (2020, 2021) extended the previous results to
show that interfacial curvature effects, viscosity and density difference
across the interface are the only mechanisms driving vorticity production.
Recently, Fuster & Rossi (2021) also demonstrated the role of interfacial
curvature and density differences across the interface with identical
dynamical viscosity via two-dimensional, non-axisymmetric numerical studies.

We note that the studies mentioned in the foregoing involve a constant surface
tension and therefore do not support the formation of Marangoni gradients.
Liquid streams, however, are invariably contaminated with surface-active-
agents (surfactants), deliberately-placed or naturally-occurring, which give
rise to surface tension gradients, and subsequently Marangoni-induced flow
(Manikantan & Squires, 2020). While the atomisation of uncontaminated liquid
jets has received significant attention in the literature (Herrmann, 2010;
Desjardins & Pitsch, 2010; Jarrahbashi & Sirignano, 2014; Jarrahbashi et al.,
2016; Zandian et al., 2018, 2019; Constante-Amores et al., 2020b, 2021a), the
effect of surfactant on their dynamics remains far less studied. The multi-
scale nature of the flow, and the complex coupling between the surfactant
concentration fields and interfacial topology complicate its experimental
scrutiny. This can be alleviated via the use of high-fidelity simulations
which can unravel the delicate interplay among the different physical
mechanisms across the relevant scales.

Through the use of state-of-the-art imaging techniques, Kooij et al. (2018),
Sijs & Bonn (2020), and Sijs et al. (2021) showed that the presence of
surfactants influences the interfacial fragmentation during atomisation and
decreases the mean-droplet size in agreement with Ellis et al. (2001) and
Ariyapadi et al. (2004). All the previous studies, however, have not reported
the role of Marangoni stresses which the present paper will address for the
case of an insoluble surfactant. Although the presence of surfactants can also
induce both shear and dilatational surface rheological effects (Wee et al.,
2021), these effects will not be considered in this study. Nonetheless, we
will use transient numerical simulations to demonstrate that the Marangoni
stresses influence the production of vorticity near the interface, and modify
the interface-vortex interactions and the three-dimensional destabilisation of
the jet. In order to focus on the role of Marangoni stresses in the jet
dynamics, we will study the case of a jet of one fluid issuing into another
characterised by equal densities and viscosities.

The rest of this paper is structured as follows: in Section 2, the problem
formulation, governing dimensionless parameters, and numerical method are
introduced. Section 3 provides a discussion of the results, and concluding
remarks are given in Section 4.

##  2 Problem formulation and numerical method

Since the aim here is to shed light on the different mechanisms that influence
the production of vorticity near the interface in the presence of surfactants,
we present a general theoretical description of vorticity and circulation in a
three-dimensional control volume enclosing an interface using Lighthill’s and
Lyman’s flux definitions Terrington et al. (2021). We also provide a brief
description of the numerical technique which is used to carry out the
computations. Finally, we provide motivation for the choice of physical and
physico-chemical parameters made in the present work.

|
---|---
(a) | (b)
Figure 1: (a) Initial interfacial shape, highlighting the computational domain
of size $(5D)^{3}$ in a three-dimensional Cartesian space
$\mathbf{x}=(x,y,z)$; (b) schematic representation of the problem in the $x-y$
($z=2.5D$) plane showing the initial ($t=0$) streamwise velocity profile,
$u_{x}$, and a representation of a monolayer of an insoluble surfactant.

###  2.1 Problem formulation and numerical method

Figure 1 shows a representation of the flow configuration considered in this
study in a three-dimensional Cartesian domain $\mathbf{x}=\left(x,y,z\right)$:
a liquid segment is initialised as a cylinder of diameter $D$, with a finite
length, i.e. $5D$, in the positive $x-$(streamwise) direction. Such an
approach has been used by Desjardins & Pitsch (2010), Jarrahbashi et al.
(2016), and Zandian et al. (2018) for planar and cylindrical jets. We will
focus on the case of insoluble surfactants, which enables us to isolate the
surfactant-induced Marangoni dynamics during the atomisation of the jet. We
acknowledge, however, that experimental studies feature soluble surfactants
which are dissolved in the liquid that issues from a nozzle to form the jet
and that the sorption kinetics control the surfactant interfacial
concentration adding extra layers of complexity.

The dimensional governing equations, which can be found in the work of Shin et
al. (2018), are rendered dimensionless using the following scalings:

$\quad\tilde{\mathbf{x}}=\frac{\mathbf{x}}{D},\quad\tilde{t}=\frac{t}{t_{r}},\quad\tilde{\textbf{u}}=\frac{\textbf{u}}{U},\quad\tilde{p}=\frac{p}{\rho
U^{2}},\quad\tilde{\sigma}=\frac{\sigma}{\sigma_{s}},\quad\tilde{\Gamma}=\frac{\Gamma}{\Gamma_{\infty}},$
(1)

where, $t$, u, and $p$ stand for time, velocity, and pressure, respectively;
here, the dimensionless variables are designated using tildes. The physical
parameters correspond to the liquid density $\rho$, viscosity, $\mu$, surface
tension, $\sigma$, surfactant-free surface tension, $\sigma_{s}$, initial jet
diameter, $D$, and injection velocity, $U$. Hence, the characteristic time
scale based on the injection velocity is $t_{r}=D/U$. The interfacial
surfactant concentration, $\Gamma$, is scaled with the saturation interfacial
concentration, $\Gamma_{\infty}$.

Using the relations in Eq. (1), the dimensionless form of the continuity and
momentum equations is respectively expressed as:

$\nabla\cdot\tilde{\textbf{u}}=0,$ (2)
$\displaystyle\tilde{\rho}\left(\frac{\partial\tilde{\textbf{u}}}{\partial\tilde{t}}+\tilde{\textbf{u}}\cdot\nabla\tilde{\textbf{u}}\right)$
$\displaystyle=$
$\displaystyle-\nabla\tilde{p}+\frac{1}{\Rey}\nabla\cdot\left[\tilde{\mu}(\nabla\tilde{\textbf{u}}+\nabla\tilde{\textbf{u}}^{T})\right]$
(3)
$\displaystyle+\frac{1}{\mbox{{We}}}\int_{\tilde{A}(\tilde{t})}\left(\tilde{\sigma}\tilde{\kappa}\hat{\mathbf{s}}+\nabla_{s}\tilde{\sigma}\right)\delta\left(\tilde{\textbf{x}}-\tilde{\textbf{x}}_{{}_{f}}\right)d\tilde{A},$

where $\tilde{\kappa}$ represents the interface curvature, $\nabla_{s}$ the
surface gradient operator, and $\hat{\mathbf{s}}$ the outward-pointing unit
normal to the interface. Here, $\tilde{\textbf{x}}_{f}$ is the parametrization
of the time-dependent interface area $\tilde{A}(\tilde{t})$, where
$\delta(\tilde{\textbf{x}}-\tilde{\textbf{x}}_{f})$ is the three-dimensional
Dirac delta function. The density, $\tilde{\rho}$, and viscosity,
$\tilde{\mu}$, are given by the following expressions

$\tilde{\rho}\left(\tilde{\textbf{x}},\tilde{t}\right)=\frac{\rho_{g}}{\rho_{l}}+\left(1-\frac{\rho_{g}}{\rho_{l}}\right)H\left(\tilde{\textbf{x}},\tilde{t}\right),~{}\tilde{\mu}\left(\tilde{\textbf{x}},\tilde{t}\right)=\frac{\mu_{g}}{\mu_{l}}+\left(1-\frac{\mu_{g}}{\mu_{l}}\right)H\left(\tilde{\textbf{x}},\tilde{t}\right),$
(4)

where $H\left(\tilde{\textbf{x}},\tilde{t}\right)$ represents a smoothed
Heaviside function; this is zero in the gas phase and unity in the liquid
phase, while the subscripts $l$ and $g$ designate the individual liquid and
gas phases, respectively.

The dimensionless surfactant transport is given by:

$\frac{\partial\tilde{\Gamma}}{\partial\tilde{t}}+\nabla_{s}\cdot(\tilde{\Gamma}\tilde{\textbf{u}}_{\text{t}})=\frac{1}{Pe_{s}}\nabla^{2}_{s}\tilde{\Gamma},$
(5)

where
$\tilde{\textbf{u}}_{\text{t}}=(\tilde{\textbf{u}}_{\text{s}}\cdot\textbf{t})\textbf{t}$
is the tangential velocity vector in which $\tilde{\textbf{u}}_{\text{s}}$ is
the surface velocity and ${\mathbf{t}}$ is the unit tangent to the interface.

The scaling results in the following dimensionless groups:

$Re=\frac{\rho UD}{\mu},~{}~{}~{}\quad We=\frac{\rho
U^{2}D}{\sigma_{s}},~{}~{}~{}Pe_{s}=\frac{UD}{\mathcal{D}_{s}},~{}~{}~{}\beta_{s}=\frac{\Re\mathcal{T}\Gamma_{\infty}}{\sigma_{s}},$
(6)

where $Re$, $We$, and $Pe_{s}$ denote the Reynolds, Weber, and (interfacial)
Peclet numbers, respectively, while $\beta_{s}$ is a surfactant elasticity
number which represents a measure of the sensitivity of $\sigma$ to $\Gamma$;
here, $\Re$ is the ideal gas constant value $8.314$ J K-1 mol-1, $T$ denotes
temperature and $\mathcal{D}_{s}$ refers to the diffusion coefficient.

To describe the relation between $\tilde{\sigma}$ and $\tilde{\Gamma}$, we use
the non-linear Langmuir equation:

$\tilde{\sigma}=1+\beta_{s}\ln{(1-\tilde{\Gamma})}.$ (7)

Surface tension gradients are expressed as a function of $\tilde{\Gamma}$ as

$\nabla_{s}\tilde{\sigma}/\mbox{{We}}=-\mbox{{Ma}}/(1-\tilde{\Gamma})\nabla_{s}\tilde{\Gamma},$
(8)

where $Ma=\beta_{s}/We=\Re T\Gamma_{\infty}/\rho U^{2}D$ is a Marangoni
parameter.

The three-dimensional numerical simulations were performed by solving the two-
phase Navier-Stokes equations in the Cartesian domain
$\mathbf{x}=\left(x,y,z\right)$. A hybrid front-tracking/level-set method was
used to treat the interface where surfactant transport was resolved in the
plane of the interface (Shin et al., 2018). The simulations are initialised
with a turbulent velocity profile in the liquid jet segment (i.e.,
$u(r)=15/14~{}U(1-(r/(D/2))^{28})$ (Constante-Amores et al., 2021a). Solutions
are sought subject to Neumann boundary conditions on all variables at the
lateral boundaries, and periodic boundary conditions in the $x-$(streamwise)
direction. The computational domain is a cube with dimensions $(5D)^{3}$
globally resolved by a uniform grid of $(786)^{3}$ cells; see Appendix of
Constante-Amores et al. (2021a) for details of mesh-refinement studies and
validation of the numerical method. This method has also been widely tested
for surfactant-laden flows (Shin et al., 2018; Constante-Amores et al., 2020a,
2021b; Batchvarov et al., 2021) and the numerical simulations in this study
conserve fluid volume and surfactant mass with a relative error of less than
$10^{-3}\%$.

Next, we motivate the values of material properties by looking into the
sources for vorticity production at an interface in a three-dimensional
framework. These sources are due to differences in density (i.e., baroclinic
effect) and viscosity, surface tension forces (due to gradients of curvature
along the interface), and Marangoni stresses. Thus, to unravel the importance
of the surfactant-induced Marangoni stresses on the vortex-surface-surfactant
interactions, we focus on situations in which surface tension forces and
Marangoni stresses are the only physical mechanisms responsible for vorticity
production at the interface, i.e., the jump in material properties across the
interface is zero (Fuster & Rossi, 2021). This is a realistic assumption for
immiscible liquid-liquid systems exemplified by the silicone oil-water pairing
used by Ibarra (2017) and Ibarra et al. (2020) in their two-phase, stratified
pipe flow experiments.

The values of the dimensionless quantities are consistent with experimentally-
realisable systems and are chosen to ensure a full coupling between
surfactant-induced Marangoni stresses and interfacial diffusion, and inertia.
We set $Re=5000$ to ensure a rich dynamics (Constante-Amores et al., 2021a)
and focus on the range $50<We<1000$ to account for realistic values of
$\sigma_{s}$, i.e.
$\mathcal{O}(10^{-3})<\sigma_{s}<\mathcal{O}(10^{-1})~{}{\rm N}~{}{\rm
m}^{-1}$. The parameter $\beta_{s}$ is related to $\Gamma_{\infty}$ and
therefore the critical micelle concentration (CMC), i.e.
$\Gamma_{\infty}\sim\mathcal{O}(10^{-6})$ mol m-2 for NBD-PC
(1-palmitoyl-2-12-[(7-nitro-2-1,3-benzoxadiazol-4-yl)amino]dodecanoyl-sn-
glycero-3 -phosphocholine) (Strickland et al., 2015); thus, we have explored
the range of $0.1<\beta_{s}<0.9$ which corresponds to CMC in the range
$\mathcal{O}(10^{-7})<$ CMC $<\mathcal{O}(10^{-6})$ mol m-2, for typical
values of $\sigma_{s}$. We have set $Pe_{s}=10^{2}$ following Batchvarov et
al. (2020) and Constante-Amores et al. (2020a) who showed that the interfacial
dynamics are weakly-dependent on $Pe_{s}$ beyond this value.

###  2.2 Vorticity and circulation

Figure 2: Schematic showing a volume $V$ with a surface $\partial V$ which
encloses two fluids separated by an interface surface $I$

In order to examine the effect of the surfactant on the vorticity near the
interface, we consider a fixed, three-dimensional (3D) control volume $V$
bounded by a closed surface of area $\partial V$ with an outward-pointing unit
normal $\hat{\mathbf{n}}$ (see figure 2). This volume encloses regions of the
incompressible fluids 1 and 2, of volumes, $V_{1}$ and $V_{2}$, separated by
an interfacial surface $I$, which the intersection with $V$ defines the curve
$\partial I$. $\hat{\mathbf{s}}$ is the outward-pointing unit normal to the
surface $I$ while $\hat{\mathbf{t}}$ and $\hat{\mathbf{b}}$ are two orthogonal
unit tangent vectors to the interface. We proceed below using dimensional
variables and then apply the scalings in equation 1 to render the final
equations dimensionless.

For fluid ‘i’, it is possible to write down expressions for $\omega_{b,i}$ and
$\omega_{t,i}$, which represent the components of the vorticity
$\boldsymbol{\omega}_{i}$ in the $\hat{\mathbf{b}}$ and $\hat{\mathbf{t}}$
directions, respectively:

$\displaystyle\omega_{b,i}=(\hat{\mathbf{s}}\times\hat{\mathbf{t}})\cdot(\nabla\times\mathbf{u}_{i}),$
(9)
$\displaystyle\omega_{t,i}=(\hat{\mathbf{s}}\times\hat{\mathbf{b}})\cdot(\nabla\times\mathbf{u}_{i}),$
(10)

where $\mathbf{u}_{i}$ denotes the velocity fields. These expressions may be
recast as follows 111Using
$(\mathbf{a}\times\mathbf{b})\cdot(\mathbf{c}\times\mathbf{d})=(\mathbf{a}\cdot\mathbf{c})(\mathbf{b}\cdot\mathbf{d})-(\mathbf{a}\cdot\mathbf{d})(\mathbf{b}\cdot\mathbf{c})$,
valid for any vector $\mathbf{a}$, $\mathbf{b}$, $\mathbf{c}$, and
$\mathbf{d}$.

$\displaystyle\omega_{b,i}=\hat{\mathbf{s}}\cdot\nabla\mathbf{u}_{i}\cdot\hat{\mathbf{t}}-\hat{\mathbf{t}}\cdot\nabla\mathbf{u}_{i}\cdot\hat{\mathbf{s}},$
(11)
$\displaystyle\omega_{t,i}=\hat{\mathbf{s}}\cdot\nabla\mathbf{u}_{i}\cdot\hat{\mathbf{b}}-\hat{\mathbf{b}}\cdot\nabla\mathbf{u}_{i}\cdot\hat{\mathbf{s}}.$
(12)

In the presence of interfacial stresses arising from gradients of surface
tension $\sigma$ due to surfactant concentration gradients, the interfacial
shear stress conditions are given by

$\displaystyle\left[\left[\hat{\mathbf{t}}\cdot\mathbf{T}\cdot\hat{\mathbf{s}}\right]\right]$
$\displaystyle=$ $\displaystyle-\hat{\mathbf{t}}\cdot\nabla\sigma,$ (13)
$\displaystyle\left[\left[\hat{\mathbf{b}}\cdot\mathbf{T}\cdot\hat{\mathbf{s}}\right]\right]$
$\displaystyle=$ $\displaystyle-\hat{\mathbf{b}}\cdot\nabla\sigma,$ (14)

$[[q]]=q_{2}-q_{1}$ represents the jump across the interface of a quantity
$q$, $\mathbf{T}_{i}=-p_{i}+\mu_{i}\mathbf{D}_{i}$ is the total stress in
fluid ‘i’ in which $p_{i}$ is the pressure,
$\mathbf{D}_{i}=(\nabla\mathbf{u}_{i}+\nabla\mathbf{u^{T}}_{i})/2$ is the rate
of deformation tensor, and $\mu_{i}$ denote the viscosities, whence

$\displaystyle\left[\left[\mu\left(\hat{\mathbf{t}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}+\hat{\mathbf{s}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{t}}\right)\right]\right]=-2\hat{\mathbf{t}}\cdot\nabla\sigma,$
(15)
$\displaystyle\left[\left[\mu\left(\hat{\mathbf{b}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}+\hat{\mathbf{s}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{b}}\right)\right]\right]=-2\hat{\mathbf{b}}\cdot\nabla\sigma.$
(16)

Substitution of these results into Eqs. (11) and (12) yields

$\displaystyle\left[\left[\mu\left(\omega_{b}+2\hat{\mathbf{t}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}\right)\right]\right]=-2\hat{\mathbf{t}}\cdot\nabla\sigma,$
(17)
$\displaystyle\left[\left[\mu\left(\omega_{t}+2\hat{\mathbf{b}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}\right)\right]\right]=-2\hat{\mathbf{b}}\cdot\nabla\sigma.$
(18)

For the case $[[\mu]]=0$, which is the focus of this paper, we obtain

$\displaystyle\left[\left[w_{b}\right]\right]$ $\displaystyle=$
$\displaystyle-\frac{2}{\mu}\nabla\sigma\cdot\hat{\mathbf{t}}-2\left[\left[\hat{\mathbf{t}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}\right]\right],$
(19) $\displaystyle\left[\left[w_{t}\right]\right]$ $\displaystyle=$
$\displaystyle-\frac{2}{\mu}\nabla\sigma\cdot\hat{\mathbf{b}}-2\left[\left[\hat{\mathbf{b}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}\right]\right],$
(20)

where $\mu_{2}=\mu_{1}=\mu$. Noting that
$\hat{\mathbf{t}}\cdot\nabla=\partial/\partial s$ and
$\hat{\mathbf{b}}\cdot\nabla=\partial/\partial b$, it can be shown that

$\left[\left[\omega_{b}\right]\right]=-\frac{2}{\mu}\frac{\partial\sigma}{\partial
s}-2\left[\frac{\partial}{\partial
s}\left[\left[\mathbf{u}\cdot\hat{\mathbf{s}}\right]\right]-\kappa_{1}\left[\left[\mathbf{u}\cdot\hat{\mathbf{t}}\right]\right]\right],$
(21)
$\left[\left[\omega_{t}\right]\right]=-\frac{2}{\mu}\frac{\partial\sigma}{\partial
b}-2\left[\frac{\partial}{\partial
b}\left[\left[\mathbf{u}\cdot\hat{\mathbf{s}}\right]\right]-\kappa_{2}\left[\left[\mathbf{u}\cdot\hat{\mathbf{b}}\right]\right]\right],$
(22)

where the curvatures $\kappa_{1}$ and $\kappa_{2}$ are defined as follows

$\kappa_{1}=\hat{\mathbf{t}}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial
s},~{}~{}~{}~{}~{}\kappa_{2}=\hat{\mathbf{b}}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial
b}.$ (23)

From continuity of the normal and tangential components of the velocity at the
interface, i.e., $[[\mathbf{u}\cdot\hat{\mathbf{s}}]]=0$, and
$[[\mathbf{u}\cdot\hat{\mathbf{t}}]]=[[\mathbf{u}\cdot\hat{\mathbf{b}}]]=0$,
respectively, it is seen that the interfacial jumps in the vorticity
components are directly related to the Marangoni stresses:

$\displaystyle\left[\left[\omega_{b}\right]\right]$ $\displaystyle=$
$\displaystyle-\frac{2}{\mu}\frac{\partial\sigma}{\partial s},$ (24)
$\displaystyle\left[\left[\omega_{t}\right]\right]$ $\displaystyle=$
$\displaystyle-\frac{2}{\mu}\frac{\partial\sigma}{\partial b}.$ (25)

We now consider the circulation vector $\boldsymbol{\Omega}$ for 3D flows
given by

$\boldsymbol{\Omega}=\int_{V}\boldsymbol{\omega}dV,$ (26)

for the fixed 3D control volume $V$ shown in figure 2. The 3D vorticity
equation is given by

$\frac{\partial\boldsymbol{\omega}}{\partial
t}+\nabla\cdot(\mathbf{u}\boldsymbol{\omega})=\nabla\cdot(\boldsymbol{\omega}\mathbf{u})+\nu\nabla^{2}\boldsymbol{\omega},$
(27)

and the total rate of change of $\boldsymbol{\Omega}$ is then expressed by

$\displaystyle\frac{D\boldsymbol{\Omega}}{Dt}$ $\displaystyle=$
$\displaystyle\int_{V}\frac{D\boldsymbol{\omega}}{Dt}dV=\frac{D}{Dt}\int_{V}\boldsymbol{\omega}dV=\int_{V}\nabla\cdot\left(\boldsymbol{\omega}\mathbf{u}+\nu\nabla\boldsymbol{\omega}\right)dV$
(28) $\displaystyle=$ $\displaystyle\int_{\partial
V}\hat{\mathbf{n}}\cdot(\boldsymbol{\omega}\mathbf{u})dS+\int_{\partial
V}\hat{\mathbf{n}}\cdot(\nu\nabla\boldsymbol{\omega})dS.$

The first term on the RHS of Eq. (4) corresponds to vortex stretching/tilting
and is present only in 3D. We now write

$\displaystyle\frac{D}{Dt}\int_{V_{1}UV_{2}}\boldsymbol{\omega}dV=$
$\displaystyle\oint_{\partial
V_{1}}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}\mathbf{u}+\nu\nabla\boldsymbol{\omega}\right)dS+\oint_{\partial
V_{2}}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}\mathbf{u}+\nu\nabla\boldsymbol{\omega}\right)dS$
(29) $\displaystyle+$ $\displaystyle\oint_{\partial
V_{1}^{\prime}}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}_{1}\mathbf{u}_{1}+\nu_{1}\nabla\boldsymbol{\omega}_{1}\right)dS+\oint_{\partial
V_{2}^{\prime}}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}_{2}\mathbf{u}_{2}+\nu_{2}\nabla\boldsymbol{\omega}_{2}\right)dS,$

and let $V_{1}UV_{2}\rightarrow V$,
$\hat{\mathbf{n}}\rightarrow\hat{\mathbf{s}}$ from fluid 1,
$\hat{\mathbf{n}}\rightarrow-\hat{\mathbf{s}}$ from fluid 2, and $(\partial
V_{1},\partial V_{2})\rightarrow I$, then it follows that

$\frac{D}{Dt}\int_{V}\boldsymbol{\omega}dV=\oint_{\partial
V}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}\mathbf{u}+\nu\nabla\boldsymbol{\omega}\right)dS-\left[\oint_{I}[[\hat{\mathbf{s}}\cdot(\boldsymbol{\omega}\mathbf{u})]]dS+\oint_{I}[[\nu\hat{\mathbf{s}}\cdot\nabla\boldsymbol{\omega}]]dS\right].$
(30)

It is important to establish a connection between
$\oint_{I}[[\nu\hat{\mathbf{s}}\cdot\nabla\boldsymbol{\omega}]]dS$, which
represents the jump across the plane of the interface of the vorticity flux,
and the momentum conservation equation given by

$\frac{D\mathbf{u}}{Dt}=-\frac{\nabla
p}{\rho}-\nu\nabla\times\boldsymbol{\omega}.$ (31)

In order to relate this term to the $\nu\nabla\times\boldsymbol{\omega}$ term
in Eq. (31), we first write down the following general result 222We have used
the vector identity $\int_{V}\nabla\times\mathbf{A}dV=-\oint_{\partial
V}\mathbf{A}\times d\mathbf{S}=-\oint_{\partial
V}\mathbf{A}\times\mathbf{n}dS=\oint_{\partial
V}\mathbf{n}\times\mathbf{A}dS,$ for any vector $\mathbf{A}$, and volume $V$
enclosed by a surface $\partial V$ with a unit normal $\mathbf{n}$.

$\displaystyle-\oint_{\partial
V}\hat{\mathbf{s}}\cdot\nabla\boldsymbol{\omega}dS$ $\displaystyle=$
$\displaystyle-\int_{V}\nabla^{2}\boldsymbol{\omega}dV=-\int_{V}\left(\nabla(\nabla\cdot\boldsymbol{\omega})-\nabla\times\nabla\times\boldsymbol{\omega}\right)dV=\int_{V}\nabla\times\nabla\times\boldsymbol{\omega}dV$
(32) $\displaystyle=$ $\displaystyle-\oint_{\partial
V}\left(\nabla\times\boldsymbol{\omega}\right)\times\hat{\mathbf{s}}dS=\oint_{\partial
V}\hat{\mathbf{s}}\times\nabla\times\boldsymbol{\omega}dS.$

Note that this relation links Lighthill’s vorticity flux to Lyman’s flux, the
latter being another form of the former (see Terrington et al. (2021) and
references therein).

Inspired by the form of Lyman’s flux, the natural way to proceed is to take
the cross product of $\hat{\mathbf{s}}=\hat{\mathbf{t}}\times\hat{\mathbf{b}}$
with the LHS of Eq. (31) and its pressure gradient term 333We have exploited
the fact that
$\hat{\mathbf{t}}\times\hat{\mathbf{b}}\times\mathbf{c}=\hat{\mathbf{b}}(\hat{\mathbf{t}}\cdot\mathbf{c})-\mathbf{c}(\hat{\mathbf{t}}\cdot\hat{\mathbf{b}})=\hat{\mathbf{b}}(\hat{\mathbf{t}}\cdot\mathbf{c})$
since $\hat{\mathbf{t}}\cdot\hat{\mathbf{b}}=0$. and a cross product of
$\hat{\mathbf{s}}$ with its $\nu\nabla\times\boldsymbol{\omega}$ term to
arrive at

$\displaystyle-\nu\hat{\mathbf{s}}\times\nabla\times\boldsymbol{\omega}$
$\displaystyle=\hat{\mathbf{b}}\hat{\mathbf{t}}\cdot\frac{D\mathbf{u}}{Dt}-\hat{\mathbf{b}}\hat{\mathbf{t}}\cdot\nabla\left(\frac{p}{\rho}\right)$
$\displaystyle=\hat{\mathbf{b}}\left[\left(\frac{D}{Dt}(\mathbf{u}\cdot\hat{\mathbf{t}})-\mathbf{u}\cdot\frac{D\hat{\mathbf{t}}}{Dt}\right)+\hat{\mathbf{t}}\cdot\nabla\left(\frac{p}{\rho}\right)\right]$
$\displaystyle=\nu\hat{\mathbf{s}}\cdot\nabla\boldsymbol{\omega};$

here, we note that the sources of vorticity are due to acceleration in the
plane of the interface, which we can think of as a vortex sheet, and
interfacial pressure gradients. Making use of this relation in Eq. (30), we
arrive at

$\displaystyle\frac{D}{Dt}\left[\int_{V}\boldsymbol{\omega}dV+\hat{\mathbf{b}}\oint_{I}[[\mathbf{u}\cdot\hat{\mathbf{t}}]]dS\right]$
$\displaystyle=$ $\displaystyle\oint_{\partial
V}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}\mathbf{u}+\nu\nabla\boldsymbol{\omega}\right)dS$
$\displaystyle-\oint_{I}[[\hat{\mathbf{s}}\cdot(\boldsymbol{\omega}\mathbf{u})]]dS+\oint_{I}\hat{\mathbf{b}}[[\mathbf{u}\cdot\frac{D\hat{\mathbf{t}}}{Dt}]]dS-\oint_{I}\hat{\mathbf{b}}\frac{\partial}{\partial
s}[[\frac{p}{\rho}]]dS,$

where we have set
$\hat{\mathbf{t}}\cdot\nabla(p/\rho)=\partial(p/\rho)/\partial s$. An
expression for $\mathbf{u}\cdot(D\hat{\mathbf{t}}/Dt)$ can be developed given
by (the details are in A)

$\mathbf{u}\cdot\frac{D\hat{\mathbf{t}}}{Dt}=\frac{1}{2}\frac{\partial}{\partial
s}\left[(\mathbf{u}\cdot\hat{\mathbf{s}})^{2}+(\mathbf{u}\cdot\hat{\mathbf{b}})^{2}\right]+\frac{1}{2}\frac{\partial}{\partial
b}\left[(\mathbf{u}\cdot\hat{\mathbf{s}})^{2}+(\mathbf{u}\cdot\hat{\mathbf{b}})^{2}\right]-\kappa_{1}(\mathbf{u}\cdot\hat{\mathbf{t}})(\mathbf{u}\cdot\hat{\mathbf{s}}).$
(34)

Furthermore, for $[[\rho]]=0$, the remaining term required to close equation
LABEL:eq:circulation_jump_3D_2 is one for $[[p]]$ (the details are in B):

$[[p]]=-\sigma(\kappa_{1}+\kappa_{2})-2[[\mu\left(\frac{\partial}{\partial
s}\left[(\mathbf{u}\cdot\hat{\mathbf{t}})+(\mathbf{u}\cdot\hat{\mathbf{b}})\right]+(\kappa_{1}+\kappa_{2})(\mathbf{u}\cdot\hat{\mathbf{s}})\right)]].$
(35)

To collapse these equations to their two-dimensional (2D) equivalents, we
first note that
$\hat{\mathbf{s}}\cdot\boldsymbol{\omega}=\hat{\mathbf{n}}\cdot\boldsymbol{\omega}=\mathbf{u}\cdot\hat{\mathbf{b}}=0$
in 2D, and set $\partial/\partial b=0$; the latter leads to $\kappa_{2}=0$. We
then take a dot product of Eq. (LABEL:eq:circulation_jump_3D_2) with
$\hat{\mathbf{b}}$ (and convert the volume and area integrals to area and line
integrals, respectively) to arrive at a 2D analogue involving the vorticity
scalar $\omega$. Moreover, in the case studied here, characterised by
$[[\mu]]=0$, $[[\mathbf{u}\cdot\hat{\mathbf{s}}]]=0$,
$[[\mathbf{u}\cdot\hat{\mathbf{t}}]]=0$, and
$[[\mathbf{u}\cdot\hat{\mathbf{b}}]]=0$, equation
(LABEL:eq:circulation_jump_3D_2) reduces to

$\frac{D}{Dt}\left[\int_{V}\boldsymbol{\omega}dV\right]=\oint_{\partial
V}\hat{\mathbf{n}}\cdot\left(\boldsymbol{\omega}\mathbf{u}+\nu\nabla\boldsymbol{\omega}\right)dS-\oint_{I}[[\hat{\mathbf{s}}\cdot(\boldsymbol{\omega}\mathbf{u})]]dS+\frac{1}{\rho}\oint_{I}\hat{\mathbf{b}}\frac{\partial}{\partial
s}\left(\sigma\left[\kappa_{1}+\kappa_{2}\right]\right)dS.$ (36)

We note that the term involving
$[[\hat{\mathbf{s}}\cdot(\boldsymbol{\omega}\mathbf{u})]]$ on the right-hand-
side of this equation is zero. To see this, we first note that
$[[\hat{\mathbf{s}}\cdot\boldsymbol{\omega}\mathbf{u}]]$ can be re-expressed
as

$\displaystyle[[\hat{\mathbf{s}}\cdot\boldsymbol{\omega}\mathbf{u}]]$
$\displaystyle=$
$\displaystyle(\hat{\mathbf{s}}\cdot\boldsymbol{\omega}_{2})\mathbf{u}_{2}-(\hat{\mathbf{s}}\cdot\boldsymbol{\omega}_{1})\mathbf{u}_{1}$
(37) $\displaystyle=$
$\displaystyle(\hat{\mathbf{s}}\cdot\boldsymbol{\omega}_{2}-\hat{\mathbf{s}}\cdot\boldsymbol{\omega}_{1})\mathbf{u}_{1}$
$\displaystyle=$
$\displaystyle(\hat{\mathbf{s}}\cdot\boldsymbol{\omega}_{2}-\hat{\mathbf{s}}\cdot\boldsymbol{\omega}_{1})\mathbf{u}_{2}$
$\displaystyle=$
$\displaystyle[[\hat{\mathbf{s}}\cdot\boldsymbol{\omega}]]\mathbf{u}_{1}=[[\hat{\mathbf{s}}\cdot\boldsymbol{\omega}]]\mathbf{u}_{2},$

since $[[\mathbf{u}]]=0$. We also note that
$\hat{\mathbf{s}}\cdot\boldsymbol{\omega}=(\hat{\mathbf{b}}\times\hat{\mathbf{t}})\cdot(\nabla\times\mathbf{u})$,
which can be re-written as

$\displaystyle\hat{\mathbf{s}}\cdot\boldsymbol{\omega}$ $\displaystyle=$
$\displaystyle\hat{\mathbf{b}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{t}}-\hat{\mathbf{t}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{b}}$
(38) $\displaystyle=$
$\displaystyle\hat{\mathbf{b}}\cdot\frac{\partial\mathbf{u}}{\partial
s}-\hat{\mathbf{t}}\cdot\frac{\partial\mathbf{u}}{\partial b}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
s}(\hat{\mathbf{b}}\cdot\mathbf{u})-\frac{\partial}{\partial
b}(\hat{\mathbf{t}}\cdot\mathbf{u}),$

since $\hat{\mathbf{b}}\neq\hat{\mathbf{b}}(s)$ and
$\hat{\mathbf{t}}\neq\hat{\mathbf{t}}(b)$. Thus, we can write

$\displaystyle[[\hat{\mathbf{s}}\cdot\boldsymbol{\omega}]]$ $\displaystyle=$
$\displaystyle[[\frac{\partial}{\partial
s}(\hat{\mathbf{b}}\cdot\mathbf{u})]]-[[\frac{\partial}{\partial
b}(\hat{\mathbf{t}}\cdot\mathbf{u})]]$ (39) $\displaystyle=$
$\displaystyle\frac{\partial}{\partial
s}[[\hat{\mathbf{b}}\cdot\mathbf{u}]]-\frac{\partial}{\partial
b}[[\hat{\mathbf{t}}\cdot\mathbf{u}]]=0,$

since $[[\hat{\mathbf{b}}\cdot\mathbf{u}]]=0$ and
$[[\hat{\mathbf{t}}\cdot\mathbf{u}]]=0$, whence
$[[\hat{\mathbf{s}}\cdot\boldsymbol{\omega}\mathbf{u}]]=0$. Inspection of the
terms remaining in equation 36 suggests that circulation is influenced by
vorticity diffusion, vortex tilting/stretching, and gradients of curvature and
interfacial tension.

The dimensionless versions of equations (25) and (24) are then expressed by

$\displaystyle\left[\left[\tilde{\omega}_{t}\right]\right]$ $\displaystyle=$
$\displaystyle-2ReMa\frac{1}{(1-\tilde{\Gamma})}\frac{\partial\tilde{\Gamma}}{\partial
b},$ (40) $\displaystyle\left[\left[\tilde{\omega}_{b}\right]\right]$
$\displaystyle=$
$\displaystyle-2ReMa\frac{1}{(1-\tilde{\Gamma})}\frac{\partial\tilde{\Gamma}}{\partial
s},$ (41)

and the dimensionless equation (36) reads

$\frac{D}{D\tilde{t}}\left[\int_{\tilde{V}}\tilde{\boldsymbol{\omega}}d\tilde{V}\right]=\oint_{\partial\tilde{V}}\hat{\mathbf{n}}\cdot\left(\tilde{\boldsymbol{\omega}}\tilde{\mathbf{u}}+\frac{1}{Re}\nabla\tilde{\boldsymbol{\omega}}\right)d\tilde{S}+\frac{1}{We}\oint_{I}\hat{\mathbf{b}}\frac{\partial}{\partial\tilde{s}}\left(\tilde{\sigma}\left[\tilde{\kappa}_{1}+\tilde{\kappa}_{2}\right]\right)d\tilde{S},$
(42)

and the tildes are dropped henceforth.

We note that in the limit of small variations of surfactant concentration
around its initial value (i.e., so-called diluted systems),
$\hat{\Gamma}=\Gamma_{0}+\delta\hat{\Gamma}$, with
$\delta\hat{\Gamma}\ll\Gamma_{0}$, leads to
$\hat{\Gamma}=1+\delta\hat{\Gamma}$, and the equation of state can be
linearized to result in $\hat{\sigma}=1-\beta_{s}\hat{\Gamma}$. Note that, in
the case of non-isothermal systems, $\hat{\sigma}$ has a linear dependency on
the local temperature (T), and a linear equation of state describes
$\hat{\sigma}(\hat{T})$ (see for example Williams et al. (2021)). Therefore,
surface tension gradients in equation (42) can also arise due to thermal
gradients.

##  3 Results

Figure 3: Regime map of the interfacial morphology in the $\beta_{s}-We$ space
for $Re=5000$, $Pe_{s}=100$, and $\Gamma_{o}=\Gamma_{\infty}/2$. The
capillary-dominated and inertia-dominated regimes, and their boundaries are
clearly demarcated; the symbols represent simulations carried out at the
transition lines separating these regimes. Three-dimensional representations
of the interface for both regimes are also shown.  Figure 4:  Spatio-temporal
representation of the interfacial dynamics and the coherent vortical
structures for $Re=5000$ and $We=500$ at
$t=(23.28,~{}28.12,~{}31.25,~{}31.50)$ shown in (a)-(d), respectively. The
three-dimensional coherent structures are visualised by the $Q$-criterion with
values of $Q=(3,~{}220,~{}320,~{}320)$, where the colour represents the
streamwise vorticity field, $\omega_{t}$. In each panel, we also show
$\omega_{x}$ in the $y$–$z$ plane for each sampling location.

Figure 3 shows a flow regime map for $Re=5000$ that depicts the interfacial
morphology associated with various regions of the $\beta_{s}-We$ parameter
space generated by over 100 transient simulations performed in the ranges
$100<We<1000$ and $0.1<\beta_{s}<0.9$. We have divided the map into two
distinct regions depending on the morphology: for small $We$, capillary forces
control the interfacial dynamics preventing the development of lobes which
could result in the formation of large droplets; for large $We$, inertial
forces dominate the dynamics triggering the formation of interfacial lobes
whose thinning eventually results in the generation of holes and eventually
droplets. Inspection of Figure 3 reveals that the presence of Marangoni
stresses counteracts the transition from the low- to high-$We$ regimes as the
critical $We$ increases with $\beta_{s}$.

To assess the effect of Marangoni-induced flow, we have analysed the flow
physics of the surfactant-free and surfactant-laden flows characterised by
$Re=5000$ and $We=500$. We start with the surfactant-free case depicted in
Figure 4 which shows the spatio-temporal interfacial dynamics for the
surfactant-free case through the $Q$-criterion (e.g., a measure of the
dominance of vorticity $\boldsymbol{\omega}$ over strain $\bf{s}$, i.e.,
$Q=(||\boldsymbol{\omega}||^{2}-||\bf{s}||^{2})/2$ (Hunt et al., 1988)). At
early times, we observe the formation of a periodic array of quasi-symmetric
Kelvin-Helhomltz (KH)-driven vortex rings as a result of the difference in
velocity in the shear layer located under the interface (see figure 4a). With
increasing time, the three-dimensional instability starts with the deformation
of the vortex-rings leading to a mutual-induction between two consecutive
vortex rings resulting in their ‘knitting’ (see figure 4b); similar vortex-
pairing has been reported by Broze & Hussain (1996) and da Silva & Métais
(2002). With increasing time, we observe the formation of inner and outer
hairpin vortices whose pairing brings about a region where both overlap. The
cascade mechanism resulting in the formation of hairpin-vortices from KH-rings
is triggered by the magnitude of the streamwise vorticity, $\omega_{x}$, which
becomes comparable to its azimuthal counterpart, $\omega_{y}$, in agreement
with Jarrahbashi et al. (2016) and Constante-Amores et al. (2021a), as shown
in figure 4b.

To provide more conclusive evidence of the existence of inner/outer hairpin
vortices in the jet dynamics, a careful study of the distribution of vortex
signs shows the assembling into counter-rotating vortex pairs (see
$\omega_{x}$ in the $y$-$z$ plane for each sampled location of the panels in
figure 4). By analysing the distribution of streamwise vorticity between the
ring and braid regions of the jet core (see figure 4a), we observe that their
distribution is $\pi$-out-of-phase. The arrangement of the vorticity comes
from vortex induction arguments, similar to those explained by Jarrahbashi et
al. (2016), Zandian et al. (2018) and Constante-Amores et al. (2021a), i.e.,
the upstream hairpin vortex from the ring overtakes the upstream hairpin
vortex from the braid as the mutual induction takes place. Finally, the
vortex-surface interaction triggers the formation of the interfacial structure
as the interface adopts the shape of the vortex which is in its vicinity (see
figure 4b-d). The mutual induction between outer and inner hairpin vortices
eventually leads to the thinning of the lobes to ultimately form inertia-
induced holes whose capillary-driven expansion gives rise to the formation of
droplets (Jarrahbashi et al., 2016; Zandian et al., 2018; Constante-Amores et
al., 2021a).

---
(a) (b)
(c)
(d)
(e) (f)
Figure 5:  Effect of surfactants on the early interfacial dynamics for
$Re=5000$, $We=500$, $\beta_{s}=0.5$, $Pe_{s}=100$ and
$\Gamma_{o}=\Gamma_{\infty}/2$ at $t=32.03$, (a), and 32.81, (b). The top and
bottom panels represent the interface coloured by $\Gamma$ and the coherent
vortical structures visualised via $Q$-criterion with $Q=10$. Panel (c) shows
a 2D representation of $\Gamma$, and $\tau$, with respect to the arc-length
$s$ (see inset) at $t=32.03$. Panel (d) shows a 2D representation of the
magnitude of vorticity $|\boldsymbol{\omega}|$ in the $x$–$z$ plane
($z=2.875$) at $t=32.81$; interface location and $\Gamma$, and
$[[\omega_{b}]]$ and $[[\omega_{t}]]$ vs the arc-length $s$, shown in (e) and
(f), respectively.

---
(a)

|  |  |
---|---|---|---
(b) | (c) | (d) | (e)

|
---|---
(f) | (g)

Figure 6: Surfactant-driven hairpin-vortical structures: temporal development
of the HV1 and HV2 hairpin-like vortical structures via $Q$-criterion with
$Q=200$ together with the interfacial location coloured by $\Gamma$ at time
$t=(32.81,~{}33.59,~{}34.37,~{}36.71,~{}39.06)$, corresponding to panels
(a)-(e), respectively. In panel (e), we show the two transversal slices
displayed in panels (f) and (g) which depict the streamwise vorticity
$\omega_{t}$ through the legs and head of $HV1$, respectively; arrows of in-
plane velocity vectors have been added. The parameter values are the same as
in figure 5.

Next, we turn our attention to the effect of surfactants on the flow dynamics.
Figure 5 shows the early interfacial surfactant concentration together with
the three-dimensional coherent vortical structures via the $Q$-criterion.
Similarly to the surfactant-free case, we observe the formation of a periodic
array of quasi-axisymmetric KH-vortex rings. These rings induce the formation
of interfacial waves that are characterised by regions of radially converging
and diverging motion that lead to higher and lower interfacial areas, and
subsequently to lower and higher surfactant concentration regions,
respectively; accumulation of $\Gamma$ is observed in the vicinity of the KH
rings (see figure 5a). Figure 5c presents the interfacial concentration
$\Gamma$, and Marangoni stresses $\tau$ along an arc-length, $s$,
corresponding to $t=32.03$. We observe that the non-uniform distribution of
$\Gamma$ gives rise to Marangoni-induced flow, which drives fluid motion from
ring-1, ‘VR1’, ($\tau>0$) to ring-2, ‘VR2’, and vice versa (i.e, flow from VR2
to VR1, $\tau<0$). This flow is therefore accompanied by the retardation of
the development of the interfacial waves and a subsequent delay of the onset
of the three-dimensional instability of the jet observed in the surfactant-
free case in figure 4.

Additionally, these Marangoni stresses promote jumps in the vorticity across
the interface which we can calculate using equations 21 and 22 in the location
which coincides with the formation of vortex $SV1$ and $SV2$ from figure 6 at
$t=32.81$. Figure 5d shows a three-dimensional representation of the interface
together with an $x$-$z$ plane at $y=2.875$ colored by the the magnitude of
vorticity, $|\boldsymbol{\omega}|$. Figure 5e,f show respectively the
variation of the interface location and the $\Gamma$ profiles, and of the
distribution of $[[\omega_{b}]]$ and $[[\omega_{t}]]$, along the arc length,
$s$, in the plane cutting the interface shown in figure 5d. From figure 5e, it
is seen that the surfactant accumulates in the down-sloping region immediately
downstream of an interfacial wave peak; here, the gradients in $\Gamma$, and
therefore in $\sigma$, are smallest corresponding to the weakest vorticity
jumps, while the largest such jumps are in the wave peak and trough regions
where the $\Gamma$ (and $\sigma$) gradients are highest, as shown in figure
5f. Inspection of figure 5f also shows that $[[\omega_{b}]]\gg[[\omega_{t}]]$,
that is, near-interface vorticity production in the azimuthal direction is
dominant. This acts to disrupt the dynamics of vortex-pairing relative to the
surfactant-free case as the ‘knitting process’ is promoted by streamwsie
rather than azimuthal vorticity production and the vortex-ring deformation is
replaced by vortex-reconnection and merging in the azimuthal direction in the
surfactant-laden case.

|
---|---
(a) | (b)
|
(c) | (d)
|
(e) | (f)
Figure 7:  Panels (a-b) and (c-d) show the location of the interface together
with the surfactant concentration and the jumps of the vorticity across the
interface for $t=36.81$ and $t=44.68$, respectively. Panels (e-f) show a three
dimensional representation of the interface location coloured by $\Gamma$ at
$t=44.68$, and vortex knitting visualised via $Q$-criterion with value of
$Q=1600$ where the colour represents $\omega_{t}$, respectively. The center of
the jet core corresponds to $z=2.5$. The parameter values are the same as in
figure 5.

For increasing time, figure 6 shows the formation of surfactant-induced inner
hairpin-like vortical structures. The shear stress, which is generated to
balance the gradients in $\sigma$ gives rise to counter-rotating streamwise
vortices of similar magnitude to the KH rings (labelled ‘SV1’ and ‘SV2’ in
figure 6b). These structures grow in the $x-$direction into a combination of
streamwise vortices close to the interface, i.e. legs, and a hairpin-like head
close to the center-plane of the jet (see figure 6d). The hairpin-legs extend
from the regions of high-to-low values of $\Gamma$ on the surface, while the
hairpin-head points down in the positive $x-$direction (labelled ‘HV1’ and
‘HV2’ in figure 6e). To complete the presentation of these hairpin-like
vortical structures, figure 6f,g show the direction of flow rotation of the
legs and head for HV1. For comparison, we have added arrows to show velocity
direction and to prove that this coherent vortical structure exhibits the same
qualitative behaviour as the HV proposed by Theodorsen (1952) for near-wall
turbulence. To the best of our knowledge, the formation of hairpin-like
vortical structures induced by surfactant effects has not been reported yet.
We have also observed surfactant-driven outer hairpin-like vortical structures
(not shown) whose heads are in the negative $x-$direction (in the frame of
reference of the legs).

Surfactant-laden case | Surfactant-free case
---|---
|
(a) | (b)
|
(c) | (d)
|
(e) | (f)
|
(g) | (h)
Figure 8:  Total rate of change of circulation, $\boldsymbol{\Omega}$,
according to equation 42: $D\boldsymbol{\Omega}/Dt$, vortex diffusion ($I_{\rm
diff}$), vortex tilting ($I_{\rm tilt}$), and surface tension ($I_{\rm curv}$)
are shown in rows one to four, respectively; see equation (44) for the
definitions of $I_{\rm tilt}$, $I_{\rm diff}$, and $I_{\rm curv}$. Surfactant-
laden and surfactant-free cases correspond to left and right panels,
respectively. For panel (g), we represent the contributions that arise from
the gradients of curvature (solid lines) and the gradients of surface tension
(dashed lines) to underscore the relative importance of the Marangoni
stresses. Red, blue, and black colored lines represent component $x$, $y$ and
$z$ of $D\boldsymbol{\Omega}/Dt$, $I_{\rm tilt}$, diff, and $I_{\rm curv}$.
The parameters are $Re=5000$, $We=500$ (and for the surfactant-laden case)
$\beta_{s}=0.5$, $Pe_{s}=100$ and $\Gamma_{o}=\Gamma_{\infty}/2$.

At later times, figure 7a-d shows the variation with arc length of the
interfacial location, $\Gamma$, and $[[\omega_{t}]]$ and $[[\omega_{b}]]$ at
$t=36.51$ and $t=44.68$; corresponding three-dimensional representations of
the interface are also shown in figure 7e,f for $t=44.68$ coloured by the
magnitude of $\Gamma$ and the $Q$-criterion, respectively. The flow is
accompanied by radially-converging and diverging motion due to vortex-surface-
interaction; interfacial convection drives surfactant towards the inner lobes
(interfacial contraction), and away from the outer lobes (interfacial
expansion). Vorticity jumps are highest in the interfacial regions with the
largest gradients in $\Gamma$. As time evolves, the ratio of these Marangoni-
driven $[[\omega_{t}]]$ to $[[\omega_{b}]]$ reduces and this results in large
coherent structures which merge to form counter-rotating streamwise vortical
rings that eventually ‘knit’ with the adjacent vortex ring located in the
$x-$direction (labelled ‘VR1-VR4’ in figure 7f); this pairing is similar to
the surfactant-free case (in agreement with Urbin & Métais (1997) and da Silva
& Métais (2002)).

|  |  |
---|---|---|---
(a) | (b) | (c) | (d)
|  |  |
(e) | (f) | (g) | (h)
Figure 9:  Spatio-temporal evolution of the control volume $V$ with a surface
$\partial V$ which encloses an interface plane $I$ used to calculate the rate
of change of the circulation in figure 8. Panels (a)-(d) and (e)-(h)
correspond to the surfactant-laden and surfactant-free cases, respectively, at
$t=(25.20,35.50,40.37,43.75)$ and $t=(27.12,30.43,34.37,40.55)$, and the same
parameters as in figure (8).

|
---|---
(a) | (b)

---
(c)

Figure 10: Temporal evolution of (a) surface area, $A$, (b) kinetic energy,
$E_{k}=\rho\int_{V}\mathbf{u}^{2}/2dV$, and (c) enstrophy,
$\varepsilon=\int_{V}|\boldsymbol{\omega}^{2}|dV$, scaled by the initial
interfacial area, $A_{0}$, kinetic energy, $E_{k0}$, and enstrophy,
$\varepsilon_{0}$, respectively. The parameter values are the same as in
figure 5.

We now examine the dynamics of the circulation $\boldsymbol{\Omega}$ by
considering equation 42 which we express as follows:

$\frac{D\boldsymbol{\Omega}}{D\tilde{t}}=I_{\rm tilt}+I_{\rm diff}+I_{\rm
curv},$ (43)

where $I_{\rm tilt}$, $I_{\rm diff}$, and $I_{\rm curv}$ are defined as

$I_{\rm
tilt}\equiv\oint_{\partial\tilde{V}}\hat{\mathbf{n}}\cdot\tilde{\boldsymbol{\omega}}\tilde{\mathbf{u}}d\tilde{S},~{}~{}I_{\rm
diff}\equiv\frac{1}{Re}\oint_{\partial\tilde{V}}\hat{\mathbf{n}}\cdot\nabla\tilde{\boldsymbol{\omega}}d\tilde{S},~{}~{}I_{\rm
curv}\equiv\frac{1}{We}\oint_{I}\hat{\mathbf{b}}\frac{\partial}{\partial\tilde{s}}\left(\tilde{\sigma}\left[\tilde{\kappa}_{1}+\tilde{\kappa}_{2}\right]\right)d\tilde{S},$
(44)

which correspond to vortex tilting/stretching, diffusion of vorticity, and
circulation variation due to gradients in curvature and interfacial tension
(in the case of surfactant-laden systems). Figure 8 shows the temporal
evolution of $D\boldsymbol{\Omega}/Dt$, $I_{\rm tilt}$, $I_{\rm diff}$, and
$I_{\rm curv}$ which allows us to identify the dominant physical mechanisms
that contribute to the creation and dissipation of circulation. In figure 9 we
also show snapshots of the three-dimensional representation of the interface
corresponding to the volume used to carry out the computations necessary to
calculate $D\boldsymbol{\Omega}/Dt$ and its constituent terms for the
surfactant-laden and surfactant-free cases; this allows one to pinpoint the
mechanisms primarily responsible for the interfacial structures observed. It
is clearly seen from figure 8 that during the early stages of the flow,
$\boldsymbol{\Omega}$ remains approximately constant. Inspection of panels
(c)-(h) of figure 8 shows clearly that the rate of change of circulation is
dominated by the mechanisms related to vortex diffusion $I_{\rm diff}$ and
curvature $I_{\rm curv}$, with vortex tilting/shielding playing a relatively
minor role. It is also clear that in the surfactant-laden jet case, the
Marangoni contribution to $I_{\rm curv}$ dominates that associated with
curvature derivatives. This observation further bolsters the claim that
Marangoni stresses drive vorticity generation in the jet dynamics.

The snapshots depicted in figure 9 for the surfactant-laden (panels (a)-(d))
and surfactant-free (panels (e)-(h)) cases have been chosen carefully so as to
link the various stages of jet destabilisation to the prominent changes in the
temporal variation of $I_{\rm diff}$, $I_{\rm tilt}$, $I_{\rm curv}$, and
$D\boldsymbol{\Omega}/Dt$. Given the dominance of $I_{\rm curv}$ over the time
range considered ($0\leq t\leq 40$), we focus on the variations in this
quantity and its signature effects on the interfacial shape. Inspection of
figures 8(g) and 9(a) reveals that the relatively gentle interfacial
undulations are linked to variations of the Marangoni contribution to $I_{\rm
curv}$ in the $x-y$ plane. The development of the more complex interfacial
shapes, on the other hand, is accompanied by a concomitant rise in three-
dimensionality of $I_{\rm curv}$ (in addition to significant contributions
from the $x-$component of $I_{\rm curv}$). In the surfactant-free case,
inspection of figures (8)(d) and (h), and (9)(e)-(h) shows that the
interfacial jet evolution is accompanied by large variations in the
$x-$component of $I_{\rm curv}$ and vorticity diffusion characterised by
$I_{\rm diff}$.

Lastly, we plot in figure 10 the effect of surfactants on the interfacial
area, kinetic energy, defined as $E_{k}=\rho\int_{V}{\bf u}^{2}/2dV$, and the
enstrophy, $\varepsilon=\int_{V}|\boldsymbol{\omega}^{2}|dV$, normalised by
their initial values, $A_{0}$, $E_{k0}$, $\varepsilon_{0}$, respectively.
After the onset of destabilization (defined when the interfacial surface has
reached $A=1.025$), we observe that the surfactant-induced effects discussed
above, which include the interfacial vorticity jumps brought about by
Marangoni stresses, and their effect on the production of circulation, and jet
destabilisation mechanisms associated with vortex formation and spanwise
reconnection, promote the delay in increase and subsequent reduction in
interfacial area; these effects also lead to a delay in the decay of the jet
kinetic energy as well as its enstrophy.

##  4 Concluding remarks

Three-dimensional numerical simulations of jet destabilisation and atomisation
in the presence of a monolayer of insoluble surfactants have been carried out
for the first time. A phase diagram in the space of dimensionless surfactant
elasticity and Weber number in the inertia-dominated region is presented in
the limiting case where there is no vorticity production associated with jumps
in material properties such as fluid density and viscosity; in the present
work, surface tension forces and Marangoni stress give rise to variations in
vorticity and circulation in addition to the vortex tilting/shielding and
diffusion mechanisms. We have derived formulae for the vorticity jumps across
the interface due to Marangoni stresses, and equations that provide a
breakdown of the rate of production of circulation within the jet into
constituent terms which we associate with vortex tilting/shielding, diffusion,
and gradients in interfacial curvature and surface tension.

Then, we have analysed in details the vortex-interface-surfactant interactions
in the flow dynamics. At early times, the presence of surfactants induces
spanwise vortex reconnections brought about Marangoni-induced flow resulting
in the delay of the onset of destabilisation to the three-dimensional
interfacial instabilities. We also show that surfactant-induced Marangoni-
stresses trigger the formation of hairpin-like structures whose head and legs
extend in the streamwise direction. Lastly, we have attempted to link the
changes in interfacial topology to the mechanisms that influence the
production of vorticity and circulation demonstrating a balance between
curvature gradients and diffusion for surfactant-free jets, and the dominance
of Marangoni stresses in the surfactant-laden cases.

The present results have been obtained for insoluble surfactants, and we
acknowledge that experimental studies feature soluble surfactants which are
dissolved in the liquid that issues from a nozzle to form the jet. Therefore,
future studies will be focused on examining the role of solubility on the flow
dynamics.

Declaration of Interests. The authors report no conflict of interest.

This work is supported by the Engineering and Physical Sciences Research
Council, United Kingdom, through the EPSRC MEMPHIS (EP/K003976/1) and PREMIERE
(EP/T000414/1) Programme Grants. O.K.M. acknowledges funding from PETRONAS and
the Royal Academy of Engineering for a Research Chair in Multiphase Fluid
Dynamics. We acknowledge HPC facilities provided by the Research Computing
Service (RCS) of Imperial College London for the computing time. D.J. and J.C.
acknowledge support through HPC/AI computing time at the Institut du
Developpement et des Ressources en Informatique Scientifique (IDRIS) of the
Centre National de la Recherche Scientifique (CNRS), coordinated by GENCI
(Grand Equipement National de Calcul Intensif) Grant 2022 A0122B06721.

##  Appendix A Kinematics

We first develop an expression for $D\hat{\mathbf{t}}/Dt$. We consider the
motion of an infinitesimal fluid parcel in the plane of the interface, which
is treated as a material surface. The position vector is
$\mathbf{x}=\mathbf{x}(s,b,t)$ where $s$ and $b$ represent arclength distances
along the $\hat{\mathbf{t}}$ and $\hat{\mathbf{b}}$ directions, respectively.
At time, $t+\delta t$, to leading order in $\delta t$, we can write the
following expression for the tangent to the interface at the fluid parcel
which at time $t$ was located at $\mathbf{x}(0,0,t)$

$\tilde{\mathbf{t}}(t+\delta t)=\frac{\partial\mathbf{x}}{\partial
s}(0,0,t)+\nabla\mathbf{u}\cdot\left(\frac{\partial\mathbf{x}}{\partial
s}+\frac{\partial\mathbf{x}}{\partial b}\right)\delta t.$ (45)

Noting that $\hat{\mathbf{t}}=\partial\mathbf{x}/\partial s$,
$\mathbf{b}=\partial\mathbf{x}/\partial b$,
$\hat{\mathbf{t}}\cdot\nabla\mathbf{u}=\partial\mathbf{u}/\partial s$, and
$\hat{\mathbf{b}}\cdot\nabla\mathbf{u}=\partial\mathbf{u}/\partial b$, this
equation can be re-expressed as

$\tilde{\mathbf{t}}(t+\delta
t)=\hat{\mathbf{t}}+\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta t.$ (46)

The magnitude of $\tilde{\mathbf{t}}(t+\delta t)$ is given by

$|\tilde{\mathbf{t}}|=1+\hat{\mathbf{t}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta t+O(\delta t)^{2},$ (47)

and normalisation of $\tilde{\mathbf{t}}(t+\delta t)$ by this magnitude gives
the following tangent unit vector $\hat{\mathbf{t}}(t+\delta t)$

$\displaystyle\hat{\mathbf{t}}(t+\delta t)$ $\displaystyle=$
$\displaystyle\frac{\hat{\mathbf{t}}+\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta
t}{|\tilde{\hat{\mathbf{t}}}(t+\delta
t)|}=\frac{\hat{\mathbf{t}}+\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta
t}{1+\hat{\mathbf{t}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta t}$ (48)
$\displaystyle=$
$\displaystyle\left(\hat{\mathbf{t}}+\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta
t\right)\left(1-\hat{\mathbf{t}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\delta t+O(\delta
t)^{2}\right)$ $\displaystyle=$
$\displaystyle\hat{\mathbf{t}}+\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial
b}-\hat{\mathbf{t}}\left(\hat{\mathbf{t}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\right)\right)\delta t+O(\delta
t)^{2}.$

From this expression, we can arrive at an approximate formula for
$D\hat{\mathbf{t}}/Dt$:

$\frac{D\hat{\mathbf{t}}}{Dt}\sim\frac{\hat{\mathbf{t}}(t+\delta
t)-\hat{\mathbf{t}}(t)}{\delta t}=\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial
b}-\hat{\mathbf{t}}\left(\hat{\mathbf{t}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right)\right)+O(\delta t).$ (49)

We now insert the following expression for $\mathbf{u}$ into $\mathbf{u}\cdot
D\hat{\mathbf{t}}/Dt$

$\mathbf{u}=(\mathbf{u}\cdot\hat{\mathbf{s}})\hat{\mathbf{s}}+(\mathbf{u}\cdot\hat{\mathbf{t}})\hat{\mathbf{t}}+(\mathbf{u}\cdot\hat{\mathbf{b}})\hat{\mathbf{b}},$
(50)

which yields

$\mathbf{u}\cdot\frac{D\hat{\mathbf{t}}}{Dt}=(\mathbf{u}\cdot\hat{\mathbf{s}})\hat{\mathbf{s}}\cdot\frac{D\hat{\mathbf{t}}}{Dt}+(\mathbf{u}\cdot\hat{\mathbf{b}})\hat{\mathbf{b}}\cdot\frac{D\hat{\mathbf{t}}}{Dt}.$
(51)

Substitution of Eq. (49) into $\hat{\mathbf{s}}\cdot D\hat{\mathbf{t}}/Dt$ and
$\hat{\mathbf{b}}\cdot D\hat{\mathbf{t}}/Dt$ gives

$\displaystyle\hat{\mathbf{s}}\cdot\frac{D\hat{\mathbf{t}}}{Dt}$
$\displaystyle=$
$\displaystyle\hat{\mathbf{s}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right),$ (52)
$\displaystyle\hat{\mathbf{b}}\cdot\frac{D\hat{\mathbf{t}}}{Dt}$
$\displaystyle=$
$\displaystyle\hat{\mathbf{b}}\cdot\left(\frac{\partial\mathbf{u}}{\partial
s}+\frac{\partial\mathbf{u}}{\partial b}\right),$ (53)

where we have made use of $\hat{\mathbf{s}}\cdot\hat{\mathbf{t}}=0$ and
$\hat{\mathbf{b}}\cdot\hat{\mathbf{t}}=0$. We can re-express the RHS of Eqs.
(52) and (53) as follows

$\displaystyle\hat{\mathbf{s}}\cdot\frac{\partial\mathbf{u}}{\partial s}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
s}(\hat{\mathbf{s}}\cdot\mathbf{u})-\mathbf{u}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial
s},$ (54)
$\displaystyle\hat{\mathbf{s}}\cdot\frac{\partial\mathbf{u}}{\partial b}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
b}(\hat{\mathbf{s}}\cdot\mathbf{u})-\mathbf{u}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial
b},$ (55)
$\displaystyle\hat{\mathbf{b}}\cdot\frac{\partial\mathbf{u}}{\partial s}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
s}(\hat{\mathbf{b}}\cdot\mathbf{u})-\mathbf{u}\cdot\frac{\partial\hat{\mathbf{b}}}{\partial
s},$ (56)
$\displaystyle\hat{\mathbf{b}}\cdot\frac{\partial\mathbf{u}}{\partial b}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
b}(\hat{\mathbf{b}}\cdot\mathbf{u})-\mathbf{u}\cdot\frac{\partial\hat{\mathbf{b}}}{\partial
b}.$ (57)

Inserting Eq. (50) into the second term on the RHS of Eqs. (54)-(57), we
obtain

$\displaystyle\mathbf{u}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial s}$
$\displaystyle=$ $\displaystyle\kappa_{1}(\mathbf{u}\cdot\hat{\mathbf{t}}),$
(58) $\displaystyle\mathbf{u}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial b}$
$\displaystyle=$ $\displaystyle\kappa_{2}(\mathbf{u}\cdot\hat{\mathbf{b}}),$
(59) $\displaystyle\mathbf{u}\cdot\frac{\partial\hat{\mathbf{b}}}{\partial s}$
$\displaystyle=$ $\displaystyle 0,$ (60)
$\displaystyle\mathbf{u}\cdot\frac{\partial\hat{\mathbf{b}}}{\partial b}$
$\displaystyle=$ $\displaystyle-\kappa_{2}(\mathbf{u}\cdot\hat{\mathbf{s}}),$
(61)

where the curvatures $\kappa_{1}$ and $\kappa_{2}$ are defined as follows

$\kappa_{1}=\hat{\mathbf{t}}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial
s},~{}~{}~{}~{}~{}\kappa_{2}=\hat{\mathbf{b}}\cdot\frac{\partial\hat{\mathbf{s}}}{\partial
b}.$ (62)

In deriving Eqs. (58)-(61), we have noted that
$\hat{\mathbf{t}}\neq\hat{\mathbf{t}}(b)$ and
$\hat{\mathbf{b}}\neq\hat{\mathbf{b}}(s)$. Substitution of Eqs. (58)-(61) into
Eqs. (54)-(57) and the resultant relations into Eqs. (52) and (53)
respectively yields the following expressions for
$\hat{\mathbf{s}}\cdot(D\hat{\mathbf{t}}/Dt)$ and
$\hat{\mathbf{b}}\cdot(D\hat{\mathbf{t}}/Dt)$

$\displaystyle\hat{\mathbf{s}}\cdot\frac{D\hat{\mathbf{t}}}{Dt}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
s}(\mathbf{u}\cdot\hat{\mathbf{s}})+\frac{\partial}{\partial
b}(\mathbf{u}\cdot\hat{\mathbf{s}})-\kappa_{1}(\mathbf{u}\cdot\hat{\mathbf{t}})-\kappa_{2}(\mathbf{u}\cdot\hat{\mathbf{b}}),$
(63) $\displaystyle\hat{\mathbf{b}}\cdot\frac{D\hat{\mathbf{t}}}{Dt}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
s}(\mathbf{u}\cdot\hat{\mathbf{b}})+\frac{\partial}{\partial
b}(\mathbf{u}\cdot\hat{\mathbf{b}})+\kappa_{2}(\mathbf{u}\cdot\hat{\mathbf{s}}).$
(64)

Substitution of Eqs. (63) and (64) into Eq. (51) and re-arranging yields

$\mathbf{u}\cdot\frac{D\hat{\mathbf{t}}}{Dt}=\frac{1}{2}\frac{\partial}{\partial
s}\left[(\mathbf{u}\cdot\hat{\mathbf{s}})^{2}+(\mathbf{u}\cdot\hat{\mathbf{b}})^{2}\right]+\frac{1}{2}\frac{\partial}{\partial
b}\left[(\mathbf{u}\cdot\hat{\mathbf{s}})^{2}+(\mathbf{u}\cdot\hat{\mathbf{b}})^{2}\right]-\kappa_{1}(\mathbf{u}\cdot\hat{\mathbf{t}})(\mathbf{u}\cdot\hat{\mathbf{s}}).$
(65)

##  Appendix B Near-interface normal stress jump

In order to generate a 3D version of the pressure gradient term in Eq.
(LABEL:eq:circulation_jump_3D_2), we first consider the jump in the normal
stress across the plane of the interface:

$p_{2}-p_{1}=-\sigma(\kappa_{1}+\kappa_{2})+[[\mu\hat{\mathbf{s}}\cdot\mathbf{D}\cdot\hat{\mathbf{s}}],$
(66)

where $\kappa_{1}$ and $\kappa_{2}$ are given by Eqs. (62). Substitution of
Eq. (50) into $\nabla\cdot\mathbf{u}=0$ yields

$\displaystyle\hat{\mathbf{s}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}$
$\displaystyle=$
$\displaystyle-\hat{\mathbf{t}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{t}}-\hat{\mathbf{b}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{b}}$
(67) $\displaystyle=$
$\displaystyle-\hat{\mathbf{t}}\cdot\frac{\partial\mathbf{u}}{\partial
s}-\hat{\mathbf{b}}\cdot\frac{\partial\mathbf{u}}{\partial b},$

where we have set
$\hat{\mathbf{t}}\cdot\nabla\mathbf{u}=\partial\mathbf{u}/\partial s$ and
$\hat{\mathbf{b}}\cdot\nabla\mathbf{u}=\partial\mathbf{u}/\partial b$. We can
re-express $\hat{\mathbf{t}}\cdot(\partial\mathbf{u}/\partial s)$ and
$\hat{\mathbf{b}}\cdot(\partial\mathbf{u}/\partial b)$ as follows

$\displaystyle\hat{\mathbf{t}}\cdot\frac{\partial\mathbf{u}}{\partial s}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
s}(\mathbf{u}\cdot\hat{\mathbf{t}})-\mathbf{u}\cdot\frac{\partial\hat{\mathbf{t}}}{\partial
s},$ (68)
$\displaystyle\hat{\mathbf{b}}\cdot\frac{\partial\mathbf{u}}{\partial b}$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
b}(\mathbf{u}\cdot\hat{\mathbf{b}})-\mathbf{u}\cdot\frac{\partial\hat{\mathbf{b}}}{\partial
b}.$ (69)

Substitution of Eq. (50) into
$\mathbf{u}\cdot(\partial\hat{\mathbf{t}}/\partial s)$ and
$\mathbf{u}\cdot(\partial\hat{\mathbf{b}}/\partial s)$ leads to

$\displaystyle\mathbf{u}\cdot\frac{\partial\hat{\mathbf{t}}}{\partial s}$
$\displaystyle=$ $\displaystyle-\kappa_{1}(\mathbf{u}\cdot\hat{\mathbf{s}}),$
(70) $\displaystyle\mathbf{u}\cdot\frac{\partial\hat{\mathbf{b}}}{\partial b}$
$\displaystyle=$ $\displaystyle-\kappa_{2}(\mathbf{u}\cdot\hat{\mathbf{s}}),$
(71)

where, again, we have made use of the fact that
$\hat{\mathbf{t}}\neq\hat{\mathbf{t}}(b)$ and
$\hat{\mathbf{b}}\neq\hat{\mathbf{b}}(s)$. Substitution of Eqs. (70) and (71)
into Eqs. (68) and (69) and the resultant relations into Eq. (67) gives

$\hat{\mathbf{s}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}=-\frac{\partial}{\partial
s}\left[(\mathbf{u}\cdot\hat{\mathbf{t}})+(\mathbf{u}\cdot\hat{\mathbf{b}})\right]-(\kappa_{1}+\kappa_{2})(\mathbf{u}\cdot\hat{\mathbf{s}}).$
(72)

Since
$\hat{\mathbf{s}}\cdot\mathbf{D}\cdot\hat{\mathbf{s}}=2\hat{\mathbf{s}}\cdot\nabla\mathbf{u}\cdot\hat{\mathbf{s}}$,
it follows that

$\hat{\mathbf{s}}\cdot\mathbf{D}\cdot\hat{\mathbf{s}}=-2\frac{\partial}{\partial
s}\left[(\mathbf{u}\cdot\hat{\mathbf{t}})+(\mathbf{u}\cdot\hat{\mathbf{b}})\right]-2(\kappa_{1}+\kappa_{2})(\mathbf{u}\cdot\hat{\mathbf{s}}).$
(73)

Substitution of this equation into Eq. (66) yields the following expression
for the pressure jump

$[[p]]=-\sigma(\kappa_{1}+\kappa_{2})-2[[\mu\left(\frac{\partial}{\partial
s}\left[(\mathbf{u}\cdot\hat{\mathbf{t}})+(\mathbf{u}\cdot\hat{\mathbf{b}})\right]+(\kappa_{1}+\kappa_{2})(\mathbf{u}\cdot\hat{\mathbf{s}})\right)]].$
(74)

## References

  * Ariyapadi et al. (2004) Ariyapadi, S., Balachandar, R. & Berruti, F. 2004 Effect of surfactant on the characteristics of a droplet-laden jet. Chem Eng Process 43 (4), 547 – 553.
  * Batchvarov et al. (2021) Batchvarov, A., Kahouadji, L., Constante-Amores, C. R., Norões Gonçalves, G. F., Shin, S., Chergui, J., Juric, D., Craster, R. V. & Matar, O. K. 2021 Three-dimensional dynamics of falling films in the presence of insoluble surfactants. J. Fluid Mech. 906, A16.
  * Batchvarov et al. (2020) Batchvarov, A., Kahouadji, L., Magnini, M., Constante-Amores, C. R., Craster, R. V ., Shin, S., Chergui, J., Juric, D. & Matar, O. K. 2020 Effect of surfactant on elongated bubbles in capillary tubes at high reynolds number. Phys. Rev. Fluids 5, 093605.
  * Broze & Hussain (1996) Broze, G. & Hussain, F. 1996 Transitions to chaos in a forced jet: intermittency, tangent bifurcations and hysteresis. J. Fluid Mech. 311, 37–71.
  * Brøns et al. (2014) Brøns, M., Thompson, M. C., Leweke, T. & Hourigan, K. 2014 Vorticity generation and conservation for two-dimensional interfaces and boundaries. J. Fluid Mech. 758, 63–93.
  * Constante-Amores et al. (2021a) Constante-Amores, C.R., Kahouadji, L., Batchvarov, A., Shin, S., Chergui, J., Juric, D. & Matar, O.K. 2021a Direct numerical simulations of transient turbulent jets: vortex-interface interactions. J. Fluid Mech. 922, A6.
  * Constante-Amores et al. (2020a) Constante-Amores, C. R., Kahouadji, L., Batchvarov, A., Shin, S., Chergui, J., Juric, D. & Matar, O. K. 2020a Dynamics of retracting surfactant-laden ligaments at intermediate ohnesorge number. Phys. Rev. Fluids 5, 084007.
  * Constante-Amores et al. (2020b) Constante-Amores, C. R., Kahouadji, L., Batchvarov, A., Shin, S., Chergui, J., Juric, D. & Matar, O. K. 2020b Rico and the jets: Direct numerical simulations of turbulent liquid jets. Phys. Rev. Fluids 5, 110501\.
  * Constante-Amores et al. (2021b) Constante-Amores, C. R., Kahouadji, L., Batchvarov, A., Shin, S., Chergui, J., Juric, D. & Matar, O. K. 2021b Dynamics of a surfactant-laden bubble bursting through an interface. J. Fluid Mech. 911, A57.
  * Desjardins & Pitsch (2010) Desjardins, O. & Pitsch, H. 2010 Detailed numerical investigation of turbulent atomization of liquid jets. Atomiz. Sprays 20 (4), 311–336.
  * Eggers (1997) Eggers, J. 1997 Nonlinear dynamics and breakup of free-surface flows. Rev. Mod. Phys. 69 (3), 865–929.
  * Ellis et al. (2001) Ellis, M.C B., Tuck, C.R & Miller, P.C.H 2001 How surface tension of surfactant solutions influences the characteristics of sprays produced by hydraulic nozzles used for pesticide application. Colloids Surf 180 (3), 267 – 276.
  * Fuster & Rossi (2021) Fuster, D. & Rossi, M. 2021 Vortex-interface interactions in two-dimensional flows. Int. J. Multiph 143, 103757\.
  * Herrmann (2010) Herrmann, M. 2010 A parallel eulerian interface tracking/lagrangian point particle multi-scale coupling procedure. J. Comput. Phys 229 (3), 745 – 759.
  * Hoepffner & Paré (2013) Hoepffner, J. & Paré, G. 2013 Recoil of a liquid filament: escape from pinch-off through creation of a vortex ring. J. Fluid Mech. 734 (183–197).
  * Hunt et al. (1988) Hunt, J., Wray, A. & Moin, P. 1988 Eddies, streams, and convergence zones in turbulent flows. Studying Turbulence Using Numerical Simulation Databases 1, 193–208.
  * Ibarra et al. (2020) Ibarra, E., Shaffer, F. & Savaş, O. 2020 On the near-field interfaces of homogeneous and immiscible round turbulent jets. J. Fluid Mech. 889, A4.
  * Ibarra (2017) Ibarra, R. 2017 Horizontal and low-inclination oil-water flow investigations using laser-based diagnostic techniques. Imperial College London PhD Thesis.
  * Jarrahbashi & Sirignano (2014) Jarrahbashi, D. & Sirignano, W. A. 2014 Vorticity dynamics for transient high-pressure liquid injection. Phys. Fluids 26 (10), 101304.
  * Jarrahbashi et al. (2016) Jarrahbashi, D., Sirignano, W. A., Popov, P. P. & Hussain, F. 2016 Early spray development at high gas density: hole, ligament and bridge formations. J. Fluid Mech. 792, 186–231.
  * Kooij et al. (2018) Kooij, S., Sijs, R., Denn, M. M., Villermaux, E. & Bonn, D. 2018 What determines the drop size in sprays? Phys. Rev. X 8, 031019.
  * Longuet-Higgins (1992) Longuet-Higgins, Michael S. 1992 Capillary rollers and bores. J. Fluid Mech. 240, 659–679.
  * Lundgren & Koumoutsakos (1999) Lundgren, T. & Koumoutsakos, P. 1999 On the generation of vorticity at a free surface. J. Fluid Mech. 382, 351–366.
  * Manikantan & Squires (2020) Manikantan, H. & Squires, T. M. 2020 Surfactant dynamics: hidden variables controlling fluid flows. J. Fluid Mech. 892, P1.
  * Marmottant & Villermaux (2004) Marmottant, P. & Villermaux, E. 2004 On spray formation. J. Fluid Mech. 498, 73 – 111.
  * Plateau (1873) Plateau, J. 1873 Experimental and Theoretical Statics of Liquids Subject to Molecular Forces Only, , vol. 1. Gand et Leipzig: F. Clemm.
  * Shin et al. (2018) Shin, S., Chergui, J., Juric, D., Kahouadji, L., K.Matar, O. & Craster, R. V. 2018 A hybrid interface tracking - level set technique for multiphase flow with soluble surfactant. J. of Comp. Phys. 359, 409–435.
  * Sijs & Bonn (2020) Sijs, R. & Bonn, D. 2020 The effect of adjuvants on spray droplet size from hydraulic nozzles. Pest Management Science 76 (10), 3487–3494.
  * Sijs et al. (2021) Sijs, R., Kooij, S., Holterman, H. J., van de Zande, J. & Bonn, D. 2021 Drop size measurement techniques for sprays: Comparison of image analysis, phase doppler particle analysis, and laser diffraction. AIP Advances 11 (1), 015315.
  * da Silva & Métais (2002) da Silva, C. B. & Métais, O. 2002 Vortex control of bifurcating jets: A numerical study. Physics of Fluids 14 (11), 3798–3819.
  * Strickland et al. (2015) Strickland, S. L., Shearer, M. & Daniels, K. E. 2015 Spatio-temporal measurement of surfactant distribution on gravity–capillary waves. J. Fluid Mech. 777, 523–543.
  * Terrington et al. (2020) Terrington, S. J., Hourigan, K. & Thompson, M. C. 2020 The generation and conservation of vorticity: deforming interfaces and boundaries in two-dimensional flows. J. Fluid Mech. 890, A5.
  * Terrington et al. (2021) Terrington, S. J., Hourigan, K. & Thompson, M. C. 2021 The generation and diffusion of vorticity in three-dimensional flows: Lyman’s flux. J. Fluid Mech. 915, A106–6–42.
  * Theodorsen (1952) Theodorsen, T. 1952 Mechanism of turbulence. Proceedings of the Midwestern Conference Fluid Mechanics, 1-19. .
  * Urbin & Métais (1997) Urbin, G & Métais 1997 Large-eddy simulations of three-dimensional spatially-developing round jets. ERCOFTAC Series. Springer 5.
  * Wee et al. (2021) Wee, Hansol, Wagoner, Brayden W., Garg, Vishrut, Kamat, Pritish M. & Basaran, Osman A. 2021 Pinch-off of a surfactant-covered jet. J. Fluid Mech. 908, A38.
  * Williams et al. (2021) Williams, A. G. L., Karapetsas, G., Mamalis, D., Sefiane, K., Matar, O. K. & Valluri, P. 2021 Spreading and retraction dynamics of sessile evaporating droplets comprising volatile binary mixtures. J. Fluid Mech. 907, A22.
  * Wu (1995) Wu, J. Z. 1995 A theory of three-dimensional interfacial vorticity dynamics. Phys. Fluids 7 (10), 2375–2395.
  * Zandian et al. (2018) Zandian, A., Sirignano, W. A. & Hussain, F. 2018 Understanding liquid-jet atomization cascades via vortex dynamics. J. Fluid Mech. 843, 293–354.
  * Zandian et al. (2019) Zandian, A., Sirignano, W. A. & Hussain, F. 2019 Vorticity dynamics in a spatially developing liquid jet inside a co-flowing gas. J. Fluid Mech. 877, 429–470.