1. Introduction
The breakup of a dispersed fluid in a stagnant phase is a classical problem in multiphase flows, as it encompasses a multitude of interfacial singularities, which exemplify situations wherein the interface undergoes topological transitions, e.g. liquid threads breakup into drops, and merging of drops/bubbles due to coalescence. These transitions involve the development of singularities where interfacial distances vanish and velocity fields diverge, and the system is controlled by a combination of capillary, inertial and viscous forces. The complex topological nature of the jet phenomenon has fascinated the scientific community for decades which has led to numerous, comprehensive reviews; see, for example, Lin & Reitz (Reference Lin and Reitz1998), Eggers & Villermaux (Reference Eggers and Villermaux2008) and Lasheras & Hopfinger (Reference Lasheras and Hopfinger2000).
Jet breakup is influenced by a range of multi-scale physics; these include the interaction of turbulence with interfaces (creating cascades of motion featuring a large separation of scales), capillarity, potentially complex rheology, the presence of fields (e.g. gravitational, electromagnetic), as well as heat transfer and phase change. Reitz & Bracco (Reference Reitz and Bracco1986) proposed four jet breakup regimes depending on the appearance of the jet far downstream from the injection point. In the ‘Rayleigh regime’ the onset of breakup is the Rayleigh–Plateau instability, in which the growth of the linear modes on the jet surface leads to the formation of large droplets with respect to the jet nozzle. In the ‘first/second wind-induced’ regime the resulting droplets have roughly the same scale or smaller than the jet nozzle. Finally, in the ‘atomisation regime’ the generated turbulence helps to create spatio-temporal chaos resulting in droplets which are up to four orders of magnitude smaller than the size of the injection nozzle. The physics of droplet generation in turbulent jets remains partially understood despite the significant scientific attention it has received over the years (Dombrowski, Fraser & Newitt Reference Dombrowski, Fraser and Newitt1954; Hoyt & Taylor Reference Hoyt and Taylor1977; Lasheras & Hopfinger Reference Lasheras and Hopfinger2000; Marmottant & Villermaux Reference Marmottant and Villermaux2004; Villermaux, Marmottant & Duplat Reference Villermaux, Marmottant and Duplat2004; Eggers & Villermaux Reference Eggers and Villermaux2008). The rest of this introduction aims to provide an up-to-date summary of the experimental and numerical efforts to study the jet breakup phenomenon.
The ground-breaking experiments of Hoyt & Taylor (Reference Hoyt and Taylor1977) and Taylor & Hoyt (Reference Taylor and Hoyt1983) showed the complex topological features of the surface waves formed on the jet. They revealed that these waves are responsible for the transition from laminar to turbulent flow, and found underlying similarities of the occurring instabilities to inviscid linear theory. Marmottant & Villermaux (Reference Marmottant and Villermaux2004) with their pioneering experiments scrutinised the various stages of the jet dynamics: from the growth of linear modes (through a Kelvin–Helmholtz instability) that characterises the early time dynamics, to the development of nonlinearities leading to ‘primary’ and subsequent ‘secondary’ breakup events (through long filament pinchoff modulated by a Rayleigh–Taylor instability), and the formation of a cascade of droplet sizes. These authors found that mean droplet size is proportional to the wavelength selected during the Rayleigh–Taylor instability (also observed by Varga, Lasheras & Hoepfinger Reference Varga, Lasheras and Hoepfinger2003). In the same premise, Kooij et al. (Reference Kooij, Sijs, Denn, Villermaux and Bonn2018) showed that the droplet size distribution is also a function of the nozzle geometry and the surrounding pressure. More recently, Ibarra, Shaffer & Savaş (Reference Ibarra, Shaffer and Savaş2020) presented the results of an experimental study of the spatial evolution of turbulent immiscible liquid-liquid jets; however, a detailed account of the spatio-temporal development, and critical mechanisms leading to droplet generation remains outstanding.
The multi-scale nature of the flow, and the complex interfacial topology complicate the experimental scrutiny of the different physical mechanisms occurring across the scales. Thus, elucidating the fundamental physics of this problem has also relied on high-fidelity simulations exemplified by the work of Bianchi et al. (Reference Bianchi, Pelloni, Toninel, Scardovelli, Leboissetier and Zaleski2005, Reference Bianchi, Minelli, Scardovelli and Zaleski2007), Ménard, Tanguy & Berlemont (Reference Ménard, Tanguy and Berlemont2007), Desjardins, Moureau & Pitsch (Reference Desjardins, Moureau and Pitsch2008), Gorokhovski & Herrmann (Reference Gorokhovski and Herrmann2008), Desjardins & Pitsch (Reference Desjardins and Pitsch2010), Shinjo & Umemura (Reference Shinjo and Umemura2010), Herrmann (Reference Herrmann2011), Chenadec (Reference Chenadec2012), Desjardins et al. (Reference Desjardins, McCaslin, Owkes and Brady2013), Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016), Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017), Agbaglah, Chiodi & Desjardins (Reference Agbaglah, Chiodi and Desjardins2017), Zandian, Sirignano & Hussain (Reference Zandian, Sirignano and Hussain2018), Ling et al. (Reference Ling, Fuster, Tryggvason and Zaleski2019) and Zandian, Sirignano & Hussain (Reference Zandian, Sirignano and Hussain2019). The first step towards the use of numerical simulations to understand the physical mechanisms at play was conducted by Desjardins & Pitsch (Reference Desjardins and Pitsch2010), who were able to identify a sequence of essential steps during the spatio/temporal interfacial development of a planar liquid jet segment: the formation of initial corrugations on the surface, followed by the development of ligaments whose capillary instability leads to droplet formation. They also showed that the early interfacial corrugations are a consequence of the turbulent eddies which carry enough kinetic energy to overcome capillary forces.
Using a similar approach to Desjardins & Pitsch (Reference Desjardins and Pitsch2010), but for a cylindrical liquid jet segment surrounded by an outer gas phase, Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) provided a comprehensive study of the flow structures in terms of the vortex-surface interaction. They showed that ‘hole formation’ of a liquid sheet is an essential requirement to trigger the formation of droplets, and the thinning of the liquid sheet is driven by the superposition of hairpin vortices near the interface rather than by capillarity action. Similar findings have been reported for a planar liquid jet segment surrounded by an outer gas phase (Zandian et al. Reference Zandian, Sirignano and Hussain2018), and for the transient dynamics of a cylindrical liquid jet surrounded by a coaxial air phase (Zandian et al. Reference Zandian, Sirignano and Hussain2019). Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017, Reference Ling, Fuster, Tryggvason and Zaleski2019) performed simulations of a two-phase mixing layer between parallel gas and liquid streams to investigate the interfacial dynamics, and the statistics for the multiphase turbulence. They also observed that the formation of ligaments, and subsequently formed droplets, are triggered by the hole-induced perforation of the liquid sheets.
In the previous numerical studies, interface capturing capabilities were used to account for the surface tension forces in the absence of intermolecular forces (i.e. disjoining pressure). For static sheets, the disjoining pressure is neglected as the minimum computational cell is larger than the film sheet thickness in which the intermolecular forces will drive its perforation. Therefore, the hole formation is an outcome of the numerical cut-off interfacial length scale, i.e. minimum mesh size (
$O(10^{-6}$) m). Nevertheless, recent experiments (Marston et al. Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016; Kooij et al. Reference Kooij, Sijs, Denn, Villermaux and Bonn2018; Néel & Villermaux Reference Néel and Villermaux2018) demonstrate the existence of hole formation in dynamic sheets with a characteristic film thickness on the order of microns.
In this study we aim to provide a comprehensive explanation of the physical mechanisms governing the interfacial dynamics of turbulent jets focusing on the less well-studied liquid-liquid systems. We will perform high-resolution three-dimensional direct numerical simulations (DNS) using a hybrid interface-tracking/level-set approach to resolve the interfacial dynamics. We will demonstrate how the interaction between the vortical structures, which accompany the development of the flow, and the interface influence the mechanisms underlying droplet generation over a wide range of Reynolds and Weber numbers.
The rest of this paper is organised as follows. In § 2 we present the governing equations along with the numerical technique used to carry out the simulations. In § 3 we begin with the presentation of a regime map in 
$Re$–
$We$ space, classifying jet spatio-temporal development, followed by an in-depth discussion of the vortex-surface interactions linked to the topological changes; in addition, we elucidate the role of hole formation as a precursor to droplet generation. Finally, concluding remarks are provided in § 4.
2. Problem formulation and numerical techniques
Figure 1 shows a three-dimensional representation of the interface highlighting several specific features that are discussed in detail in the present work: Kelvin–Helmholtz (KH) waves close to the injection nozzle, outer lobes formed on the main body of the jet, a leading edge mushroom-like structure and droplets resulting from the atomisation process. The numerical framework is based on solving the two-phase incompressible Navier–Stokes equations in a three-dimensional Cartesian domain 
$x = (x, y, z)$. The surface tension force in the momentum equation is treated by using a hybrid front-tracking/level-set method presented by Shin & Juric (Reference Shin and Juric2002). The governing flow equations in the ‘one-fluid’ formulation are described by
\begin{equation} \left.\begin{gathered} {\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} =0},\\ \displaystyle{\rho\left(\frac{\partial \boldsymbol{u}}{\partial t} + \boldsymbol{u} \boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u}\right) ={-} \boldsymbol{\nabla} p + \boldsymbol{\nabla}\boldsymbol{\cdot} \mu ({\boldsymbol{\nabla} \boldsymbol{u}}+{\boldsymbol{\nabla} \boldsymbol{u}}^\textrm{T}) + \boldsymbol{F}_s}, \end{gathered}\right\} \end{equation}
where 
$t$, 
$\boldsymbol {u}$, 
$p$ and 
$\boldsymbol {F}_s$ stand for time, velocity, pressure and the surface tension force, respectively. The density and viscosity are expressed using the formulation
\begin{equation} \left.\begin{gathered} {\rho=\rho_{{so}} + \left(\rho_{w} - \rho_{{so}}\right) \mathcal{H}\left(\boldsymbol{x},t\right) },\\ {\mu=\mu_{{so}} + \left(\mu_{w} - \mu_{{so}}\right) \mathcal{H}\left( \boldsymbol{x},t \right)}. \end{gathered}\right\} \end{equation}
wherein 
$\mathcal {H}( \boldsymbol {x},t)$ represents a smoothed Heaviside function, which is zero in the dispersed phase (water) and unity in the stagnant phase (silicone oil), while the subscripts ‘
$w$’ and ‘
$so$’ refer to the individual phases. The surface tension force 
$\boldsymbol {F}_s$ is defined by using the hybrid formulation as in Shin & Juric (Reference Shin and Juric2009) and Shin, Chergui & Juric (Reference Shin, Chergui and Juric2017),
\begin{equation} \boldsymbol{F}_s=\sigma \kappa_\mathcal{H} \boldsymbol{\nabla} \mathcal{H}, \end{equation}
where 
$\sigma$ refers to surface tension, and 
$\kappa _\mathcal {H}$ is twice the mean interface curvature calculated from the Eulerian grid by using
\begin{equation} \kappa_\mathcal{H} = \frac{\boldsymbol{F}_L\cdot \boldsymbol{G}}{\sigma \boldsymbol{G}\cdot \boldsymbol{G}}, \end{equation}where
\begin{equation} \boldsymbol{F}_{L} = \int_{\varGamma (t) } \sigma \kappa \boldsymbol{n} \delta (\boldsymbol{x}-\boldsymbol{x}_f) \,{\textrm{d}}s,\quad \boldsymbol{G} = \int_{\varGamma (t) } \boldsymbol{n} \delta (\boldsymbol{x}-\boldsymbol{x}_f),{\textrm{d}}s. \end{equation}
Here, 
$\boldsymbol {x}_f$ is the parametrisation of the interface, 
${\varGamma (t) }$, and 
$\delta (\boldsymbol {x}-\boldsymbol {x}_f)$ is a Dirac delta function which vanishes everywhere except at the interface; 
$\boldsymbol {n}$ is the outward-pointing unit normal vector to the interface and 
$ds$ is the length of the surface element.
Injection of a turbulent water jet into a stagnant silicone oil phase. Three-dimensional representation of the interfacial shape for case 4 in table 2 (i.e. 
$Re=6530$ and 
$We=303$) at 
$t=28.97$ with definitions of particular regions and specific features discussed in this work.
The incompressible Navier–Stokes equations (2.1) are solved by using a second-order finite-difference method on a staggered grid (Harlow & Welch Reference Harlow and Welch1965; Temam Reference Temam1968). The computational domain is then discretised by a fixed, regular, Eulerian grid, and the spatial derivatives are approximated by standard centred difference discretisation, except for the nonlinear term, which makes use of a second-order essentially non-oscillatory scheme (Sussman, Smereka & Osher Reference Sussman, Smereka and Osher1994). As for the viscous term, we use a second-order centred difference scheme. We have used the projection method to handle the incompressibility condition combined with a multigrid iterative method for solving the elliptic pressure Poisson equation (Chorin Reference Chorin1968; Kwak, do, Lee Reference Kwak, do, Lee2004). The numerical method uses an additional adaptive Lagrangian grid based on a hybrid front-tracking/level-set method (Shin & Juric Reference Shin and Juric2007; Shin et al. Reference Shin, Chergui and Juric2017) to track the interface location. The geometrical information, pressure and velocity variables are exchanged between the adaptive Lagrangian mesh and the fixed Eulerian grid following the immersed boundary method of Peskin (Reference Peskin1977). The physical elements that form the Lagrangian interface are advected according to 
$\rm {d}\boldsymbol {x}_f/ \rm {d}t =\boldsymbol {V}$, where 
$\boldsymbol {V}$ is the interface velocity interpolated numerically from the fixed Eulerian grid, using a second-order Runge–Kutta method. Finally, the numerical method uses a domain-decomposition technique for its parallelisation and message passing interface for exchanging information between adjacent subdomains. More information regarding the full implementation of the numerical method can be found in Shin et al. (Reference Shin, Chergui and Juric2017, Reference Shin, Chergui, Juric, Kahouadji, Matar and Craster2018).
2.1. Numerical configuration and physical parameters
Figure 1 shows the three-dimensional computational domain, which is a rectangular box of size 
$20D \times 4D \times 4D$, where 
$D$ stands for the inner diameter of the nozzle (e.g. 4 mm). The jet is produced when water leaves the cylindrical nozzle to enter progressively into the stagnant silicone oil. The physical properties of the fluids are given in table 1. An inflow boundary condition is applied to the nozzle on the left of the domain, i.e. at 
$(x = 0)$, which follows a simplified power-law turbulent velocity profile,
\begin{equation} u\left(r,t\right) = \frac{15}{14} {U} \left(1-\left(\frac{r}{D/2}\right)^{28} \right) \left(1+A \sin\left(2{\rm \pi} f t\right) \right). \end{equation}
Here, 
${U}$, 
$A$ and 
$f$ stand for the average injection velocity, amplitude and frequency, respectively, of the external pulsatile perturbation. The radial distance, 
$r$, within the jet measured from its centreline 
$(y_0, z_0)$ is 
$r=\sqrt {(y-y_0)^{2} + (z-z_0)^{2}}$. The values for 
$A$ and 
$f$ are informed by the previous work of Ling, Zaleski & Scardovelli (Reference Ling, Zaleski and Scardovelli2015) (e.g. 
$A=0.05$ m s
$^{-1}$ and 
$f=20$ Hz), and the same values have been used throughout the entire paper.
Density and viscosity of the fluids used throughout this work (Ibarra Reference Ibarra2017).
The numerical set-up closely follows other computational studies; for example, we impose a free boundary condition on the walls of the computational domain to let the fluid freely enter or leave the boundaries (Taub et al. Reference Taub, Lee, Balachandar and Sherif2013; Ling et al. Reference Ling, Zaleski and Scardovelli2015, Reference Ling, Fuster, Tryggvason and Zaleski2019). A pressure outflow boundary condition is applied on the right surface of the domain to allow the fluid to exit the domain. The solid nozzle is treated as a no-slip surface (Asadi, Asgharzadeh & Borazjani Reference Asadi, Asgharzadeh and Borazjani2018).
The distance 
$\boldsymbol {x}$, velocity 
$\boldsymbol {u}$, time 
$t$ and pressure 
$p$ in (2.1) are rendered dimensionless using the following characteristic scales, 
$D$, 
$U$, 
$D/{U}$ and 
$\rho _{w} U^{2}$, respectively. Hence, the dimensionless control parameters governing the phenomena we will study are given by
\begin{equation} Re =\frac{\rho_wU D}{\mu_w} , \quad We =\frac{\rho_w U^{2} D}{\sigma}, \end{equation}
where 
$Re$ stands for the Reynolds number (i.e. the ratio of inertial to viscous forces), and 
$We$ represents the Weber number (i.e. the ratio of inertial to capillary forces). All the variables appearing in the equations and boundary conditions are rendered dimensionless using the aforementioned scalings, unless stated otherwise.
The first part of the results section corresponds to a discussion of a regime map of jet dynamics in the 
$Re$–
$We$ space. In this instance, we have used the M2 mesh (see table 5 in Appendix A) to perform the fully three-dimensional simulations. The selection of this mesh is dictated by the need to map out parameter space relatively rapidly before focusing on elucidating the details of the dynamics using the M3 mesh, which provides higher resolution, for four cases; the 
$Re$–
$We$ combinations for these cases are listed in table 2. Information regarding the mesh-refinement study, resolution considerations and the validation for this work are detailed in Appendix A.
The Reynolds-Weber number combinations for the four cases studied in detail with the M3 mesh (see table 5 in Appendix A).
3. Results
3.1. Interfacial dynamics: phase diagram in 
$Re$–$We$ space
We start the discussion of the results by presenting a phenomenological picture of the interfacial dynamics in a phase diagram in 
$Re$–
$We$ space. The Reynolds and Weber numbers range between 
$10^{3}$–
$10^{4}$ and 7–900, respectively. Figure 2 shows the regime map in terms of the interfacial dynamics predicted from our numerical simulations. Several features emerge from this figure, which we have divided into four phenomenological regions based on the appearance of the interfacial structures. We aim to quantify the different jet behaviours by close inspection of the interplay between the vorticity field, 
$\boldsymbol {\omega } =\boldsymbol {\nabla } \times \boldsymbol {u}$, and the interface. Such in-depth analysis is carried out following the introduction of each region, utilising the higher resolution M3 mesh examples, as shown in table 2.
Regime map of the phenomenological interfacial dynamics in the 
$Re$–
$We$ space using the 
$M2$ mesh (see table 5 in Appendix A for mesh details). Four different regimes and their boundaries are identified. A snapshot of the flow corresponding to the three-dimensional representation of the interface for each simulated point (i.e. black marker) is shown. The red squares refer to the cases presented in table 2.
Region ‘A’ in figure 2 is characterised by low Reynolds and Weber numbers and defined by the dominance of the capillary over inertial forces, yielding an axisymmetric behaviour of the jet. For small 
$We$, the dominance of the surface tension forces results in the entrapment of the initial vortex head inside the discharged liquid edge. Then, the formation of an interfacial leading mushroom-like structure is not observed. Figure 3 shows the interplay between the azimuthal and streamwise vorticity components, 
$\omega _{\theta }$ and 
$\omega _x$, respectively, for case 1 (see table 2). We observe that 
$\omega _{\theta }$ exceeds 
$\omega _{x}$ by two orders of magnitude for the entirety of the jet, which explains the lack of deformation of the jet core, and its axisymmetric shape.
Two-dimensional representation of the interfacial location together with 
$\omega _{\theta }$ (a) and 
$\omega _{x}$ (b) in the 
$x$–
$z$ plane (
$y=2$) for case 1 in table 2 (
$Re=1000$ and 
$We=7$) at 
$t=24.38$. The colour represents the respective vorticity field, where appropriate scales are shown in each panel.
Region ‘B’ is defined by low Reynolds and high Weber numbers. As shown in figure 2, the snapshots of the interface in this region of parameter space reveal the development of interfacial waves as well as the formation of a mushroom-like structure at the jet leading edge. This is due to the rolling of the boundary layer at the edge of the solid nozzle as a result of the initial vortex ring (more details in § 3.2). Figure 4 shows that although the magnitude of 
$\omega _{\theta }$ exceeds that of 
$\omega _x$ by almost two orders of magnitude for case 2, the relative significance of 
$\omega _x$ has increased in comparison with case 1. This is correlated with the corrugations which are observed on the main body of the jet in case 2 that appear to be largely absent in case 1. These corrugations are related to the development of KH instabilities that arise from the velocity difference between the injected jet and the surrounding, initially stagnant phase. As shown in figure 4, the mushroom structure is also accompanied by the formation of droplets within it. More information regarding the mechanisms which induce droplet formation is provided in § 3.3.
Two-dimensional representation of the interfacial location together with 
$\omega _{\theta }$ (a) and 
$\omega _{x}$ (b) in the 
$x$–
$z$ plane (
$y=2$) for case 2 in table 2 (
$Re=1000$ and 
$We=100$) at 
$t=26.87$. The colour represents the respective vorticity fields, where appropriate scales are shown in each panel.
Region ‘C’ is characterised by low Weber and high Reynolds numbers. To elucidate the mechanisms at play in this region, we refer to case 3 (see table 2), as shown in figure 5. It is seen clearly that the magnitude of 
$\omega _{\theta }$ is an order of magnitude larger than that of 
$\omega _{x}$. The higher levels of inertia in case 3, in comparison to case 2, leads to formation of pronounced KH instabilities with vortical rings close to the interface, yielding spanwise interfacial deformations. Close inspection of figure 5 reveals that the mushroom structure, which is also present in case 3, undergoes significant deformation leading to the formation of a toroidal sheet, which envelopes that main body of the jet.
Two-dimensional representation of the interfacial location together with 
$\omega _{\theta }$ (a) and 
$\omega _{x}$ (b) in the 
$x$–
$z$ plane (
$y=2$) for case 3 in table 2 (
$Re=3260$ and 
$We=75$) at 
$t=24.65$. The colour represents the respective vorticity fields, where appropriate scales are shown in each panel.
In region ‘D’ the flow is characterised by both high Reynolds and Weber numbers. Here, it is seen from the snapshots of the interface in figure 2 that the interfacial dynamics in this region of 
$Re$–
$We$ space are the most complex. The mushroom structures suffer severe deformation as does the main body of the jet which also features the formation of lobes arising from the KH-induced corrugations. It is also evident that flow in region D is also accompanied by droplet generation. In § 3.2 below we focus on case 4 in figure 2 and provide an extensive explanation of the mechanisms linking droplet formation to vortex-interface interaction.
With the purpose of providing a better understanding of the vortex-surface interaction for the development of the three-dimensional instabilities, we have analysed the rate of change of the vorticity production by taking the curl of the momentum equation (i.e. (2.1)), which can be written as
\begin{equation} \frac{{\textrm{D}}\boldsymbol{\omega}}{{\textrm{D}}x}=\left( \boldsymbol{\omega} \boldsymbol{\cdot} \boldsymbol{\nabla} \right )\boldsymbol{u} + \boldsymbol{\nabla} \times\left ( \frac{\boldsymbol{\nabla}\boldsymbol{\cdot} \boldsymbol{ \tau }}{\rho} \right ) +\frac{1}{\rho^{2}} \boldsymbol{\nabla} \rho \times \boldsymbol{\nabla} P + \boldsymbol{\nabla} \times \left (\frac{\boldsymbol{F}_s}{\rho}\right), \end{equation}
where 
$\boldsymbol {\omega }$ stands for the vorticity and 
$\boldsymbol { \tau }$ represents the viscous stress tensor. The right-hand side represents the vortex stretching, the vorticity generation as a result of the viscous diffusion, the effect of density variation (‘baroclinic torque’) and the surface tension forces, respectively. Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) have shown that the baroclinic term is responsible for the three-dimensional instabilities for large density ratios 
$O(10^{-2})$, meanwhile, at lower density ratios 
$O(10^{-1})$, the streamwise vorticity generation is dominated by the azimuthal tilting and radial tilting of the vortex rings. For our study, the contribution of the baroclinic term to the vorticity generation is negligible as the density ratio is of the same order of magnitude (i.e. 
${\sim }O(10^{0})$).
Table 3 collects the results of a simple dimensional analysis to study the dominant terms of the vorticity generation equation for the selected cases shown in table 2. On this basis, the viscous diffusion term scales as 
$\mu _w U/\rho {\rm \Delta} x^{3}$, and the surface tension term scales as 
$\sigma \kappa / \rho {\rm \Delta} x^{2}$ (similar to what was presented by Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016; Zandian et al. Reference Zandian, Sirignano and Hussain2019). Finally, the vorticity stretching term 
$( \boldsymbol {\omega } \boldsymbol {\cdot } \boldsymbol {\nabla }) \boldsymbol {u}$ has been computed numerically. For regions A–C, the vortex stretching term is negligible in comparison to the other terms, which explains the non-development of the three-dimensional interfacial instabilities on the surface of the jet core. Those regions are characterised by a balance/competition between the surface tension and viscous terms giving rise to the different interfacial instabilities explained above. Finally, in the inertia-dominated region the three terms on the right-hand side of the vortex generation equation are balanced, and their competition determines the interfacial topology. Thus, the vortex stretching term from the vorticity equation plays the major role in the development of the three-dimensional interfacial destabilisation for turbulent jets with density ratios of the same order of magnitude. Similar conclusions were drawn by Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016).
Scalings for the terms of the vorticity transport equation (non-dimensional values).
3.2. Interfacial dynamics explained through vortex-surface interaction
This section focuses on the inertia-controlled region D and case 4 as a principal example of the flow in this region. In figure 6 we examine the spatio-temporal interfacial dynamics overlaid with the magnitude of 
$\omega _\theta$ and 
$\omega _x$ in the 
$x$–
$z$ plane. During the early stages of the injection, a vortex ring is initially formed by the rolling of the boundary layer at the edge of the solid surface. The physical mechanism which results in the formation of the leading vortex ring is in agreement with previous work, such as Gharib, Rambod & Shariff (Reference Gharib, Rambod and Shariff1998), Marugán-Cruz, Rodríguez-Rodríguez & Martínez-Bazán (Reference Marugán-Cruz, Rodríguez-Rodríguez and Martínez-Bazán2013) and Asadi et al. (Reference Asadi, Asgharzadeh and Borazjani2018). Once the jet enters the domain, the rotation of the head vortex drives the entrainment of the stagnant phase as it is transported radially outwards. The initial vortex has a profound effect on the interfacial dynamics with the formation of a mushroom-like structure. As time evolves, the entrainment and formation of a toroidal liquid sheet of stagnant phase are observed inside the mushroom structure (see, for example, figure 6a,b). Upstream of the mushroom, it is seen that the KH instability develops (see figure 6c) amplified by the pulsatile injection into waves which cause a local adverse pressure gradient by virtue of the local interfacial curvature (see figure 6d). With increasing time, we observe the formation of outer lobes as a result of the entrainment of the stagnant phase in the jet core (see figure 6e).
Two-dimensional representation of the spatio-temporal evolution of the interface for case 4 in table 2 (
$Re=6530$ and 
$We=303$) together with 
$\omega _{\theta }$ (top panel of each subfigure) and 
$\omega _{x}$ (bottom panel of each subfigure) in the 
$x$–
$z$ plane (
$y=2$) for the dimensionless times shown in each panel. The colour represents the respective vorticity fields, where appropriate scales are shown in each panel.
Figure 6 also highlights the spatio-temporal development of the azimuthal and streamwise components of the vorticity. At early times, vorticity generation coincides with the velocity boundary layer attached to the interface, which corresponds to strong tangential flow near the interface. The roll-up of the shear layers, and subsequently, the vortex roll-up gives rise to the formation of KH vortex rings close to the interface (see figure 6c). As time increases, the vorticity boundary layer is convected towards the jet core (see figure 6e). Additionally, vorticity dissipates strongly in the stagnant phase due to the damping effect of the viscosity.
Attention is now turned towards the competition between 
$\omega _{\theta }$ and 
$\omega _{x}$. In the region adjacent to the injection point, the azimuthal vorticity component dominates over its streamwise counterpart by two orders of magnitude; this dominance is reflected by the fact that the KH vortex rings shown in figure 6(c) are essentially quasi-axisymmetric. As the flow develops downstream (see figure 6d,e), 
$\omega _{x}$ becomes comparable in magnitude to 
$\omega _{\theta }$ leading to streamwise stretching of the KH vortex rings to form hairpin-shaped vortices. Hairpin vortical or ‘horseshoe’ structures were proposed by Theodorsen (Reference Theodorsen1952) to describe the features of turbulence dynamics related to the existence of a shear or boundary layer near a wall. A hairpin vortex is made from a ‘vortex head’ which is the arched region farther away from the boundary layer. The head is connected to the free surface or wall by its ‘vortex legs’. The hairpin-head orientation results from a balance between the effect of shear flow and the local velocity. Interestingly, the streamwise alignment of hairpin vortices near the surface (see figure 6e, f) has been reported previously by Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) and Zandian et al. (Reference Zandian, Sirignano and Hussain2018) for coaxial round and planar jets, respectively. In spite of the absence of a coaxial phase in our case, we observe the same phenomenon. Additionally, we observe the successive alignment of opposite-signed vortices, where the angle between each vortex pair is determined by the induction effect of the opposite-signed ‘neighbour’ vortex. A similar arrangement of vortices has been reported previously by Jeong et al. (Reference Jeong, Hussain, Schoppa and Kim1997) for boundary layer turbulence, and Davoust, Jacquin & Leclaire (Reference Davoust, Jacquin and Leclaire2012) for compressible homogeneous jets. Further details about the specific alignment of the vortical structure is provided below.
Figure 7(a) presents a three-dimensional representation of the jet at 
$t=28.97$ alongside five spatial locations which will be used to show the evolution of vorticity and velocity across the jet core. The specific locations were chosen as follows: panel ‘A’ is representative of the dynamics near the injection point, where there is little disturbance to the interface; panel ‘B’ displays the cross-section of the jet in a jet ring (i.e. the part of the interface where we observe the KH-waves crest), whereas panel ‘C’ depicts a jet braid (i.e. the KH-waves trough); panels ‘D’ and ‘E’ portray the jet core dynamics where lobe formation is present. Figure 7(b) shows the velocity and vorticity line profiles at the 
$y=2$ plane alongside the streamwise vorticity in the 
$y$–
$z$ plane for each of those locations. Near the injection point (location ‘A’), 
$\omega _{\theta } \gt \gt \omega _{x}$, which is in agreement with the lack of deformation of the vortical structures in the streamwise direction (see also figure 6f). The tangential motion of the fluid close to the surface results in the emergence of the velocity boundary layer. In turn, this causes the formation of two large peaks of opposite signs for 
$\omega _{\theta }$. As the flow evolves downstream, 
$\omega _{\theta }$ still dominates the physics, although, its value has been reduced in favour of an increase in 
$\omega _{x}$, and consequently, the local formation of corrugations on the jet core (location ‘B’). Additionally, the 
$\omega _{\theta }$ peaks have widened in their base due to the growth of the velocity boundary layer as the jet moves downstream (Liepmann & Gharib Reference Liepmann and Gharib1992). Further downstream (location ‘D’ and ‘E’), the interface has lost its cylindrical shape due to the entrainment of the stagnant phase as 
$\omega _{x}$ becomes comparable in magnitude to 
$\omega _{\theta }$. Several peaks can be seen in the vorticity profiles as the probe line passes through several lobes. By inspecting the instantaneous velocity fields in locations ‘D’ and ‘E’, we observe that the main component of velocity is associated with the direction of injection, where there is a decay of its mean centreline value as flow evolves downstream, which is in agreement with Pope (Reference Pope2000) for single-phase jets. The velocity decay leads to a reduction in the velocity difference between the injected and stagnant phases, which attenuates the interfacial shearing, supporting the development of the vortical rings close to the surface of the jet. Notice that the 
$y$ and 
$z$ velocity components are positive or negative close to the surface, showing the entrainment of injected fluid towards the stagnant viscous phase or vice versa.
Three-dimensional representation of the interface for case 4 (
$Re=6530$ and 
$We=303$) at 
$t=28.97$, showing the spatial locations ‘A’–‘E’, (a). The streamwise location of the ‘A’–‘E’ probe lines correspond to 
$x=(1.50, 2.50, 2.75, 4.25, 5.25)$, respectively. (b) Velocity (a) and vorticity profiles (b) in the 
$y=2$ plane for each probe location. (c) Streamwise vorticity in the 
$y$–
$z$ plane for each sampling location. The arrows show examples of identified hairpin vortex legs. Solid and dashed arrows correspond to inner and outer hairpin vortex legs, respectively. The colour represents the streamwise vorticity field, 
$\omega _x$.
Next our attention is turned towards the physical mechanisms which lead to the alignment of the hairpin-like vortical structures along the free surface. Figure 7(c) shows cross-section panels with 
$\omega _{x}$ of the jet at the streamwise locations shown in figure 7(a). As depicted in panel ‘A’, 
$\omega _{x}$ is mainly confined inside of the injected stream, where hairpin vortex legs are observed (see first panel of figure 7c). In panel ‘B’ the outer layer of vorticity comes from the braid located immediately downstream, which starts to roll-up around the ring. Further analysis of panel ‘B’ shows the existence of additional pairs of vortex legs which are 
$180^{\circ }$ out-of-phase with respect to the inner layer, and subsequently, they correspond to a second hairpin vortex layer with opposite direction (e.g. upstream direction). Similarly, in panel ‘C’ we observe that the inner layer of the vorticity comes from the ring located upstream, whereas the outer layer comes from the ring located immediately downstream. Therefore, we can conclude that the alignment of the hairpin vortices is a result of the existence of a vortex-induction mechanism which causes the reorientation of vortical structure within the ring and braid skeleton. This phenomenon is in agreement with Brancher, Chomaz & Huerre (Reference Brancher, Chomaz and Huerre1994) for homogeneous jets, Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016) for the coaxial atomisation of a round liquid jet, Zandian et al. (Reference Zandian, Sirignano and Hussain2018) for coaxial atomisation of planar jets and Bernal & Roshko (Reference Bernal and Roshko1986) for plane mixing layers. The analysis between panels ‘B’ and ‘C’ can be extended to other locations of the jet where the core has undergone further interfacial development. For example, in panels ‘D’ and ‘E’ we observe the same distribution of inner and outer layer of hairpin vortex legs. The alignment of vortical structures has a detrimental effect on the interfacial dynamics, which will be shown in § 3.3.
Additionally, we have used the Q-criterion to present a three-dimensional visualisation of the vortices in the present study. The Q-criterion was described by Hunt, Wray & Moin (Reference Hunt, Wray and Moin1988) as a quantity which measures the dominance of vorticity 
$\omega$ over strain 
$s$, 
$Q=1/2 ( \| \omega \|^{2} - \| s \|^{2} )$. Figure 8 shows the spatial development of the coherent structures for 
$Q=0.1$. Near the injection point, quasi-axisymmetric KH vortex rings are located close to the surface (a magnified view is shown). As explained previously, when 
$\omega _{x}$ and 
$\omega _{\theta }$ are of comparable magnitude, the KH vortices are stretched downstream or upstream. The topological shape of these vortical structures resembles the instantaneous hairpin-like vortical structures reported in experiments and numerical simulations by Head & Bandyopadhyay (Reference Head and Bandyopadhyay1981), Zhou et al. (Reference Zhou, Adrian, Balachandar and Kendall1999) and Zandian et al. (Reference Zandian, Sirignano and Hussain2018). Outer hairpin vortices (a magnified view is also shown) are observed clearly and the inner hairpin vortices not as clearly as they are localised underneath the interface. Further downstream, we observe the ‘head vortex’ covering the mushroom-like structure. Inside of this structure, the vortices are unstable and break down.
Illustration of the coherent vortical structures close to the interface for case 4 (
$Re=6530$ and 
$We=303$) at 
$t=28.97$. The coherent vortical structures are visualised by the Q-criterion with a value of 
$Q=0.1$, where the colour represents the streamwise vorticity field, 
$\omega _w$. Magnified views of the KH vortex rings near the injection point, and the hairpin vortical structures in the jet core are also shown.
To conclude, we have identified three different zones in the transient flow field which are associated with different vortical dynamics. We present the different regions at 
$t=28.97$ (see figure 8). Zone 1 starts from the injection point and extends up to the deformation of the free surface (
$x \sim 1.8$). This zone is characterised by the dominance of 
$\omega _{\theta }$. Zone 2 extends from 
$x \sim 1.8$ up to 
$x \sim 11.5$ (e.g. behind the mushroom-like structure). This region is characterised by the interfacial deformation of the jet core by KH vortex rings and their posterior deformation by the competition between 
$\omega _{x}$ and 
$\omega _{\theta }$. As the flow moves downstream, 
$\omega _{x}$ becomes responsible for the entrainment of the stagnant phase to form interfacial lobes (this agrees with Liepmann & Gharib (Reference Liepmann and Gharib1992) for homogeneous jets). In this region vortex rings pair up and merge together. Additionally, the ending of the visible potential core of the jet and the beginning of the mixing region is also observed. Zone 3 expands from the end of zone 2 up to the leading edge of the mushroom-like structure. A large vortex ring dominates the dynamics in this region (so-called ‘cap vortex’, Zandian et al. Reference Zandian, Sirignano and Hussain2019), wherein its interaction with upstream vortex structures via means of velocity induction leads to the formation of the neck, connecting the mushroom-like structure with the cylindrical body of the jet. The narrowing of the neck as a result of the streamwise convection of vortical structures is in agreement with the work of Asadi et al. (Reference Asadi, Asgharzadeh and Borazjani2018).
3.3. Cascade mechanism for droplet formation: hole formation genesis
In the previous section we discussed the vortex dynamics of the transient jet. This section focuses on the effect of vortices on the formation of inner/outer lobes, the thinning of which leads to genesis of droplets. Figure 9 shows a close view of the interfacial jet dynamics where outer lobes are present. The formation of lobes is observed along the spanwise direction of the jet. Close inspection of the interface shows it has assumed the shape of the surrounding hairpin vortical structures, where a change of the vorticity direction is observed for subsequent lobes (similar observations have been made by Jarrahbashi et al. Reference Jarrahbashi, Sirignano, Popov and Hussain2016; Zandian et al. Reference Zandian, Sirignano and Hussain2018). Consequently, the alignment of the hairpin vortices induces the thinning of the outer lobes, leading to the formation of holes, which expand and eventually give rise to injected water droplets.
Three-dimensional representation of the surface of the jet for case 4 (
$Re=6530$ and 
$We=303$) at 
$t=28.97$ showing the spatial formation of outer lobes. The colour represents the streamwise vorticity field, 
$\omega _x$.
Figure 10 shows the temporal stretching of an outer lobe to form a ligament, and eventually droplets. The ligament orientation is linked to the vortical structures, as suggested above. A bulbous tip is formed at the edge of the ligament driven by capillarity. The ligament detaches from the jet core; then it retracts driven by capillarity to form droplets from both ends according to the so-called ‘end-pinching mechanism’ (Notz & Basaran Reference Notz and Basaran2004). When the stagnant phase has enough momentum, the ligament is dragged into the stagnant phase, setting its characteristic radial length (similarly to what has been reported by Lasheras, Villermaux & Hopfinger Reference Lasheras, Villermaux and Hopfinger1998). The ligament is characterised by having a small thickness prior to its detachment from the base. However, the droplets formed after the breakup process have a larger length scale (in agreement with Villermaux Reference Villermaux2007).
Three-dimensional representation of the interface showing the spatio-temporal evolution of the formation of injected water droplets via the stretching of the outer lobe for case 4 (
$Re=6530$ and 
$We=303$). Dimensionless times are shown in each panel.
Figure 11 shows the formation of droplets inside the jet core, and the mushroom structure. Specifically, figure 11(b–k) depicts the formation and elongation of the inner lobe in the streamwise direction. As before, the stretching of the interface gives rise to a thin sheet as a result of the mutual induction of neighbouring hairpin vortices. The hairpin vortex located above the liquid sheet induces flow downwards, whereas the vortex located under the sheet induces flow upwards. Under the joint action of the vortical structures, the liquid sheet thins until it is perforated (see figure 11e). Following the formation of the hole, retraction of the liquid sheet is radially driven by capillarity (see figures 11( f) and 11(g)), and as shown in figure 11(i), it subsequently gives rise to the formation of a curved liquid thread or ligament. Once more, ligament retraction is driven by capillarity where bulbous ends that form initially undergo ‘end pinching’ to form smaller droplets, as shown in figure 11(j,k). The inertia-induced mechanism for the sheet thinning presented here has been previously reported by Jarrahbashi et al. (Reference Jarrahbashi, Sirignano, Popov and Hussain2016).
(a) Three-dimensional representation of the interface for case 4 (
$Re=6530$ and 
$We=303$) at 
$t=7.33$ showing a selected region in the jet core where inner lobes will form, together with two probe locations in the leading-edge structure, where the arrows indicate the direction of the view. (b–k) Cascade mechanism for the formation of entrapped oil droplets within the jet core at 
$t=(8.84, 9.96, 10.10, 10.92, 11.36, 11.41, 12.06, 13.61, 16.38)$, respectively. (l) Illustration of the entrapped oil droplets inside of the mushroom-like structure at 
$t=7.33$ from the back (probe ‘A’) and front (probe ‘B’) of the structure, respectively. A magnified region of the structure is also shown to illustrate the hole formation behind the rim edge.
The analysis will now focus on the dynamics underneath the mushroom-like structure. Figure 11 shows entrapped oil droplets inside of the leading-edge structure. The formation of these droplets is also linked to the interfacial rupture of the toroidal liquid sheet. The liquid sheet retracts driven by surface tension to accumulate liquid in a rim at its edge. The rim experiences a spanwise destabilisation (i.e. Rayleigh–Plateau type), which leads to a non-uniform rim radius. The film retraction is also accompanied by the formation of interfacial capillary waves that precede the rim. The capillary waves vary the film thickness, and consequently, induce the perforation of the film adjacent to the rim in regions where the film thickness is sufficiently small (in agreement with Mirjalili, Chan & Mani Reference Mirjalili, Chan and Mani2018). The radial expansion of multiple holes yields the formation of liquid threads which experience a capillary instability to produce droplets. Following their formation, the entrapped oil droplets rotate inside the leading structure due to their interaction with the vorticity field (not shown here). The rotation leads to complex interfacial phenomena such as coalescence or collision (similar to that reported by Desjardins & Pitsch Reference Desjardins and Pitsch2010). The coalescence has been observed not only between droplets, but also between droplets and ligaments or droplets and the jet core (see figure 12).
Three-dimensional representation of the interface for case 4 (
$Re=6530$ and 
$We=303$) at dimensionless times shown in each panel. A magnified region of the structure is also shown to illustrate droplet-droplet coalescence and the droplet rotation inside of the leading-edge structure.
Next, we aim to explain the validity of the hole formation mechanism as a result of liquid film thinning. Our numerical method does not consider the destabilising effect brought about by the van der Waals attractive forces as the film thickness tends to zero. Churaev, Derjaguin & Muller (Reference Churaev, Derjaguin and Muller1987) reported that while the thinning takes place, the dynamics of the film enters into an asymptotic behaviour leading to its rupture, depending only on a balance between viscous, van der Waals and surface tension forces. After the film puncture, the circular hole expansion is driven by capillarity, and fluid is accumulated in a circular rim which grows in size as it moves away from the puncture point. The retraction speed 
$V_{TC}$ can be estimated by the classical ‘Taylor–Culick’ theory, 
$V_{TC}= (2 \sigma /\rho _o h)^{1/2}$, where 
$h$ is the film thickness. Prior to the puncture 
$h \sim 26\ \mathrm {\mu }$m, giving an estimate of the retraction velocity as 
$V_{TC} \sim 1.8$ m s
$^{-1}$. The measured retraction speed of the holes in our simulations ranges from 
$0.75$ to 
$1.8$ m s
$^{-1}$. This ensures that although our simulations cannot predict the exact location of the puncture, and their formation is mesh dependent, its expansion is well predicted, ensuring that the physics is fully resolved following the formation of the hole. Additionally, the film thickness of the sheet prior to its rupture agrees with the findings of Marston et al. (Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016), Lhuissier & Villermaux (Reference Lhuissier and Villermaux2009) and Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017) who observed the formation of holes in dynamics sheets of the same order of magnitude (e.g. Marston et al. (Reference Marston, Truscott, Speirs, Mansoor and Thoroddsen2016) and Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017) reported film rupture in the range of 
$9 - 16\ \mathrm {\mu }$m and about 
$22\ \mathrm {\mu }$m, respectively).
3.3.1. Droplet size distribution
This section draws attention to the size distribution of the droplets, formed as a result of the ruptures of the liquid threads. Attention is focused on case 4, which is in the inertia-dominated regime characterised by high 
$Re$ and 
$We$; as discussed above, the flow in this case is accompanied by significant drop creation. During the early stages, the formation of droplets is only observed inside the leading-edge structure owing to film rupture via hole formation. At later times, both entrapped oil droplets and injected water droplets coexist in the computational domain. We have identified all the droplets inside the entire domain, and each droplet diameter is calculated through knowledge of its volume, and the assumption of a spherical shape.
Figure 13 shows the probability density function (p.d.f.) of the entrapped oil droplets and injected water droplets normalised by the mean droplet diameter at different temporal stages: 
$t=9.37,\ 20.37$ and 
$28.52$. The shape of the distribution for the entrapped droplets (see figure 13a) remains largely unchanged with increasing time except for the development of a tail; this represents the creation of larger numbers of smaller droplets due to ligament breakup as explained previously. It is possible to obtain information regarding the characteristics of the droplets in the domain; for instance, the average mean diameters of the entrapped and injected droplets are 338 and 
$291\ \mathrm {\mu }$m, respectively. Figure 13(b) shows that the p.d.f. for the injected droplets has a similar shape to that of the entrapped ones depicted in figure 13(a), and the distribution remains essentially unaltered for 
$t=20.37$ and 
$28.52$.
Probability density function for the entrapped oil droplets and injected water droplets for case 4 (
$Re=6530$ and 
$We=303$) at 
$t=9.37, 20.37$ and 
$28.52$ shown in (a,b), respectively. The distributions are normalised by the calculated mean diameter for each time, 
$d_m$.
4. Conclusions
Three-dimensional numerical simulations using a hybrid interface-tracking/level-set approach were carried out for turbulent water jets entering a stagnant and more viscous silicon oil phase. Particular attention has been paid to the temporal interfacial dynamics that arise as a result of the vortex-surface interaction. Using less computationally expensive simulations, we explored the 
$Re$–
$We$ space, identifying four distinct regions based on the appearance of the interfacial structures. On this basis, we carried out high-fidelity simulations representative of each region which showed that the streamwise vorticity plays a major role in the development of the three-dimensional instabilities on the jet surface.
At low 
$Re$ and 
$We$ numbers (i.e. in the capillary-controlled regime), the jet interface remains axisymmetric and no surface corrugations are observed due to strong capillarity. Moreover, this region is also characterised by the fact that the streamwise vorticity, 
$\omega _x$, never becomes comparable to the azimuthal vorticity, 
$\omega _{\theta }$: the azimuthal vorticity always remains two orders of magnitude larger than the streamwise vorticity. At low 
$Re$ and high 
$We$ numbers, the reduction of capillarity enhances the formation of a mushroom-like structure at the leading edge of the jet, and the formation of corrugations in the jet core; however, 
$\omega _{x}$ is still almost two orders of magnitude smaller than 
$\omega _{\theta }$. At high 
$Re$ and low 
$We$ numbers, streamwise vorticity gains in importance resulting in further deformation of the jet core.
In the inertia-controlled regime (i.e. for high 
$Re$ and 
$We$ numbers), the streamwise vorticity becomes comparable to the azimuthal vorticity triggering the streamwise deformation of KH vortices. We show the formation of hairpin vortices near the interface, which are aligned in the streamwise direction forming layers of inner and outer hairpin vortices. The formation of inner and outer lobes in the jet core are closely linked to their neighbouring hairpin vortices. The alignment of vortices enhances the perforation of the lobes to form holes which expand radially by capillarity, and ultimately give rise to the formation of droplets. Another main feature of the inertia-controlled regime is the creation of a thin toroidal sheet inside of the leading-edge structure. The thinning of this sheet leads to its perforation to form droplets. The droplets rotate as a result of their interaction with the vorticity field, and further topological transitions occur (e.g. coalescence).
While in the present study we have focused on the spatial and temporal development of the interfacial dynamics as a result of the vortex-surface interaction, further investigations should be carried out on the statistical response of the turbulent multiphase flow, similar to the work presented by Ling et al. (Reference Ling, Fuster, Tryggvason and Zaleski2019) for a mixing layer between parallel gas and liquid streams. Additionally, our study assumes a constant value of surface tension, but it is known that streams are usually contaminated with deliberately placed or naturally occurring surfactants, which reduce the surface tension and give rise to surface tension gradients and Marangoni stresses. As shown recently by Constante-Amores et al. (Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020), surfactants are capable of inhibiting capillary singularities and rigidifying the interface, subsequently changing the fate of the atomisation. Hence, the presence of surfactants may change the present results, providing us with an exciting avenue of future research.
Acknowledgements
We also acknowledge HPC facilities provided by the Research Computing Service (RCS) of Imperial College London for the computing time. The numerical simulations were performed with code BLUE (Shin et al. Reference Shin, Chergui and Juric2017) and the visualisations have been generated using ParaView.
Funding
This work is supported by the Engineering and Physical Sciences Research Council, United Kingdom, through a studentship for C.R C-A. in the Centre for Doctoral Training on Theory and Simulation of Materials at Imperial College London funded by the EPSRC (grant no. EP/L015579/1), and through the EPSRC MEMPHIS (grant no. EP/K003976/1) and PREMIERE (EP/T000414/1) Programme Grants. The authors would like to acknowledge the funding and technical support from BP through the BP International Centre for Advanced Materials (BP-ICAM), which made this research possible. O.K.M. also acknowledges funding from PETRONAS and the Royal Academy of Engineering for a Research Chair in Multiphase Fluid Dynamics, and the PETRONAS Centre for Engineering of Multiphase Systems. D.J. and J.C. acknowledge support through 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 no. 2020A0082B06721.
Declaration of interests
The authors report no conflict of interest.
Appendix A
A.1. Mesh study and resolution considerations
Solving the small scales of the atomisation phenomenon is a challenging process, and in order to ensure the physical validity of the results, we have assessed the grid dependence nature of our results by performing a mesh study for case 4 (see table 5). The appropriate choice for the minimum grid length scale will depend on either the vanishing interfacial singularity or the smallest turbulence length scale (i.e. Kolmogorov length scale).
For single-phase jets, Pope (Reference Pope2000) suggested that the Kolmogorov length scale 
$\eta$ is resolved if 
${\rm \Delta} x / \eta \leq 2.1$, where 
${\rm \Delta} x$ represents the minimum computational cell size, and 
$\eta$ is estimated by 
$\eta = Re^{-3/4} l$, where 
$l$ is the length scale. Previous computational studies of temporal and spatial multiphase jets are presented in table 4, showing the current state-of-the-art in terms of Pope's criterion. It is evident that Pope's criterion is difficult to meet due to the high computational costs of the simulations. In their recent work, Ling et al. (Reference Ling, Fuster, Zaleski and Tryggvason2017, Reference Ling, Fuster, Tryggvason and Zaleski2019) performed simulations of a two-phase mixing layer between parallel gas and liquid streams to investigate the interfacial dynamics and the statistics of the multiphase turbulence, estimating the Kolmogorov length scale to be 
$\eta \sim 0.945\ \mathrm {\mu }$m, which leads to 
${\rm \Delta} x / \eta \sim 3$. However, through their numerical simulations, they estimated that 
$\eta$ is larger in size (e.g. 
$\eta = 3 - 4.5\ \mathrm {\mu }$m), and they also showed that their lower resolution mesh, i.e. 
${\rm \Delta} x / \eta \sim 6$ (see their figure 19d), was capable of predicting similar results in terms of turbulence dissipation.
List of computational studies of atomisation showing compliance with the Pope criterion (Pope Reference Pope2000).
Characteristics of different mesh sizes used to study the jet dynamics in this study.
For our highest Reynolds number, our simulation does not meet Pope's criterion. But as shown by Ling et al. (Reference Ling, Fuster, Tryggvason and Zaleski2019), the actual 
$\eta$ could be larger. Additionally, the atomisation of the injected phase in a stagnant viscous phase alleviates the range of relevant physical scales.
The second biggest challenge of computational atomisation is ‘numerical breakup’ of liquid threads. A coarse grid would trigger the formation of thicker numerical threads, and consequently, the formation of larger droplets (this problem has been previously reported by Shinjo & Umemura (Reference Shinjo and Umemura2010), Herrmann (Reference Herrmann2011), Gorokhovski & Herrmann (Reference Gorokhovski and Herrmann2008) and Jarrahbashi & Sirignano (Reference Jarrahbashi and Sirignano2014)). Shinjo & Umemura (Reference Shinjo and Umemura2010) stated that the mesh should be refined up to the point where the dynamics of the thread are solely governed by surface tension forces. These forces would trigger the formation of capillary waves during the retraction of these ligaments, giving rise to the onset of the Rayleigh–Plateau instability (i.e. the ‘end-pinching’ mechanism). After those mechanisms are initialised, the thread dynamics enter in an asymptotic behaviour towards the interfacial singularity (i.e. a refined mesh would not affect the size of the resulting droplets). On this basis, the numerical resolution regarding the interfacial length scales will be assessed following the methodology proposed by Ménard et al. (Reference Ménard, Tanguy and Berlemont2007) and Desjardins & Pitsch (Reference Desjardins and Pitsch2010), who used a ‘grid-based Weber number’, defined by 
$We_{{\rm \Delta} x_{min}} = \rho _{w} {U^{2}}{\rm \Delta} x_{min} / \sigma$. This equation provides us with the smallest interfacial length scale that the simulation is capable of resolving by assuming that the smallest interfacial structure is equal to the minimum mesh size of the computational domain. Ménard et al. (Reference Ménard, Tanguy and Berlemont2007) suggested that no further breakup is observed for values under 
$10$. For a Reynolds number corresponding to 
$Re \sim 10^{4}$ and the M3 type mesh, the grid-based Weber number is 
$We_{{\rm \Delta} x_{min}} \sim 4.12$, which meets the above criterion, suggesting that all capillary singularities would be resolved.
Therefore, we have proved that the M3 mesh is capable of resolving turbulent scales and interfacial singularities, and consequently, detailed analysis of interfacial and vortical structures is performed using a M3 mesh type (unless stated otherwise). Figure 14 shows the temporal evolution of the kinetic energy, 
$E_k=\int _V (\rho \boldsymbol {u}^{2})/2 \,\rm {d} V$, the interfacial area and the maximum axial location of the jet tip (i.e. leading edge) for different meshes. Additionally, the vorticity profiles were checked between the M2 and M3 mesh, and no significant differences were found.
Mesh study for case 4 (
$Re=6530$ and 
$We=303$). The panels highlight the temporal evolution of the kinetic energy 
$E_k$ (a), the interfacial area (b) and the maximum axial location of the jet tip (c).
A.2. Linear stability analysis
Following Plateau (Reference Plateau1873) and Rayleigh (Reference Rayleigh1879), we will show that our numerical method is capable of predicting the growth rates in relation to the capillary breakup of jets (i.e. Rayleigh–Plateau instability). On this basis, we have considered an inviscid cylindrical jet of radius 
$R_o$ and density 
$\rho$, surrounded by a dynamically passive gas. We consider the evolution of infinitesimal perturbations to linearise the equations. Thus, the perturbed jet is represented by
\begin{equation} R(\omega,\varepsilon)=R_o \left (1+\varepsilon \cos\left ( \omega x\right )\right )= R_o \left (1+\varepsilon \cos\left ( \frac{2{\rm \pi} }{\lambda } x\right )\right ), \end{equation}
where 
$\varepsilon$, 
$\lambda$, 
$\omega$ stand for the amplitude of infinitesimal perturbation, wavelength and the growth rate, respectively. After the linearisation of the governing equations and retaining terms only to order of 
$\varepsilon$ results in a dispersion relation, which indicates the dependence of the growth rate on the wavenumber (i.e. 
$k=2 {\rm \pi}/\lambda$) of the inviscid jet, expressed as
\begin{equation} \omega^{2}=\frac{\sigma \left(kR_o\right) }{\rho R_o^{3}}\frac{I_1(kR_0)}{I_0(kR_0)} (1-k^{2}R_o^{2})=\omega_o^{2} x \frac{I_1(x)}{I_0(x)} (1-x^{2}), \end{equation}
where 
$x=kR_o$ stands for a reduced wavenumber and 
$\omega _o^{2}=\sigma /(\rho R_o^{3})$. When 
$x>1$, the initial perturbation presents stable oscillatory solutions; whereas in the other case (i.e. 
$x<1$), the initial perturbation grows exponentially with a growth rate given by (A2).
Figure 15 shows the growth rates obtained through the numerical simulations in comparison to the theoretical dispersion relation. The growth rates from DNS shows a discrepancy under 
$2\,\%$ in all cases in comparison to the theoretical values. Based on this, we can conclude the capability of our numerical method for capturing accurately the dynamics of capillary jets.
Comparison of non-dimensional growth rates between numerical simulations (black markers) and linear theory (solid line).
A.3. Scalings laws for the capillary breakup of liquid threads
As we have presented in the introduction, DNS must be capable of predicting the developing two-phase fluid interfacial dynamics featuring interface breakup, and droplet coalescence. In light of this, we aim to show the capabilities of our numerical framework in the prediction of the scaling laws for the capillary singularity of liquid threads. As the point of singularity approaches, the system is driven locally by the large interfacial curvature, and its interfacial dynamics depends solely upon the physical properties of the liquid. Lister & Stone (Reference Lister and Stone1998) have suggested that the pinchoff of a viscous thread (of radius 
$r(z,t)$, density 
$\rho$, viscosity 
$\mu$ and surface tension 
$\sigma$) surrounded by another viscous fluid transits between different dynamical regimes as 
$r(t) \rightarrow 0$ towards the breakup time 
$\tau$ (see figure 1 of Lister & Stone Reference Lister and Stone1998). The thinning of a water thread, surrounded by air, transitions from an inertial-capillary regime (
$r\sim \tau ^{2/3}$ and 
$u\sim \tau ^{-1/3}$) to an inertial-viscous regime (
$r\sim \tau$ and 
$u\sim \tau ^{-1/2}$), more details can be found in Day, Hinch & Lister (Reference Day, Hinch and Lister1998), Eggers (Reference Eggers1993) and Lister & Stone (Reference Lister and Stone1998).
Therefore, we consider the capillary singularity of a water thread (of initial radius 
$R$) surrounded by air in absence of gravity as the interfacial dynamics is locally driven by a balance of viscous capillary and inertial forces, which can be expressed by the Ohnesorge number 
$Oh=\mu /(\rho \sigma R)$. As 
$Oh <1$, at early thinning stages there is a competition between the fluid inertia and the opposing capillary pressure, and the predicted thread radius towards the singularity agrees with the theoretical scalings (see figure 16) as the thread transits between different regimes (in agreement with Castrejón-Pita et al. Reference Castrejón-Pita, Castrejón-Pita, Thete, Sambath, Hutchings, Hinch, Lister and Basaran2015). After the singularity point, the formation of a satellite ligament is observed which has an initially cylindrical shape (with 
$Oh=2.62 \times 10^{-3}$ and 
$L_o=11.55$); however, it undergoes retraction driven by capillarity to form more spherical droplets. This process known as ‘end pinching’ has been previously well described by Notz & Basaran (Reference Notz and Basaran2004). By performing this analysis we have proved the capability of our numerical technique to predict the dynamics of the capillary singularity of a liquid thread and its post-breakup events (table 6 shows a summary of the different meshes evaluated for this study).
Scaling laws for the capillary singularity of a water thread surrounded by air for different numerical resolutions (see table 6). Here, the Ohnesorge number is 
$1.178 \times 10^{-3}$. The minimum thread radius, (a), and the maximum streamwise velocity, (b), vs the remaining time to breakup 
$\tau$ agree with the inertial-capillary and inertial-viscous regimes presented by Eggers (Reference Eggers1993) and Day et al. (Reference Day, Hinch and Lister1998). Additionally, a three-dimensional representation of the interface is shown to depict the retraction of the satellite ligament to produce daughter droplets via the ‘end-pinching’ mechanism (in agreement with Notz & Basaran Reference Notz and Basaran2004).
Characteristics of different mesh sizes used to study the capillary breakup of a water thread.
Additionally, the accuracy and validation of the numerical method has been previously addressed to other complex interfacial phenomena. These phenomena include breakup and recoiling of liquid threads (Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2020), falling film flows (Batchvarov et al. Reference Batchvarov, Kahouadji, Constante-Amores, Gonçalves, Shin, Chergui, Juric and Matar2020a), propagation of elongated bubbles in channels (Batchvarov et al. Reference Batchvarov, Kahouadji, Magnini, Constante-Amores, Craster, Shin, Chergui, Juric and Matar2020b), bubbles undergoing bursting (Constante-Amores et al. Reference Constante-Amores, Kahouadji, Batchvarov, Shin, Chergui, Juric and Matar2021) and drops coalescing partially or completely with deformable interfaces.
A.4. Time step
Finally, the temporal integration scheme is based on a second-order Gear method, with implicit solution of the viscous terms of the velocity components. The time step 
${\rm \Delta} t$ is adaptive to ensure stability, and it is defined through the criterion
\begin{equation} {\rm \Delta} t = \min \left\{{\rm \Delta} t_{cap},{\rm \Delta} t_{vis}, {\rm \Delta} t_{CFL},{\rm \Delta} t_{int}\right\}, \end{equation}
where 
${\rm \Delta} t_{cap}$, 
${\rm \Delta} t_{vis}$, 
${\rm \Delta} t_{CFL}$, 
${\rm \Delta} t_{int}$ stand for the capillary time step, the viscous time step, the Courant–Friedrichs–Lewy (CFL) time step and interfacial CFL time step, respectively. These terms are defined by
\begin{equation} \left.\begin{gathered} {\rm \Delta} t_{vis} = \min \left ( \frac{\rho_{w}}{\mu_{{w}}}, \frac{\rho_{{so}}}{\mu_{{so}}}\right ) \frac{ {\rm \Delta} x^{2}_{\min} }{6 },\quad {\rm \Delta} t_{cap} =\frac{1}{2} \left ( \frac{(\rho_{w} +\rho_{{so}}) {\rm \Delta} {x_{\min}}^{3}}{{\rm \pi} \sigma }\right )^{1/2},\\ {\rm \Delta} t_{CFL} = \underset{j}{\min} \left ( \underset{domain}{\min} \left (\frac{ {\rm \Delta} x_{j} }{u_j} \right )\right ),\quad {\rm \Delta} t_{int} = \underset{j}{\min} \left ( \underset{\varGamma (t)}{\min} \left (\frac{ {\rm \Delta} x_{j} }{\left \| \boldsymbol{V} \right \|}\right )\right ), \end{gathered}\right\} \end{equation}
where 
${\rm \Delta} {x_{\min }} = \min _j ({\rm \Delta} x_j)$. In our simulations the adaptive time step is controlled by 
${\rm \Delta} t_{int}$ which is 
$O(10^{-5})$s.
