# Averaged equations for disperse two-phase flow with interfacial properties and their closures for dilute suspension of droplets

Nicolas Fintzi, Jean-Lou Pierson

► To cite this version:

Nicolas Fintzi, Jean-Lou Pierson. Averaged equations for disperse two-phase flow with interfacial properties and their closures for dilute suspension of droplets. International Journal of Multiphase Flow, 2026, 194, pp.105424. 10.1016/j.ijmultiphaseflow.2025.105424. hal-05491443

HAL Id: hal-05491443

https://ifp.hal.science/hal-05491443v1

Submitted on 3 Feb 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_2.png" /></figure>

Research Paper

# Averaged equations for disperse two-phase flow with interfacial properties and their closures for dilute suspension of droplets

Nicolas Fintzi<sup>a,b,□</sup>,*, Jean-Lou Pierson<sup>a,□</sup>

<figure><img src="image_3.png" /></figure>

## ARTICLE INFO

**Keywords:**

Dispersed Two-phase flow

Averaged equations

Lagrangian formulation

Interfacial transport

Hybrid modeling

Moment equations

## ABSTRACT

This article provides a derivation of the averaged equations governing the motion of dispersed two-phase flows with interfacial transport. We begin by revisiting the two-fluid formulation, as well as the distributional form of the interfacial transport equation which holds on the entire domain. Following this, a general Lagrangian model is introduced, which accounts for the effects of both internal and interfacial properties of the dispersed inclusions (bubbles, droplets, or particles) within a continuous phase. This is achieved by derivation of conservation laws for particle surface and volume-integrated properties. By summing the internal and interfacial conservation laws, we derive a conservation equation for an arbitrary Lagrangian property associated with the inclusion. We then proceed by deriving the lesser-known conservation equations for the moments of the volume and surface distribution of an arbitrary Lagrangian property. Next, the averaged equations for the dispersed phase are derived through two distinct approaches: the particle-averaged (or Lagrangian-based) formalism, and the phase-averaged method. One important conclusion of this work is the demonstration of the relationship between the particle-averaged and phase-averaged equations. We show that the dispersed phase-averaged equations can be interpreted as a series expansion of the particle-averaged moment equations. We then present a "hybrid" set of equations, consisting of phase-averaged equations for the continuous fluid phase, complemented by an arbitrary number of moment conservation equations for the dispersed phase. To further illustrate the methodology, we derive the mass, momentum, second moment of mass and first moment of momentum equations for droplets or bubbles suspended in a Newtonian fluid. In particular, we highlight the role of the second-order moment of mass equation and first-order moment of momentum equation, which link droplet deformation to the stresslet. We then derive closure laws in the dilute, viscous-dominated regime, with particular emphasis on the effects of surface tension gradients. Additionally, we discuss several covariance closure terms that emerge in the averaged equations. Finally we demonstrate how the leading order deformation of the droplets can be obtained thanks to the second-order mass moment and first moment of momentum equation.

## 1. Introduction

Dispersed multiphase flows are encountered across a wide range of chemical engineering applications. These include gas-solid interactions in fluidized bed reactors, liquid-liquid flows in extractors, and gas-liquid bubbly flows in flotation processes. These systems exhibit a wide range of scales, from the size of individual inclusions (as small as a few micrometers) to the size of the reactor (often exceeding one meter), making fully resolved simulations computationally impractical. As a result, the current engineering practice relies on averaged equations of motion for both the dispersed and continuous phases. Therefore, developing a robust set of averaged equations that accurately captures the complexities of dispersed multiphase flows is essential. In this study, we aim to overcome the limitations of existing models by proposing a comprehensive set of averaged equations that can be applied to various types of dispersed multiphase flows, including interfacial transport processes. Indeed, in many chemical

<footnote><sup>a</sup> IFP Energies Nouvelles, Rond-point de l'échangeur de Solaize, 69360 Solaize, France</footnote>

<footnote><sup>b</sup> Sorbonne Université, Institut Jean le Rond d'Alembert, 4 place Jussieu, 75252 PARIS CEDEX 05, France</footnote>

<footnote>* Corresponding author at: IFP Energies Nouvelles, Rond-point de l'échangeur de Solaize, 69360 Solaize, France.</footnote>

<footnote>E-mail address: nicolas.fintzi@ifpen.fr (N. Fintzi).</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 1. Topology of dispersed two-phase flows.</figcaption></figure>

engineering processes involving liquid or gas inclusions, mass transfer occurs between phases. This interphase mass transfer can generate gradients in surface tension, which in turn give rise to various phenomena such as Marangoni convection (Wegener et al., 2014) or capillarity-induced motion of a droplet due to concentration gradient (Subramanian, 1985).

The majority of averaged two-fluid models are based on the framework proposed by Drew (1983), where both the dispersed and continuous phases are treated similarly (see, for example, Hu et al. (2021)). Despite its versatility-allowing applications beyond dispersed flows, such as stratified or slug flow, the physical interpretation of the closure terms remains difficult (Drew, 1983), and the mathematical well-posedness of the entire system is still a matter of debate (Panicker et al., 2018; Lhuillier et al., 2013). An alternative approach treats the dispersed phase using a Lagrangian framework. This approach originates from the kinetic theory of gases, where the motion of particles is described by a single-particle distribution function governed by the Boltzmann equation (Chapman and Cowling, 1990). It has proven particularly successful in predicting the dynamics of dry granular materials (Rao et al., 2008) and gas-solid flows (Simonin, 1996). However, the fluid phase and the dispersed phase are handled through two distinct formalisms: phase averaging for the fluid and Lagrangian averaging for the dispersed phase. In a pioneering study, Buyevich and Shchelchkova (1979) demonstrated that these two averaging methods could be connected by performing a Taylor expansion of the closure terms around the particle's center of mass. This "hybrid" formalism was later revisited by Lhuillier (1992) and by Zhang and Prosperetti (1994a,b). The latter provided closures for the momentum transport equations in the inviscid flow limit. Subsequently, Jackson (1997) and Zhang and Prosperetti (1997) provided appropriate closures for the Stokes flow regime for dilute suspensions involving spherical solid particles and spherical droplets, respectively.

Many previous studies that used the “hybrid” formalism have focused on solid particles (Buyevich and Shchelchkova, 1979; Jackson, 1997). Notable exceptions include Zhang and Prosperetti (1994b), who investigated spherical bubbles with time-varying radii, and Zhang and Prosperetti (1997), who examined spherical droplets. Despite these advancements, the question remains: how can different inclusion shapes, internal fluid motion, and surface transport equations be incorporated into such Lagrangian models? To the best of our knowledge, important surface properties, such as variable or constant surface tension, the arbitrary shape of fluid particles, and the distribution of surfactants, have not yet been fully incorporated into current averaged hybrid models. While the Lagrangian formulation of the interfacial area equation is well established (Lhuillier et al., 2000), its application to define general Lagrangian quantities for a dispersed phase is still not fully explored.

Several authors have addressed the question of equivalence between Lagrangian (or particle-averaged) and phase-averaged equations in various contexts (Zhang and Prosperetti, 1997; Lhuillier et al., 2000; Nott et al., 2011). In Zhang and Prosperetti (1997, Appendix A), the authors demonstrated that these two frameworks are equivalent for spherical inclusions when only considering the first-order moments. Later, Nott et al. (2011, Appendix A) extended the proof of Zhang and Prosperetti (1997), showing that the Lagrangian and phase-averaged equations are strictly equivalent for suspensions of solid spherical particles by considering an infinite number of higher-order terms. In addition, Lhuillier et al. (2010) argued that the phase-averaged equations applied to the dispersed phase are, in fact, a Taylor series expansion of the particle-averaged moment equations. Similarly, in the context of interfacial area balance equations, Lhuillier et al. (2000) reached a comparable conclusion when comparing the phase-averaged and particle-averaged area density equations for spherical particles. These studies focus on monodisperse suspensions of spherical particles and all arrive at the same conclusion: the particle-averaged equation is rigorously equivalent to the phase-averaged equation for the dispersed phase.

Following the above review, there is a need for a comprehensive hybrid model that encompasses all relevant physical aspects, including surface properties and arbitrarily shaped particles. Although this problem has been studied by several research groups including Lhuillier (1992) and Zhang and Prosperetti (1997); our contribution lies in presenting a unified and simple theoretical framework for deriving the hybrid set of governing equations. Here the Lagrangian balance equations are derived without any simplifying assumptions, thus allowing for: (1) the incorporation of interfacial properties while retaining a Lagrangian framework, and (2) the representation of internal gradients within droplets through the use of distribution moments. Additionally we provide a general relationship between the phase-averaged and particle-averaged frameworks for inclusions of arbitrary shapes, including interfacial transport phenomena. We rigorously demonstrate the relation between the Lagrangian-averaged and Eulerian-averaged formulations, thereby extending earlier results by Nott et al. (2011). A central feature of our approach is the emphasis on moment equations for the dispersed phase, as previously evidenced by Lhuillier and Nadim (2009) in the context of solid particles. We argue that, regardless of the specific problem under consideration — including those involving fluid particles of complex shape — the hybrid formulation offers a more physically grounded alternative to the traditional two-fluid model for describing dispersed flows. To illustrate our methodology, we derive the mass, momentum, and two first-order moment equations for droplets or bubbles suspended in a Newtonian fluid. In particular, we present the first-order moment of momentum equation, which connects droplet deformation to the stresslet, extending the classical result of Batchelor (1970).

Despite the extensive work on dispersed two-phase flows, only a few theoretical results exist that fully account for the whole closure problem. In the regime of Stokes flow and for very dilute suspensions, Jackson (1997) and Zhang and Prosperetti (1997) independently derived the momentum closure for spherical solid particles, with Zhang and Prosperetti (1997) extending the work to spherical drops. The inclusion of the symmetric part of the first moment of hydrodynamic forces, namely the stresslet, leads to Einstein's viscosity formula (Lhuillier, 1996). Even in this regime, the averaged fluid phase exhibits non-Newtonian behavior, induced by the relative motion of the particles due to the incorporation
of the second moment of hydrodynamic forces (Lhuillier, 1996; Zhang and Prosperetti, 1997). In this work we provide closure laws in the dilute, viscous-dominated regime, with special attention to the role of surface tension gradients and the closure of higher-order moment equations.

In Section 2, we begin by presenting the two-fluid formulation and the distributional form of the interfacial transport equation. Next, in Section 3, we introduce a Lagrangian-based model that describes individual fluid particles with arbitrary shapes and surface properties. In Section 4, we highlight the link between the particle-averaged equations and the phase-averaged equation for any type of dispersed inclusions. Building on this demonstration, we introduce a hybrid set of conservation equations: one equation for the fluid phase and multiple equations for the dispersed phase, representing the conservation of moments of the distribution of the conserved quantity. We next derive the conservation equations for mass, momentum, and higher-order moments for both the dispersed and fluid phases, as presented in Section 5. Subsequently, in Section 6, we focus on deriving the closure relations in the limit of dilute, viscosity-dominated flows. Finally, we discuss the findings of this investigation in Section 7.

## 2. The two-fluid formulation

In this section, we derive the conservation equations using a two-fluid formulation. While this approach has been employed in numerous studies, including those by Kataoka (1986), Lhuillier et al. (2010), Ishii and Hibiki (2010), Morel (2015), and Bothe (2022), our method enables us to establish general results that encompass previous work.

We consider a system consisting of two phases, separated by a sharp interface $Γ(t, \mathcal{F})$, where $t$ represents the current time and $\mathcal{F}$ a flow realization over a set of possible dispersed multiphase flow realization. A realization can be represented by a list of variables such as particle positions, velocities and other factors that provide a complete and unique description of the flow (Zhang, 2021). The phase subdomains are labeled as $Ω_f(t, \mathcal{F})$ and $Ω_d(t, \mathcal{F})$, corresponding to the continuous (fluid) phase ($f$) and the dispersed phase ($d$), respectively (see Fig. 1). The complete domain, $Ω$, is composed of the union of $Ω_f$, $Ω_d$, and $Γ$. To track the position of phase $k$ and the interfaces, we introduce the phase indicator function $χ_k$ defined as follows

$$χ_k(\mathbf{x}, t, \mathcal{F}) = \begin{cases} 1 & \text{if } \mathbf{x} \in Ω_k(t, \mathcal{F}) \\ 0 & \text{if } \mathbf{x} \notin Ω_k(t, \mathcal{F}) \end{cases} \quad \text{for } k = f, d. \qquad (2.1)$$

Additionally, we define $Ω_α(t, \mathcal{F})$ as the region occupied by particle $α$ at time $t$ within the configuration $\mathcal{F}$ (Fig. 1). The collective union of $Ω_α(t, \mathcal{F})$ for $α = 1, ..., N$, where $N$ denotes the number of particles present in the flow, corresponds to $Ω_d(t, \mathcal{F})$. Similarly, we assume that $Γ(t, \mathcal{F})$ can be partitioned into $N$ subregions, labeled $Γ_α(t, \mathcal{F})$, each representing the surface of particle $α$.

### 2.1. Topological equations

Using the distribution formalism, one may show that $χ_k(\mathbf{x}, t, \mathcal{F})$ obeys the following relations (Drew, 1983).

$$\partial_t χ_k + \mathbf{u}_Γ^0 \cdot \nabla χ_k = 0, \qquad (2.2)$$

$$\nabla χ_k = -δ_Γ n_k, \qquad (2.3)$$

where $\mathbf{u}_Γ^0(\mathbf{x}, t, \mathcal{F})$ is the velocity of the interface and $δ_Γ(\mathbf{x}, t, \mathcal{F})$ denotes the Dirac function localized on the interface, also called the interfacial indicator function (Drew, 1983; Junqua-Moullet, 2003). More specifically the distribution $δ_Γ$ is defined as (Appel, 2007)

$$\langle δ_Γ, ϕ \rangle = \int_Γ ϕ dΓ \qquad (2.4)$$

where $ϕ$ is a function with compact support. In this work, we use the subscript $Γ$ to denote any quantity inherently defined on the interfaces, such as $\mathbf{u}_Γ^0(\mathbf{x}, t, \mathcal{F})$ and $δ_Γ(\mathbf{x}, t, \mathcal{F})$. Additionally, we employ the superscript 0 on a field to indicate that it is defined at the local, non-averaged level, meaning it is a function of the flow configuration $\mathcal{F}$, the position vector $\mathbf{x}$, and time $t$. For clarity and readability, we omit the arguments $\mathcal{F}$, $\mathbf{x}$, and $t$ for any local function denoted by 0 (as well as for $δ_Γ$ and $χ_k$), as their dependence on $\mathbf{x}$, $t$, and $\mathcal{F}$ is implied.

To derive a conservation equation for $δ_Γ$, one can take the gradient of Eq. (2.2) and then compute the dot product of the resulting expression with $n_k$. This approach yields the following result (Marle, 1982; Drew, 1990; Lhuillier et al., 2000; Junqua-Moullet, 2003)

$$\partial_t δ_Γ + \nabla \cdot [(\mathbf{u}_Γ^0 \cdot \mathbf{n}) n] δ_Γ = (\mathbf{u}_Γ^0 \cdot \mathbf{n})(\nabla \cdot \mathbf{n}) δ_Γ, \qquad (2.5)$$

where we make use of the relation $n_k \cdot \partial_t n_k = \frac{1}{2} \partial_t (n_k \cdot n_k) = 0$. We did not specify the index of the normal $\mathbf{n}$ in Eq. (2.5) to emphasize that this equation is independent of the index $k$. This is because $\mathbf{n}$ appears twice in each term of this equation. Likewise, we will omit the index of the normal vector in subsequent discussions unless it is required. Additionally, note that in Eq. (2.5) only the normal component of the surface velocity is present, and the term $\nabla \cdot \mathbf{n}$ represents twice the mean curvature of the surface (Aris, 2012). A second equation is obtained by taking the gradient of the distribution $δ_Γ$ (Appendix A), namely,

$$\nabla δ_Γ = \nabla (n δ_Γ) \cdot \mathbf{n}. \qquad (2.6)$$

Eq. (2.6) has been derived straightforwardly, yet it differs from the expression recently obtained by Orlando et al. (2023). The source of this discrepancy is unclear to the current authors. It might be due to differences in their notation, with their operator gradient potentially representing a normal derivative. Eqs. (2.2), (2.3), (2.5) and (2.6) are commonly referred to as the topological equations of two-phase flows. These equations describe the spatiotemporal evolution of the interfaces topology.

### 2.2. Local conservation equations

We now introduce the local conservation laws that govern the fluid inside bulk phases (inside $Ω_d$ and $Ω_f$) and on the interfaces (on $Γ$). Let $f_k^0$ denote a volumetric quantity of arbitrary tensorial order defined in $Ω_k$. Similarly, let $f_Γ^0$ represent an arbitrary surface property defined on $Γ$. Note that $f_Γ^0$ can be defined as a volumetric quantity integrated over the thickness of the interfaces Ishii and Hibiki (2010, Chapter 2). Therefore, if $f_k^0$ has units $[x]$, then $f_Γ^0$ has the unit $[x \cdot L]$, where $L$ is the unit of a length. The subscript $‖$ is used to denote the projection of a field onto the plane
tangential to the surface Γ. Specifically, for an arbitrary quantity g defined on Γ, its tangential projection is given by $g_{||} = \delta_{||} \cdot g$ where $\delta_{||} = (\delta - \mathbf{n}\mathbf{n})$ is a tangential projection operator and $\delta$ is the identity tensor. This definition also extends to the gradient operator, where $\nabla_{||} = \delta_{||} \cdot \nabla$ is referred to as the surface gradient operator. In addition we may define the surface divergence operator as $\nabla_{||} \cdot () = (\delta_{||} \cdot \nabla) \cdot ()$.

Following the strategy outlined in R. Gatignol and Prud'homme (2001), Ishii and Hibiki (2010) and Bothe (2022), the local conservation equations for $f_k^0$ and $f_r^0$ can be written as follows,

$$ \partial_t f_k^0 + \nabla \cdot (f_k^0 \mathbf{u}_k^0 - \mathbf{\Phi}_k^0) = s_k^0 \qquad \text{in } \Omega_k, \qquad (2.7) $$

$$ \partial_t f_r^0 + f_r^0 (\mathbf{u}_r^0 \cdot \mathbf{n})(\nabla \cdot \mathbf{n}) + \nabla_{\|} \cdot (f_r^0 \mathbf{u}_{\|}^0 - \mathbf{\Phi}_{\|}^0) = s_{\|}^0 - \llbracket f_k (\mathbf{u}_{\|}^0 - \mathbf{u}_k^0) + \mathbf{\Phi}_k^0 \rrbracket \cdot \mathbf{n}_k \quad \text{on } \Gamma. \qquad (2.8) $$

The tensors $\mathbf{\Phi}_k^0$ and $\mathbf{\Phi}_{\|}^0$ represent the non-convective fluxes corresponding to the quantities $f_k^0$ and $f_r^0$, respectively. Similarly, $s_k^0$ and $s_{\|}^0$ are the source terms corresponding to the quantities $f_k^0$ and $f_r^0$, respectively. We assume that only the tangential component of the tensor $\mathbf{\Phi}_k^0$ plays a role in the surface balance equation (Bothe, 2022). Consequently, the flux terms of the form $\nabla_{\|} \cdot [(\dots)_{\|}]$ in Eq. (2.8) represent two-dimensional convective and non-convective fluxes arising from tangential motions or diffusive processes at the interface. Conversely, the term involving the local curvature $(\nabla \cdot \mathbf{n})$ is related to the interface expansion or contraction resulting from the normal velocity of the interface $\mathbf{u}_r^0 \cdot \mathbf{n}$. For practical uses, note that the advecting term in Eq. (2.8) can be written in the more compact form $\nabla_{\|} \cdot (f_r^0 \mathbf{u}_r^0)$. Indeed, one can show that $f_r^0 (\mathbf{u}_r^0 \cdot \mathbf{n})(\nabla \cdot \mathbf{n}) + \nabla_{\|} \cdot (f_r^0 \mathbf{u}_{\|}^0) = \nabla_{\|} \cdot (f_r^0 \mathbf{u}_r^0)$, by noticing that $\mathbf{n} \cdot \nabla_{\|} \cdot (\dots) = 0$ and $\nabla_{\|} \cdot \mathbf{n} = \nabla \cdot \mathbf{n}$ (Nadim, 1996). In Eq. (2.8) we have introduced the notation $\llbracket \dots \rrbracket \cdot \mathbf{n}_k$, where $\llbracket \dots \rrbracket \cdot \mathbf{n}_k = \sum_{k=1}^2 \dots \cdot \mathbf{n}_k$ signifies the jump across the interface. Thus, the last term on the right-hand side of Eq. (2.8) represent a source of $f_r^0$ due to the discontinuity of bulk properties on either side of the interface. It is important to note that Eqs. (2.7) and (2.8) are uniquely defined within the domains $\Omega_k$ and $\Gamma$, respectively. Consequently, these equations are referred to as local conservation equations.

## 2.3. The two-fluid formulation

The phase indicator function $\chi_k$, and the Dirac delta function $\delta_r$, allow the extension of Eqs. (2.7) and (2.8) to the entire flow domain $\Omega$. This extension is achieved by employing the methodology introduced by Drew (1983) and Kataoka (1986) for the conserving laws inside the volume (Eq. (2.7)). The two-fluid formulation may be obtained by multiplying Eq. (2.7) by $\chi_k$. Using Eqs. (2.2) and (2.3) we obtain

$$ \partial_t (\chi_k f_k^0) + \nabla \cdot (\chi_k f_k^0 \mathbf{u}_k^0 - \chi_k \mathbf{\Phi}_k^0) = \chi_k s_k^0 + \delta_r \llbracket f_k^0 (\mathbf{u}_r^0 - \mathbf{u}_k^0) + \mathbf{\Phi}_k^0 \rrbracket \cdot \mathbf{n}_k. \qquad (2.9) $$

This yields an equation defined over $\Omega$, to be solved for the quantity $\chi_k f_k^0$ instead of $f_k^0$. Moreover, in contrast to Eq. (2.7) we observe the emergence of the interfacial term $\delta_r \llbracket \dots \rrbracket$ on the right-hand side of Eq. (2.9). To the best of the author's knowledge, the general form of the interfacial transport equation, expressed using the distribution formalism, has not yet been derived. However, some specific forms applicable to a given $f_r^0$ can be found in Marle (1982) and Teigen et al. (2009). The derivation of this equation is straightforward but requires some algebra which is detailed in Appendix B. This yields

$$ \partial_t (\delta_r f_r^0) + \nabla \cdot (\delta_r f_r^0 \mathbf{u}_r^0 - \delta_r \mathbf{\Phi}_{\|}^0) = \delta_r s_{\|}^0 - \delta_r \llbracket f_k^0 (\mathbf{u}_{\|}^0 - \mathbf{u}_k^0) + \mathbf{\Phi}_k^0 \rrbracket \cdot \mathbf{n}_k, \qquad (2.10) $$

which corresponds to the conservation equation of $\delta_r f_r^0$, which is defined over $\Omega$. Although Eqs. (2.8) and (2.10) appear similar in many respects, they possess a fundamental distinction: Eq. (2.8) is defined exclusively on $\Gamma$, it is a two-dimensional partial differential equation with a time-dependent metric, whereas Eq. (2.10) is a three-dimensional partial differential equation. This distinction holds significant importance, as the formulation of Eq. (2.10) proves to be more practical for the numerical modeling of quantities on deformable interfaces (Teigen et al., 2009). The set of equations formed by Eq. (2.9) for $k = f, d$ and the surface transport equation or jump condition (Eq. (2.10)) is commonly known as the two-fluid formulation of multiphase flows (Morel, 2015; Tryggvason et al., 2011; Drew, 1983; Kataoka, 1986).

According to Eq. (2.4) the distribution $\delta_r$ has the dimension of the inverse of a length, $[L^{-1}]$ (Appel, 2007). Consequently, if the dimensionality of $f_k^0$ is assumed to be $[x]$, it follows that the dimensionality of $f_r^0$ must be $[x \cdot L]$. Therefore, the product $\delta_r f_r^0$ carries the dimensionality $[x]$. Consequently, we can define a bulk property $f^0$ as $f^0 = \sum_k \chi_k f_k^0 + \delta_r f_r^0$ where $f^0$ represents any property of the flow of arbitrary tensorial order. Then by summing Eq. (2.9) for $k = f, d$ and Eq. (2.10), one obtains the bulk formulation of multiphase flows, namely,

$$ \partial_t f^0 + \nabla \cdot (f^0 \mathbf{u}^0 - \mathbf{\Phi}^0) = s^0. \qquad (2.11) $$

where $f^0 \mathbf{u}^0 = \sum_k \chi_k f_k^0 \mathbf{u}_k^0 + \delta_r f_r^0 \mathbf{u}_r^0$, $\mathbf{\Phi}^0 = \sum_k \chi_k \mathbf{\Phi}_k^0 + \delta_r \mathbf{\Phi}_{\|}^0$ and $s^0 = \sum_k \chi_k s_k^0 + \delta_r s_{\|}^0$. In the literature, the single-fluid formulation is characterized by summing only the quantities within the two phases, while the interfacial components are treated as source terms (Tryggvason et al., 2011; Morel, 2015). Nevertheless, we want to point out here that our definition of $f^0$ yields a straightforward transport equation for $f^0$ thereby ensuring the consistency of the entire system of equations.

## 2.4. The averaged conservation equations

In this study, we employ the ensemble average technique to establish the averaged conservation equations. This method is just one of several averaging approaches, including the volume average method (Jackson, 1997) and time averaging (Ishii and Hibiki, 2010). Although these techniques differ, they ultimately produce the same set of averaged equations, as evidenced by comparing volume averaging theory with ensemble averaging theory (Jackson, 1997; Zhang and Prosperetti, 1997). However, for slightly non-uniform suspensions, ensemble averaging is better suited for deriving the averaged equations because it does not require assumptions about the characteristic length scale of the spatial filter (Lhuillier, 1992). In the following we recall some properties of the ensemble average operator. Let, $P(\mathcal{F})$ be the probability density function that describes the probability of finding the flow in the configuration $\mathcal{F}$. We note $d\mathcal{P} = P(\mathcal{F})d\mathcal{F}$ the probable number of flows located in the incremental region of the phase space $d\mathcal{F}$ around the point $\mathcal{F}$. It follows from this definition, that the ensemble average of an arbitrary local property $f^0(\mathbf{x}, t; \mathcal{F})$ defined on the whole space $\Omega$, is,

$$ f(\mathbf{x}, t) = \langle f^0 \rangle (\mathbf{x}, t) = \int f^0(\mathbf{x}, t; \mathcal{F}) d\mathcal{P}. \qquad (2.12) $$
Note that we dropped the super script 0 on $f(x,t)$ to indicate that this is an averaged quantity. The macroscopic variables are averaged over all $\mathcal{F}$, and therefore depend only on $x$ and $t$. Thus, we omit the arguments of the averaged fields, as this notation eliminates any potential ambiguity. The ensemble average operator is assumed to satisfy the following properties (Drew, 1983)

$$ \begin{aligned} \langle f^0 + h^0 \rangle &= f + h, \\ \langle \partial_t f^0 \rangle &= \partial_t f, \end{aligned} \qquad \langle \langle f^0 \rangle h^0 \rangle = fh, \qquad (2.13) $$

were $f^0$ and $h^0$ are two arbitrary Eulerian fields defined over $\Omega$. The first two relations are called the Reynolds' rules, the third one is the Leibniz' rule and the last one, is the Gauss' rule (Drew, 1983). Additionally, for any phase quantity defined in $\Omega_k$ we introduce the definition,

$$ \phi_k f_k(x, t) = \langle \chi_k f_k^0 \rangle, \qquad (2.14) $$

where $\phi_k(x, t) = \langle \chi_k \rangle$ is the probability of finding the phase $k$ at the location $x$ and time $t$. And $f_k$ is the average of the field $f_k^0$ conditioned on the presence of the phase $k$ in the configuration $\mathcal{F}$ at $x$ and time $t$. Equally, for interface quantities we have

$$ \phi_{\Gamma} f_{\Gamma}(x, t) = \langle \delta_{\Gamma} f_{\Gamma}^{0} \rangle, \qquad (2.15) $$

with $\phi_{\Gamma} = \langle \delta_{\Gamma} \rangle$ the interface area per unit volume, also called the specific area at the point $x$ at time $t$. Here, $f_{\Gamma}$ is the average of $f_{\Gamma}^0$ conditioned on the presence of an interface in the configuration $\mathcal{F}$ at $x$ and time $t$. Additionally, we define the field of fluctuation of a given quantity around its mean as,

$$ f'(x, t, \mathcal{F}) = f^0(x, t, \mathcal{F}) - f(x, t). \qquad (2.16) $$

This relation applies to phase averaged quantities such that $f_k' = f_k^0 - f_k$ and $f_{\Gamma}' = f_{\Gamma}^0 - f_{\Gamma}$.

Applying the ensemble average on Eqs. (2.9) and (2.10) and considering the properties from Eqs. (2.13) to (2.16), yields the general form of the averaged equations of multiphase flows, namely,

$$ \partial_t(\phi_k f_k) + \nabla \cdot (\phi_k f_k \mathbf{u}_k + \mathbf{\Phi}_k^{\text{eq}}) = \phi_k s_k + \langle \delta_{\Gamma} [\mathbf{\Phi}_k^0 + f_k^0 (\mathbf{u}_R^0 - \mathbf{u}_k^0)] \rangle \cdot \mathbf{n}_k, \qquad (2.17) $$

$$ \partial_t(\phi_{\Gamma} f_{\Gamma}) + \nabla \cdot (\phi_{\Gamma} f_{\Gamma} \mathbf{u}_{\Gamma} + \mathbf{\Phi}_{\Gamma}^{\text{eq}}) = \phi_{\Gamma} s_{\Gamma} - \langle \delta_{\Gamma} \|\| \mathbf{\Phi}_{\Gamma}^{0} + f_{\Gamma}^{0} (\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{k}^{0}) \|\| \mathbf{n}_{\Gamma} \rangle, \qquad (2.18) $$

$$ \text{with,} \\ \mathbf{\Phi}_k^{\text{eq}} = \langle \chi_k f_k' \mathbf{u}_k' \rangle - \phi_k \mathbf{\Phi}_k, \qquad \mathbf{\Phi}_R^{\text{eq}} = \langle \delta_R f_R' \mathbf{u}_R' \rangle - \phi_R \mathbf{\Phi}_R. \qquad (2.19) $$

These equations are to be solved for the averaged field $f_f$, $f_d$, and $f_{\Gamma}$. The main differences between these equations and their microscale counterparts (Eqs. (2.7) and (2.8)) are: (1) The unknowns are now averaged quantities, (2) factors $\phi_k$ and $\phi_{\Gamma}$ are introduced in front of all the terms, and (3) the additional terms $\langle \chi_k f_k' \mathbf{u}_k' \rangle$ and $\langle \delta_{\Gamma} f_{\Gamma}' \mathbf{u}_{\Gamma}' \rangle$ appear, representing the covariance between the conserved quantity ($f_k$ or $f_{\Gamma}$) and the local velocities.

It is important to highlight that the two-fluid model fails to adequately distinguish between the two phases, as evidenced by the symmetry $k = f$ and $d$ in the aforementioned equations. This symmetry does not hold physically because the dispersed phase possesses a distinct topological nature compared to the continuous phase. In a dispersed two-phase flow system, the closure terms are typically expressed as functions of the Lagrangian properties of the particles (Jackson, 2000). In contrast, the current system of equations yields continuously averaged quantities, which are not directly connected to the Lagrangian properties. Therefore, in the subsequent section, we will introduce a kinetic-based model specifically devoted to the dispersed phase. As illustrated below, the equations governing the dispersed phase are more comprehensive as they bear a resemblance to the equations governing a single particle.

### 3. Lagrangian equations for the dispersed phase

Lagrangian-based modeling of the dispersed phase has been explored extensively in numerous studies (Buyevich and Shchelchkova, 1979; Lhuillier and Nozières, 1992; Simonin, 1996; Zhang and Prosperetti, 1994a,b, 1997; Jackson, 1997; Zaepffel, 2011). However, these studies predominantly focus on solid particles (Buyevich and Shchelchkova, 1979; Lhuillier and Nozières, 1992; Simonin, 1996; Zhang and Prosperetti, 1994a; Jackson, 1997) or non-deformable spherical fluid inclusions (Zhang and Prosperetti, 1994b; Zaepffel, 2011). In this work, we aim to address a more general scenario by considering fluid particles with arbitrary shapes and surface properties. Thus, we propose a Lagrangian-based model for each inclusion capable of describing the dispersed phase with arbitrary accuracy.

#### 3.1. Fundamental properties of the dispersed phase

We define the mass $m_\alpha$, the position of center of mass $x_\alpha$ and the momentum $p_\alpha$ as

$$ m_{\alpha}(t, \mathcal{F}) = \int_{\Omega_{\alpha}(\mathcal{F}, t)} \rho_{d} d\Omega, \qquad (3.1) $$

$$ \mathbf{x}_{\alpha}(t, \mathcal{F}) = \frac{1}{m_{\alpha}} \int_{\Omega_{\alpha}(\mathcal{F}, t)} \rho_{d} \mathbf{x} d\Omega, \qquad (3.2) $$

$$ \mathbf{p}_{\alpha}(t, \mathcal{F}) = \int_{\Omega_{\alpha}(\mathcal{F}, t)} \rho_{d} \mathbf{u}_{d}^{0} d\Omega. \qquad (3.3) $$

Note that while the mass of the interface might be considered in the center of mass calculation, in practice, the mass within the volume is typically much greater than the mass on the interface. Therefore, we have opted to use the more straightforward definition considering only the mass within the particle volume. Subsequently, we define the velocity of the particle center of mass as

$$ \mathbf{u}_{\alpha}(t, \mathcal{F}) = \frac{d\mathbf{x}_{\alpha}}{dt}. \qquad (3.4) $$
We also define the particle's internal relative motions or the inner velocity as $\mathbf{w}_d^0 = \mathbf{u}_d^0 - \mathbf{u}_\alpha$. Similarly we define $\mathbf{w}_\Gamma^0$ as $\mathbf{w}_\Gamma^0 = \mathbf{u}_\Gamma^0 - \mathbf{u}_\alpha$. Replacing $\mathbf{x}_\alpha$ by its definition (Eq. (3.2)) in Eq. (3.4) we obtain

$$ \mathbf{u}_\alpha(t, \mathcal{F}) = \frac{1}{m_\alpha} \frac{d}{dt} \left( \int_{\Omega_\alpha(\mathcal{F}, t)} \rho_d \mathbf{x} d\Omega \right) - \frac{1}{m_\alpha^2} \frac{d}{dt} \left( \int_{\Omega_\alpha(\mathcal{F}, t)} \rho_d d\Omega \right) \int_{\Omega_\alpha(\mathcal{F}, t)} \rho_d \mathbf{x} d\Omega. \quad (3.5) $$

As will become clear in Section 3.2, using the Reynolds transport theorem (Eq. (3.10)) for both terms in parentheses and the definition of $\mathbf{x}_\alpha(t, \mathcal{F})$ in the last term gives

$$ \mathbf{u}_\alpha(t, \mathcal{F}) = \frac{1}{m_\alpha} \int_{\Omega_\alpha(\mathcal{F}, t)} \left[ \partial_t(\mathbf{x}\rho_d) + \nabla \cdot (\mathbf{u}_d^0 \mathbf{x}\rho_d) \right] d\Omega + \frac{1}{m_\alpha} \oint_{\Gamma_\alpha(\mathcal{F}, t)} \mathbf{x}\rho_d(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma \\ - \frac{1}{m_\alpha^2} \int_{\Omega_\alpha(\mathcal{F}, t)} \partial_t(\rho_d) + \nabla \cdot (\mathbf{u}_d^0 \rho_d) d\Omega - \frac{\mathbf{x}_\alpha}{m_\alpha} \oint_{\Gamma_\alpha(\mathcal{F}, t)} \rho_d(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma \quad (3.6) $$

Making use of the conservation of mass, the first term on the second line of Eq. (3.6) vanishes. This yields,

$$ \mathbf{u}_\alpha(t, \mathcal{F}) = \frac{1}{m_\alpha} \int_{\Omega_\alpha(\mathcal{F}, t)} [\partial_t(\mathbf{x}\rho_d) + \nabla \cdot (\mathbf{u}_d^0 \mathbf{x}\rho_d)] d\Omega + \frac{1}{m_\alpha} \oint_{\Gamma_\alpha(\mathcal{F}, t)} \mathbf{x}\rho_d(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma - \frac{\mathbf{x}_\alpha}{m_\alpha} \oint_{\Gamma_\alpha(\mathcal{F}, t)} \rho_d(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma \quad (3.7) $$

Then by considering the mass conservation for the first term and noticing that $\nabla \mathbf{x} = \delta \mathbf{x}$ gives

$$ \mathbf{u}_\alpha(t, \mathcal{F}) = \frac{1}{m_\alpha(t, \mathcal{F})} \left( \mathbf{p}_\alpha(t, \mathcal{F}) + \oint_{\Gamma_\alpha(\mathcal{F}, t)} \rho_d \mathbf{r}(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma \right), \quad (3.8) $$

where $\mathbf{r}(\mathbf{x}, t, \mathcal{F}) = \mathbf{x} - \mathbf{x}_\alpha(t, \mathcal{F})$ is the distance between any point inside $\Omega_\alpha$, to the center of mass of the particle $\alpha$. In Eq. (3.8), it can be observed that the first component of the velocity represents the linear momentum divided by the mass of the particle. This corresponds to the mass-averaged velocity over the volume of the particle. The second term in Eq. (3.8) arises from the contribution of anisotropic mass transfer across the surface of the particle. This mass transfer leads to the motion of the particle's center of mass, thereby contributing to the total velocity. To illustrate this concept, let us consider a fixed drop with no momentum lying over a very hot plate. In this scenario, we assume that the plate is sufficiently hot to induce evaporation, specifically on the bottom portion of the drop. Hence, under the effect of an anisotropic evaporation flux one may expect the second term to be non-negligible. Consequently, the center of mass of the drop has a non-zero velocity in the opposite direction of the plate, even though the momentum is assumed to be zero. Note that Eq. (3.8) generalizes usual expression of the center of mass velocity whom neglect the second term. In the following, for the sake of brevity we discard the dependency on $t$ and $\mathcal{F}$ on the notations for all Lagrangian quantities denoted by the subscript $\alpha$ and in particular $\Gamma_\alpha$ and $\Omega_\alpha$. Nevertheless, the reader must understand that all Lagrangian quantities and integration domains subscribed by $\alpha$ are time and configuration-dependent.

## 3.2. Lagrangian conservation laws

We assign to a particle indexed $\alpha$, occupying the domain $\Omega_\alpha$ (see Fig. 1) an arbitrary Lagrangian property $q_\alpha$ defined by $q_\alpha = \int_{\Omega_\alpha} f_d^0 d\Omega$. The Reynolds transport theorem for the (material) particle volume can be written as (Leal, 2007),

$$ \frac{d}{dt} \int_{\Omega_\alpha} f_d^0 d\Omega = \int_{\Omega_\alpha} \partial_t f_d^0 d\Omega + \oint_{\Gamma_\alpha} f_d^0 \mathbf{u}_\Gamma^0 \cdot \mathbf{n}_d d\Gamma. \quad (3.9) $$

Adding and subtracting $f_d^0 \mathbf{u}_d^0 \cdot \mathbf{n}_d$ in the last integral of this equation and using the divergence theorem, yields

$$ \frac{d}{dt} \int_{\Omega_\alpha} f_d^0 d\Omega = \int_{\Omega_\alpha} [\partial_t f_d^0 + \nabla \cdot (f_d^0 \mathbf{u}_d^0)] d\Omega + \oint_{\Gamma_\alpha} f_d^0 (\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma. \quad (3.10) $$

By substituting Eq. (2.7) into Eq. (3.10) and using the divergence theorem we obtain the conservation law of the quantity $q_\alpha$

$$ \frac{d q_\alpha}{dt} = \int_{\Omega_\alpha} s_d^0 d\Omega + \oint_{\Gamma_\alpha} [f_d^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) + \mathbf{\Phi}_d^0] \cdot \mathbf{n}_d d\Gamma. \quad (3.11) $$

The first term on the right-hand side of Eq. (3.11) accounts for the total contribution of the source term $s_d^0$ to the particle $\alpha$, while the second term is the surface integration of the phase exchange flux $f_d^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0)$ and the non-convective flux $\mathbf{\Phi}_d^0$.

Similarly, we define $q_{\Gamma\alpha} = \oint_{\Gamma_\alpha} f_d^0 d\Gamma$ as an integrated surface property of the particle $\alpha$. To differentiate time-varying surface integrals with respect to time, we make use of the surface Reynolds transport theorem, which reads (Morel, 2015, Appendix B)

$$ \frac{d}{dt} \oint_{\Gamma_\alpha} f_\Gamma^0 d\Gamma = \oint_{\Gamma_\alpha} \left[ \frac{D\Gamma f_\Gamma^0}{Dt} + f_\Gamma^0 \nabla_\parallel \cdot \mathbf{u}_\Gamma^0 \right] d\Gamma. \quad (3.12) $$

where $\frac{D\Gamma}{Dt} = \partial_t + \mathbf{u}_\Gamma^0 \cdot \nabla_\parallel$ is the material derivative operator on the surface of the particle. Upon expanding the acceleration term we obtain the expression,

$$ \frac{d}{dt} \oint_{\Gamma_\alpha} f_\Gamma^0 d\Gamma = \oint_{\Gamma_\alpha} [\partial_t f_\Gamma^0 + \nabla_\parallel \cdot (\mathbf{u}_\Gamma^0 f_\Gamma^0)] d\Gamma. \quad (3.13) $$

Then inserting Eq. (2.8) into Eq. (3.13) we obtain

$$ \frac{d}{dt} q_{\Gamma\alpha} = \oint_{\Gamma_\alpha} \nabla_\parallel \cdot \mathbf{\Phi}_{\Gamma\parallel}^0 d\Gamma + \oint_{\Gamma_\alpha} s_\Gamma^0 d\Gamma - \oint_{\Gamma_\alpha} [f_k^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_k^0) + \mathbf{\Phi}_k^0] \cdot \mathbf{n}_k d\Gamma. \quad (3.14) $$

The surface divergence theorem applied to closed surfaces reads (Nadim, 1996)

$$ \oint_{\Gamma_\alpha} \nabla_\parallel \cdot \mathbf{F} d\Gamma = \oint_{\Gamma_\alpha} \mathbf{F} \cdot \mathbf{n} (\nabla \cdot \mathbf{n}) d\Gamma, \quad (3.15) $$
where F is an arbitrary field. This theorem demonstrates that any surface property parallel to the tangential plane of Γ, such as $\Phi_{\Gamma\parallel}$, fulfills the condition $\oint_{\Gamma_a} \nabla_{\parallel} \cdot \Phi_{\Gamma\parallel}^0 d\Gamma = 0$. Therefore Eq. (3.16) yields,

$$ \frac{d}{dt} q_{\Gamma\alpha} = \oint_{\Gamma_a} s_{\Gamma}^0 d\Gamma - \oint_{\Gamma_a} \|f_k^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_k^0) + \Phi_k^0\| \cdot \mathbf{n}_k d\Gamma. \quad (3.16) $$

In Eq. (3.11), the non-convective flux within the dispersed phase appears on the right hand side. While this is not problematic for fluid particles for which the internal stress is defined, the stress within rigid particles is not defined (Batchelor, 1970). Hence to remove the internal stress from Eq. (3.11), we treat the particle's volume and surface as a single entity and derive a conservation equation for $Q_\alpha = q_\alpha + q_{\Gamma\alpha}$. By summing Eqs. (3.11) and (3.16) we directly obtain

$$ \frac{d}{dt} Q_{\alpha} = \int_{\Omega_{\alpha}} s_{\alpha}^{0} d\Omega + \oint_{\Gamma_{\alpha}} s_{\Gamma}^{0} d\Gamma + \oint_{\Gamma_{\alpha}} [f_{f}^{0}(\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{f}^{0}) + \Phi_{f}^{0}] \cdot \mathbf{n}_{d} d\Gamma. \quad (3.17) $$

This equation is the general form of the linear conservation law for the quantity $Q_\alpha$. It applies to any particle immersed in a continuous phase following the conservation laws given by Eqs. (2.7) and (2.8). We refer to this equation as the zeroth-order conservation equation of $f_d^0$, or alternatively, the linear conservation law of $f_d^0$, for the particle $\alpha$. We would like to highlight that due to the consideration of closed surface, the non-convective flux $\Phi_{\Gamma\parallel}^0$ does not appear in Eq. (3.17). Consequently, for the conservation of linear momentum, the surface stresses do not contribute to the momentum balance of a particle. As a result, surface tension or surface viscous stresses do not directly affect the linear momentum balance.<fnref n="1" /> This property has already been demonstrated by Hesla et al. (1993) who showed that the surface tension force does not contribute to the linear and angular momentum balance. Here, we have provided the general proof that the interfacial non-convective flux $\Phi_{\Gamma\parallel}^0$, which is present at the local scale according to Eq. (2.8), does not contribute to the zeroth-order conservation law of a particle with a closed surface. It is important to note that this conclusion is valid only if the non-convective flux is parallel to the interface, which is what we assumed here.

### 3.3. Higher order moment equations

Because $f_d^0$ and $f_\Gamma^0$ are not always constant over the volumes and surfaces of the particles, it is interesting to introduce the first moment of the quantities $f_d^0$ and $f_\Gamma^0$. They are defined as $q_\alpha^{(1)} = \int_{\Omega_\alpha} \mathbf{r} f_d^0 d\Omega$ and $q_{\Gamma\alpha}^{(1)} = \oint_{\Gamma\alpha} \mathbf{r} f_\Gamma^0 d\Gamma$. In general, the first moments $q_\alpha^{(1)}$ and $q_{\Gamma\alpha}^{(1)}$ hold significant importance when considering particles with high internal gradients, i.e. when $\nabla f_d^0$ or $\nabla f_\Gamma^0$ are non-negligible at the scale of one particle. Typically, when considering the motion of a fiber, one must consider the angular momentum balance (Guazzelli and Morris, 2011), which corresponds to the antisymmetric part of the first moment of momentum. It is possible to differentiate the moments $q_\alpha^{(1)}$ and $q_{\Gamma\alpha}^{(1)}$ with respect to time to obtain their conservation laws. We use Eq. (3.10) to describe the evolution of $q_\alpha^{(1)}$ within time. This yields,

$$ \frac{d}{dt} q_{\alpha}^{(1)} = \int_{\Omega_{\alpha}} [\partial_{t}(f_{d}^{0}\mathbf{r}) + \nabla \cdot (f_{d}^{0}\mathbf{r}\mathbf{u}_{d}^{0})] d\Omega + \oint_{\Gamma_{\alpha}} f_{d}^{0}\mathbf{r}(\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{d}^{0}) \cdot \mathbf{n}_{d} d\Gamma. \quad (3.18) $$

The first term on the right-hand side may be rewritten as

$$ \int_{\Omega_{\alpha}} [\partial_{t}(\mathbf{r} f_{d}^{0}) + \nabla \cdot (\mathbf{u}_{d}^{0} \mathbf{r} f_{d}^{0})] d\Omega = \int_{\Omega_{\alpha}} \mathbf{r} [\partial_{t} f_{d}^{0} + \nabla \cdot (f_{d}^{0} \mathbf{u}_{d}^{0})] d\Omega + \int_{\Omega_{\alpha}} f_{d}^{0} [\partial_{t} \mathbf{r} + (\mathbf{u}_{d}^{0} \cdot \nabla) \mathbf{r}] d\Omega \quad (3.19) $$

Substituting Eq. (2.7) in the first integral on the right-hand side, and considering the relation, $\partial_t \mathbf{r} + (\mathbf{u}_d^0 \cdot \nabla)\mathbf{r} = -\frac{d}{dt}\mathbf{x}_\alpha + \mathbf{u}_d^0 = \mathbf{w}_d^0$, for the second integral yields

$$ \frac{d}{dt} q_{\alpha}^{(1)} = \int_{\Omega_{\alpha}} \mathbf{r} (s_{d}^{0} + \nabla \cdot \mathbf{\Phi}_{d}^{0}) d\Omega + \int_{\Omega_{\alpha}} f_{d}^{0} \mathbf{w}_{d}^{0} d\Omega + \oint_{\Gamma_{\alpha}} f_{d}^{0} \mathbf{r} (\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{d}^{0}) \cdot \mathbf{n}_{d} d\Gamma \quad (3.20) $$

Then using relation $\int_{\Omega_\alpha} \mathbf{r} \nabla \cdot \mathbf{\Phi}_d^0 d\Omega = \oint_{\Gamma_\alpha} \mathbf{r} \mathbf{\Phi}_d^0 \cdot \mathbf{n}_d d\Gamma - \int_{\Omega_\alpha} \mathbf{\Phi}_d^0 d\Omega$ we get

$$ \frac{d}{dt} q_{\alpha}^{(1)} = \int_{\Omega_{\alpha}} (\mathbf{r} s_{d}^{0} + f_{d}^{0} \mathbf{w}_{d}^{0} - \mathbf{\Phi}_{d}^{0}) d\Omega + \oint_{\Gamma_{\alpha}} \mathbf{r} [\mathbf{\Phi}_{d}^{0} + f_{d}^{0}(\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{d}^{0})] \cdot \mathbf{n}_{d} d\Gamma. \quad (3.21) $$

Eq. (3.21) is the first order moment conservation equation for the particle $\alpha$. In Eq. (3.21), we recognize the first moment of the source term $s_d^0$, the first moment of the non convective flux term $\mathbf{\Phi}_d^0 \cdot \mathbf{n}_d$ and the first moment of phase exchange term, $f_d^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d$. Additionally, two supplementary terms appear in Eq. (3.21), namely: the volume integral of the non convective flux $\mathbf{\Phi}_d^0$, and a term related to the fluctuation of the internal velocity $f_d^0 \mathbf{w}_d^0$.

Following the same procedure, and making use of Eqs. (2.8), (3.13) and (3.15), one can equally show that

$$ \frac{d}{dt} q_{\Gamma\alpha}^{(1)} = \oint_{\Gamma_{\alpha}} (\mathbf{r} s_{\Gamma}^{0} + f_{\Gamma}^{0} \mathbf{w}_{\Gamma}^{0} - \mathbf{\Phi}_{\Gamma\parallel}^{0}) d\Gamma - \oint_{\Gamma_{\alpha}} \mathbf{r} \|\mathbf{\Phi}_{k}^{0} + f_{k}^{0}(\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{k}^{0})\| \cdot \mathbf{n}_{k} d\Gamma, \quad (3.22) $$

Observations similar to those for Eq. (3.21) can also be made for the first moment of surface Eq. (3.22). In particular, it is worth noting the presence of the surface non-convective flux $\Phi_{\Gamma\parallel}^0$ in Eq. (3.22).

For similar reason than the linear conservation equations, we sum Eqs. (3.21) and (3.22) to express the conservation equation of the total first moment $Q_\alpha^{(1)} = q_\alpha^{(1)} + q_{\Gamma\alpha}^{(1)}$, this yields

$$ \frac{d}{dt} Q_{\alpha}^{(1)} = \int_{\Omega_{\alpha}} (\mathbf{r} s_{d}^{0} + f_{d}^{0} \mathbf{w}_{d}^{0} - \mathbf{\Phi}_{d}^{0}) d\Omega + \oint_{\Gamma_{\alpha}} (\mathbf{r} s_{\Gamma}^{0} + f_{\Gamma}^{0} \mathbf{w}_{\Gamma}^{0} - \mathbf{\Phi}_{\Gamma\parallel}^{0}) d\Gamma + \oint_{\Gamma_{\alpha}} \mathbf{r} [\mathbf{\Phi}_{f}^{0} + f_{f}^{0}(\mathbf{u}_{\Gamma}^{0} - \mathbf{u}_{f}^{0})] \cdot \mathbf{n}_{d} d\Gamma. \quad (3.23) $$

<footnote n="1">¹ However, Marangoni effect can induce motion through the imbalance of tangential stress (Young, 1959). Indeed the resulting imbalance alters the stress within the continuous phase, generating a hydrodynamic force.</footnote>
Likewise, conservation laws can be derived for the nth order moments of volume and surface, i.e. for

$$ \mathbf{q}_{\alpha}^{(n)} = \int_{\Omega_{\alpha}} \underbrace{\mathbf{r}\mathbf{r}\dots\mathbf{r}\mathbf{r}}_{n \text{ times}} f_{d}^{0} d\Omega, \quad \text{and} \quad \mathbf{q}_{\Gamma\alpha}^{(n)} = \oint_{\Gamma_{\alpha}} \underbrace{\mathbf{r}\mathbf{r}\dots\mathbf{r}\mathbf{r}}_{n \text{ times}} f_{\Gamma}^{0} d\Gamma, \qquad (3.24) $$

respectively. The time derivatives of $q_{\alpha}^{(n)}$ and $q_{\Gamma\alpha}^{(n)}$ do not introduce any additional terms beyond those already present in Eqs. (3.21) and Eq. (3.22). Instead, they only involve the nth order moments of the existing terms. We provide the full derivation of $\frac{d}{dt}q_{\alpha}^{(n)}$ in Appendix C. The higher-order moments characterize the distributions of the local quantities $f_{d}^{0}$ and $f_{\Gamma}^{0}$ within the regions $\Omega_{\alpha}$ and $\Gamma_{\alpha}$, respectively. To precisely reconstruct the fields $f_{d}^{0}$ and $f_{\Gamma}^{0}$ within $\Omega_{\alpha}$ and $\Gamma_{\alpha}$, an infinite number of moments would theoretically be required. Nevertheless, as pointed out by Zhang and Prosperetti (1997, Appendix A), when the degrees of freedom of the particles are limited (solid particles, spherical droplets in stokes flows...), only a finite number of moments are necessary to reconstruct $f_{d}^{0}$ and $f_{\Gamma}^{0}$.

### 3.4. Particle-averaged equations

In the preceding subsections, we have described the dispersed phase using a Lagrangian framework. However, to ensure consistency with the Eulerian conservation equations that describe the continuous phase, it is necessary to extend the Lagrangian equations to an Eulerian description. The approach presented here follows the methodology pioneered by Lhuillier (1992). We introduce the function $\delta_{\alpha}$, which is defined as follows,

$$ \delta_{\alpha}(\mathbf{x}, \mathbf{x}_{\alpha}(t, \mathcal{F})) = \delta(\mathbf{x} - \mathbf{x}_{\alpha}(t, \mathcal{F})), \qquad (3.25) $$

where $\delta$ is the Dirac function. Note that we explicitly note the arguments $(t, \mathcal{F})$ to highlight that the position of the particle $\alpha$ is a function of time and of the flow configuration $\mathcal{F}$. Taking the time derivative on $\delta_{\alpha}$ and applying the chain rule yields (Lhuillier and Nadim, 1998, 2009)

$$ \partial_t \delta_{\alpha} + \nabla \cdot (\mathbf{u}_{\alpha} \delta_{\alpha}) = 0, \qquad (3.26) $$

where we used the identity, $\frac{\partial \delta_{\alpha}}{\partial x_{\alpha}} = -\nabla \delta_{\alpha}$ and the fact that $\mathbf{u}_{\alpha}(t, \mathcal{F})$ is not a function of $\mathbf{x}$. Eq. (3.26) does not apply in scenarios where topological changes occur, such as break-up or coalescence events. In these cases, a source term can be introduced on the right-hand side of Eq. (3.26), similar to the approach used in population balance equations, to account for the birth or death of particles (Randolph, 2012). Consider a domain containing N particles. We define the particle field for a quantity $q_{\alpha}$ as the sum of $\delta_{\alpha}q_{\alpha}$ over all particles within the domain, expressed by $\sum_{\alpha=0}^{N} \delta_{\alpha}q_{\alpha}$. Note that the formula given by Eq. (3.26) remains valid for the sum of such fields, since the operations of differentiation and summation commute. To obtain the averaged equations for the dispersed phase, we define the particle average of $q_{\alpha}$ as

$$ n_{\rho}q_{\rho}(\mathbf{x}, t) = \left\langle \sum_{\alpha} \delta_{\alpha}q_{\alpha} \right\rangle, \qquad (3.27) $$

where, $n_{\rho}(\mathbf{x}, t) = \langle \sum_{\alpha} \delta_{\alpha}q_{\alpha} \rangle$ is the probability density of finding a particle center of mass in the infinitesimal volume $d\mathbf{x}$ around $\mathbf{x}$, and $q_{\rho}(\mathbf{x}, t)$ is the average of $q_{\alpha}$ conditionally on the presence of a particle's center of mass at $\mathbf{x}$ and time $t$. To simplify the notations, we consider the shorthand (Lhuillier and Nadim, 1998),

$$ \sum_{\alpha} \delta_{\alpha} \rightarrow \delta_{\rho}, $$

such that $\langle \delta_{\rho}q_{\alpha} \rangle = \langle \sum_{\alpha} \delta_{\alpha}q_{\alpha} \rangle = n_{\rho}q_{\rho}$. Note that we used the subscript $\rho$ on $q_{\rho}$ to denote that this represents a particle-averaged field, initially derived from Lagrangian quantities. Additionally, in light of Eq. (2.16) we define the fluctuating part of a particle field $q_{\rho}$ as

$$ q'_{\alpha}(\mathbf{x}, t, \mathcal{F}) = q_{\alpha}(\mathcal{F}, t) - q_{\rho}(\mathbf{x}, t). \qquad (3.28) $$

To derive the averaged equations for the particle phase, we begin by noting that multiplying a Lagrangian quantity $q_{\alpha}$ by $\delta_{\alpha}$ yields the corresponding Eulerian field $q_{\alpha}(t, \mathcal{F})\delta_{\alpha}(\mathbf{x}, t, \mathcal{F})$, which is defined throughout the entire domain $\Omega$. Similarly, for any derivative of a Lagrangian quantity, such as $\frac{d}{dt}q_{\alpha}$, the corresponding Eulerian field is defined by multiplying $\frac{d}{dt}q_{\alpha}$ with $\delta_{\alpha}$. Given that $q_{\alpha}(t, \mathcal{F})$ and $\mathbf{u}_{\alpha}(t, \mathcal{F})$ do not depend on $\mathbf{x}$, and by using Eq. (3.26), we obtain

$$ \delta_{\alpha} \frac{d}{dt} q_{\alpha} = \partial_{t}(\delta_{\alpha} q_{\alpha}) + \nabla \cdot (\delta_{\alpha} q_{\alpha} \mathbf{u}_{\alpha}). \qquad (3.29) $$

Multiplying Eqs. (3.17) and (3.23) by $\delta_{\alpha}$, using Eq. (3.29) and applying the ensemble average (Eq. (2.12)), yields

$$ \partial_t (n_p Q_p) + \nabla \cdot (n_p Q_p \mathbf{u}_p + \langle \delta_p \mathbf{u}_\alpha' Q_\alpha' \rangle) = \left\langle \delta_p \int_{\Omega_\alpha} s_d^0 d\Omega \right\rangle + \left\langle \delta_p \oint_{\Gamma_\alpha} s_\Gamma^0 d\Gamma \right\rangle + \left\langle \delta_p \oint_{\Gamma_\alpha} [\mathbf{\Phi}_f^0 + f_f^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle, \qquad (3.30) $$

$$ \partial_t (n_p Q_p^{(1)}) + \nabla \cdot (n_p Q_p^{(1)} \mathbf{u}_p + \langle \delta_p \mathbf{u}_\alpha' Q_\alpha^{(1)'} \rangle) = \left\langle \delta_p \int_{\Omega_\alpha} (\mathbf{r} s_d^0 + f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0) d\Omega \right\rangle + \left\langle \delta_p \oint_{\Gamma_\alpha} (\mathbf{r} s_\Gamma^0 + f_\Gamma^0 \mathbf{w}_\Gamma^0 - \mathbf{\Phi}_\Gamma^{0\parallel}) d\Gamma \right\rangle \\ + \left\langle \delta_p \oint_{\Gamma_\alpha} \mathbf{r} [\mathbf{\Phi}_f^0 + f_f^0(\mathbf{u}_\Gamma^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle. \qquad (3.31) $$

The derivation of the higher moment particle-averaged equations is provided in Appendix C. The only fluxes appearing in Eqs. (3.30) and (3.31) are the fluctuation tensors $\langle \delta_p \mathbf{u}_\alpha' Q_\alpha' \rangle$ and $\langle \delta_p \mathbf{u}_\alpha' Q_\alpha' \rangle$. Therefore, the non-convective fluxes $\mathbf{\Phi}_d^0$ and $\mathbf{\Phi}_\Gamma^0$ do not play the role of macroscopic fluxes, as it is the case in Eqs. (2.17) and (2.18). Instead, they act as source terms in the first moment and higher moment equations. This distinction is the main structural differences between the Kinetic-like model (Eqs. (3.30) and (3.31)) and the two-phase flow model (Eqs. (2.17) and (2.18)). In this study, Eqs. (2.17) and (2.18) are referred to as the phase-averaged equations, while Eqs. (3.30) and (3.31) are called the particle-averaged equations.
## 4. Hybrid description of the two phases

In this section, we discuss the link between the phase averaged equations for the dispersed phase presented in Section 2 and the particle or Lagrangian averaged equations presented in Section 3.

### 4.1. Link between particle-averaged and phase-averaged equations

To model the dispersed phase, there are two distinct approaches. We can either use Eq. (2.17) with $k = d$, or we can employ the particle-averaged equations Eqs. (3.30), (3.31) and potentially the higher moments equations found in Appendix C. To better understand the physical significance of this choice and determine which formalism is more appropriate, it is essential to discuss the relationship that connects these two formalisms.

It has been demonstrated in various studies (Buyevich and Shchelchkova, 1979; Lhuillier, 1992; Zhang and Prosperetti, 1994a), that phase-averaged quantities can be expressed as a Taylor series expansion of particle-averaged quantities. The aforementioned studies used the single-particle conditionally averaged approach to demonstrate this equivalence. In this work, we follow the “distributional” approach proposed by Pähtz et al. (2025), as it clarifies the relationship between phase quantities and moment expansions prior to the application of any averaging formalism. As demonstrated by Pähtz et al. (2025) (see also Appendix D for a demonstration in the sense of distribution) we may show

$$ \chi_a f_d^0 = \delta_a \int_{\Omega_a} f_d^0 d\Omega - \nabla \cdot \left( \delta_a \int_{\Omega_a} \mathbf{r} f_d^0 d\Omega \right) + \frac{1}{2} \nabla \nabla : \left( \delta_a \int_{\Omega_a} \mathbf{r} \mathbf{r} f_d^0 d\Omega \right) - \dots \quad (4.1) $$

We recall that $\chi_a = 1$ within the particle domain $\Omega_a(\mathcal{F}, t)$ and 0 otherwise. Additionally, we have extended this approach to surface quantities in Appendix D, resulting in

$$ \delta_{\Gamma a} f_{\Gamma}^{0} = \delta_{a} \oint_{\Gamma_{a}} f_{\Gamma}^{0} d\Gamma - \nabla \cdot \left( \delta_{a} \oint_{\Gamma_{a}} \mathbf{r} f_{\Gamma}^{0} d\Gamma \right) + \frac{1}{2} \nabla \nabla : \left( \delta_{a} \oint_{\Gamma_{a}} \mathbf{r} \mathbf{r} f_{\Gamma}^{0} d\Gamma \right) - \dots \quad (4.2) $$

Likewise, we recall that $\delta_{\Gamma a} = 1$ on the particle surface $\Gamma_a(\mathcal{F}, t)$ and 0 otherwise. We can identify the zeroth, first and second order moments of $f_d^0$ and $f_{\Gamma}^0$, into Eqs. (4.1) and (4.2), respectively. Then summing Eqs. (4.1) and (4.2) we obtain,

$$ \chi_a f_d^0 + \delta_{\Gamma a} f_{\Gamma}^0 = \delta_a Q_a - \nabla \cdot (\delta_a Q_a^{(1)}) + \frac{1}{2} \nabla \nabla : (\delta_a Q_a^{(2)}) - \dots \quad (4.3) $$

It is important to observe that, even before any averaging procedure is applied, Eq. (4.3) already demonstrate the link between the dispersed phase fields, even for a single particle. Notably, this relationship is valid in a distributional sense at the local level.

The dispersed phase indicator function $\chi_d$ can be expressed as a sum of phase indicator function, $\chi_d(\mathbf{x}, t, \mathcal{F}) = \sum_a \chi_a(\mathbf{x}, t, \mathcal{F})$. Likewise, the interface indicator function $\delta_R$ can be written as $\delta_R = \sum_a \delta_{\Gamma a}$. Thus, summing Eqs. (4.1) and (4.2) over all particles in the domain and averaging over all configurations, gives the general relations that link continuous-averaged and particle-averaged fields, namely (Lhuillier, 1992; Lhuillier and Nadim, 1998; Lhuillier et al., 2000),

$$ \langle \chi_d f_d^0 \rangle = \langle \delta_p q_a \rangle - \nabla \cdot \langle \delta_p q_a^{(1)} \rangle + \frac{1}{2} \nabla \nabla : \langle \delta_p q_a^{(2)} \rangle + \dots \quad (4.4) $$

$$ \langle \delta_R f_R^0 \rangle = \langle \delta_p q_{\Gamma a} \rangle - \nabla \cdot \langle \delta_p q_{\Gamma a}^{(1)} \rangle + \frac{1}{2} \nabla \nabla : \langle \delta_p q_{\Gamma a}^{(2)} \rangle - \dots \quad (4.5) $$

Summing Eqs. (4.4) and (4.5) we obtain

$$ \langle \chi_d f_d^0 + \delta_R f_R^0 \rangle = \langle \delta_p Q_a \rangle - \nabla \cdot \langle \delta_p Q_a^{(1)} \rangle + \frac{1}{2} \nabla \nabla : \langle \delta_p Q_a^{(2)} \rangle - \dots \quad (4.6) $$

When considering an infinite number of terms in Eq. (4.6) one might eventually obtain a converged approximation of $\langle \chi_d f_d^0 + \delta_R f_R^0 \rangle$. However, it is important to note that Taylor series have what is known as a *radius of convergence* beyond which adding more terms does not necessarily improve the approximation (Appel, 2007, Chapter 1). In particular, for distances beyond a certain limit **r** the series might diverge depending on the behavior of the function $\langle \chi_d f_d^0 + \delta_R f_R^0 \rangle$ near the point **x**. For the purposes of this article, we will assume that the Taylor series has an infinite radius of convergence, although this assumption warrants further investigation.

To demonstrate the relation between the two formalisms, we follow a strategy similar to Lhuillier et al. (2000) and Lhuillier and Nadim (2009). Let $C_d$ denote the phase-averaged equation of conservation (Eq. (2.17) with $k = d$) and $C_R$ the averaged surface transport Eq. (2.18). Specifically, they are defined as follows

$$ C_d = -\partial_t \langle \chi_d f_d^0 \rangle - \nabla \cdot \langle \chi_d \Phi_d^0 - \chi_d f_d^0 \mathbf{u}_d^0 \rangle + \langle \chi_d s_d^0 \rangle + \langle \delta_R [\Phi_d^0 + f_d^0 (\mathbf{u}_R^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d \rangle, \quad (4.7) $$

$$ C_R = -\partial_t \langle \delta_R f_R^0 \rangle - \nabla \cdot \langle \delta_R f_R^0 \mathbf{u}_R^0 - \delta_R \Phi_{\parallel}^0 \rangle + \langle \delta_R s_R^0 \rangle - \langle \delta_R [\Phi_k^0 + f_k^0 (\mathbf{u}_R^0 - \mathbf{u}_k^0)] \cdot \mathbf{n}_k \rangle. \quad (4.8) $$

It should be noted from Eqs. (2.17) and (2.18) that $C_d \equiv 0$ and $C_R \equiv 0$. By applying the Taylor expansion to each term of $C_d + C_R$ as described in Eq. (4.6) yields

$$ C_d + C_R = \mathcal{M}^{(0)} - \nabla \cdot \mathcal{M}^{(1)} + \frac{1}{2} \nabla \nabla : \mathcal{M}^{(2)} - \dots = 0, \quad (4.9) $$

where the expressions for $\mathcal{M}^{(0)}$ and $\mathcal{M}^{(1)}$ are given by

$$ \mathcal{M}^{(0)} = - \langle \delta_p \frac{d}{dt} Q_a \rangle + \langle \delta_p \int_{\Omega_a} s_d^0 d\Omega \rangle + \langle \delta_p \oint_{\Gamma_a} s_R^0 d\Gamma \rangle + \langle \delta_p \oint_{\Gamma_a} [\Phi_f^0 + f_f^0 (\mathbf{u}_R^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \rangle, \quad (4.10) $$

$$ \mathcal{M}^{(1)} = - \langle \delta_p \frac{d}{dt} Q_a^{(1)} \rangle + \langle \delta_p \int_{\Omega_a} (\mathbf{r} s_d^0 + f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0) d\Omega \rangle + \langle \delta_p \oint_{\Gamma_a} (\mathbf{r} s_R^0 + f_R^0 \mathbf{w}_R^0 - \mathbf{\Phi}_{\parallel}^0) d\Gamma \rangle + \langle \delta_p \oint_{\Gamma_a} \mathbf{r} [\mathbf{\Phi}_f^0 + f_f^0 (\mathbf{u}_R^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \rangle. \quad (4.11) $$
Using Eq. (4.9), we reach one of the main conclusion of this study. We observe that $M^{(0)}$ and $M^{(1)}$ correspond to the zeroth and first-order moment equations, respectively. Additionally, as demonstrated in Appendix C the coefficient $M^{(n)}$ in Eq. (4.9) represents the nth order moment in the particle-averaged conservation equation. From Eq. (4.9) we conclude that combining Eq. (2.17) for $k = d$ and Eq. (2.18) effectively captures the particle moment equations through a Taylor expansion around the particle center of mass. Therefore, it is clear that by considering an arbitrary order of particle moments equations, one can achieve a highly accurate description of the dispersed phase.

The particle-averaged equations ($M^{(0)}...M^{(n)}$) form a system with $n$ equations, one for each moment. In contrast, the dispersed phase-averaged equations ($C_d$ and $C_r$) consist of only two equations, which aggregate all the particle-averaged equations. This indicates that the particle-averaged formalism provides more information since it yields a separate equation for each moment, as opposed to the phase-averaged equations, which are limited to just two. This enhanced level of detail is achieved by considering the topology of the dispersed phase as demonstrated in the previous sections.

## 4.2. Conservation equations

Given that the aim of this work is not only to demonstrate the relationship between particle and phase-averaged formalisms but also to establish a comprehensive framework for analyzing dispersed two-phase flows, we now present the hybrid set of conservation equations. The system of equations governing a macroscopic quantity $f$ consists of one equation for the fluid phase, which ensures the conservation of $f_f$ and $n$ equations for the dispersed phase, representing the conservation of the quantities $Q_p^{(n)}$. In its most general form, the hybrid description of $f$ can be expressed as

$$ \partial_t(\phi_f f_f) + \nabla \cdot (\phi_f f_f \mathbf{u}_f + \mathbf{\Phi}_f^{\text{eff}}) = \phi_f s_f - \left\langle \delta_p \oint_{\Gamma_a} [\mathbf{\Phi}_f^0 + f_f^0(\mathbf{u}_f^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle, \quad (4.12) $$

$$ \partial_t(n_p Q_p) + \nabla \cdot (n_p Q_p \mathbf{u}_p + \langle \delta_p \mathbf{u}_a' Q_a \rangle) = \left\langle \delta_p \int_{\Omega_a} s_d^0 d\Omega \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} s_r^0 d\Gamma \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} [\mathbf{\Phi}_f^0 + f_f^0(\mathbf{u}_r^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle, \quad (4.13) $$

$$ \partial_t(n_p Q_p^{(1)}) + \nabla \cdot (n_p \mathbf{u}_p Q_p^{(1)} + \langle \delta_p \mathbf{u}_a' Q_a^{(1)} \rangle) = \left\langle \delta_p \int_{\Omega_a} (\mathbf{r} s_d^0 + f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0) d\Omega \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} (\mathbf{r} s_r^0 + f_r^0 \mathbf{w}_r^0 - \mathbf{\Phi}_r^0) d\Gamma \right\rangle \\ + \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} [\mathbf{\Phi}_f^0 + f_f^0(\mathbf{u}_f^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle, \quad (4.14) $$

...

where the effective continuous phase non-convective flux term reads

$$ \mathbf{\Phi}_f^{\text{eff}} = \langle \chi_f f_f' \mathbf{u}_f' \rangle - \langle \chi_f \mathbf{\Phi}_f^0 \rangle - \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} [\mathbf{\Phi}_f^0 + f_f^0 (\mathbf{u}_f^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle + \nabla \cdot [\dots]. \quad (4.15) $$

The only difference in the conservation equation for the continuous phase (between Eqs. (2.17) and (4.12)) is the expansion of the exchange term $\langle \delta_r [\mathbf{\Phi}_f^0 + f_f^0(\mathbf{u}_f^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d \rangle$ into a Taylor series in a similar way to Eq. (4.5). The presence of the terms $\nabla \cdot [\dots]$ in the expression for $\mathbf{\Phi}_f^{\text{eff}}$ suggests that higher-order moments of the interphase exchange term are involved. Similarly, the ellipsis below Eq. (4.14) implies that an arbitrary number of dispersed phase moment equations can be introduced. In this format, it becomes evident that the exchange term on the right-hand side of Eq. (4.12) is identical to that on the right-hand side of Eq. (4.13). Additionally, in the effective flux $\mathbf{\Phi}_f^{\text{eff}}$ the exchange term from Eq. (4.14) appears, and this property continues for higher-order moments. Consequently, the zeroth order exchange term in the equation for $Q_a^{(0)}$ plays the role of a source term for $f_f$, while the first and higher order exchange terms act as a source into the higher moment equations, and contribute to the effective non-convective fluxes for $f_f$.

The system of equations presented here offers a clear understanding of the roles of the dispersed phase non-convective flux terms, $\mathbf{\Phi}_d$ and $\mathbf{\Phi}_r$, in the particle phase conservation equation. As evidenced in Eq. (4.13), $\mathbf{\Phi}_d$ and $\mathbf{\Phi}_r$ do not influence the lowest order particle-phase averaged conservation equation, i.e. the equation of $n_p Q_p$. However, Eq. (2.17) (for $k = d$) and Eq. (2.18), show that the phase-averaged quantities $f_d$ and $f_r$, are affected by the non-convective fluxes at the particle surface and internally, since $\mathbf{\Phi}_d$ and $\mathbf{\Phi}_r$ appear in these equations. This may seem contradictory at first, but it is important to note that $\mathbf{\Phi}_d^0$ and $\mathbf{\Phi}_r^0$ serve as source terms in the conservation equations for higher moments such as $Q_p^{(1)}...Q_p^{(n)}$, which are related to $f_d$ and $f_r$ through Eq. (4.6). In summary, the non-convective fluxes $\mathbf{\Phi}_d$ and $\mathbf{\Phi}_r$ are not explicitly related to $Q_p$, regardless of particle nature or volume fraction. Instead, their influence on $Q_p$ is mediated through the closure terms in Eq. (4.13), which may depend on the higher moments $Q_p^{(1)}...Q_p^{(n)}$ or other higher-order particle-related moments. These moments, in turn, explicitly depend on $\mathbf{\Phi}_d$ and $\mathbf{\Phi}_r$ as indicated by Eq. (4.14).

The hybrid model consists of $n + 1$ equations for $n + 1$ unknowns, specifically $f_f$, $Q_p^{(1)}...Q_p^{(n)}$.<fnref n="2" /> The closure problem then involves deriving explicit expressions for all terms of Eqs. (4.12) to (4.14), of the form $\langle \dots \rangle$ in terms of the problem's unknowns, i.e. $f_f$, $Q_a$, $Q_a^{(1)}...Q_a^{(n)}$. Once these expressions are determined, the equations can be solved.

## 4.3. Discussion

The last question we would to address is: in what circumstances is the inclusion of higher moments, such as $Q_p^{(n)}$, necessary within this system of equations? There are three primary reasons why the moment $Q_p^{(n)}$ may be required:

1. First when $Q_p^{(n)}$ provides critical information sought in the analysis. For example, consider the case of fiber orientation in a flow, which corresponds to the second moment of the particle mass distribution. This information might be the central objective of the study, especially when the goal is to understand how fiber orientation evolves in industrial applications like composite materials (Advani and Tucker, 1987).

<footnote n="2">² $\phi_f$ and $n_p$ can also be treated as unknowns by adding corresponding equations, achieved by setting $f_f = 1$ and $Q_a^{(n)} = 1$ in Eqs. (4.12) and (4.13).</footnote>
<table><caption>**Table 1**

Definition of the physical quantities and local constitutive laws for phases k.</caption><thead><tr><th></th><th>Conservation law</th><th>Mass</th><th>Momentum</th></tr></thead><tbody><tr><td>Conserved quantity</td><td>$f_k^0$</td><td>$\rho_k$</td><td>$\rho_k \mathbf{u}_k^0$</td></tr><tr><td>Source term</td><td>$s_k^0$</td><td>0</td><td>$\rho_k \mathbf{g}$</td></tr><tr><td>Diffusive flux</td><td>$\Phi_k^0$</td><td>0</td><td>$\sigma_k^0 = -p_k^0 + \mu_k (\nabla \mathbf{u}_k^0 + \nabla \mathbf{u}_k^0)$</td></tr><tr><td>Surface diffusive flux</td><td>$\Phi_r^0$</td><td>0</td><td>$\sigma_r^0 = \gamma(\delta - \mathbf{n}\mathbf{n})$</td></tr></tbody></table>

2. Second $Q_{\rho}^{(n)}$ is essential for accurately modeling closure terms. Again, fiber orientation serves as a relevant example. Indeed in dilute flows of axisymmetric fibers within the Stokes regime the force acting on the fiber (which represents the closure term in the momentum equation) is dependent on the fiber orientation (Kim and Karrila, 2013).

3. Finally, in the continuous phase conservation Eq. (4.12), an infinite number of moments are involved. Hence, the last reason why one may need to obtain such higher moments is because they may be required directly by Eq. (4.12).

In the following, we focus on dispersed fluid-fluid suspensions, emphasizing the role of mass and momentum moments in their dynamics.

## 5. Averaged equations for dispersed fluid-fluid flows with surface tension

We consider a monodisperse fluid-fluid suspension of droplets (or bubbles) which are not necessarily spherical with volume $v_{\rho}$, without mass transfer. Both the dispersed and continuous phases are treated as incompressible Newtonian fluids, characterized by constant viscosities $\mu_k$ and densities $\rho_k$. The surface tension at the interface between the two fluids is denoted by $\gamma$ which is not necessarily constant. We only consider in the next section averaged mass and momentum conservation equations to describe the motions of the phases. The conserved physical quantities relevant to the problem are summarized in Table 1.

In the following, we consider the influence of Marangoni effects on the suspension dynamics. Surfactants and/or non-uniform temperature gradients typically cause variations in the surface tension coefficient at the droplet interfaces. One must solve a transport equation for either the temperature field or the surfactant concentration in order to determine the proper surface tension distribution on the droplet surface (Subramanian, 1985; Leal, 2007). This can be done using the surface transport equations derived in Section 2.2. However, as these complexities fall outside the scope of this work, we assume that the surface tension coefficient is known in advance, having been determined from a separate problem solved a priori.

After presenting the higher-order mass and momentum moments required to describe the particle rotation and deformation we derive the averaged equations for mass, momentum and first moment of momentum for a fluid-fluid suspension. The presentation is inspired by Lhuillier and Nadim (2009), while extending it beyond the scope of solid particles. Then we turn our attention to the decomposition of stress tensor in the momentum exchange terms. The section concludes with a discussion of the symmetry properties of the effective stress tensor, employing arguments analogous to those presented by Lhuillier (1996).

### 5.1. Higher-order mass and momentum moments

To describe the dispersed phase in addition to the quantities defined in Section 3 we define the second-order moment of mass and the first-order moment of momentum as

$$ \mathbf{M}_{\alpha} = \int_{\Omega_{\alpha}} \rho_{d} \mathbf{r} \mathbf{r} d\Omega \quad \text{and} \quad \mathbf{P}_{\alpha} = \int_{\Omega_{\alpha}} \rho_{d} \mathbf{r} \mathbf{u}_{d}^{0} d\Omega. \qquad (5.1) $$

respectively. Note that the tensor $\mathbf{M}_{\alpha}$ plays a role analogous to the inertia tensor $\mathbf{I}_{\alpha}$ in solid mechanics. The two are related by the expression $\mathbf{I}_{\alpha} = (\delta : \mathbf{M}_{\alpha})\delta - \mathbf{M}_{\alpha}$. The second order tensor $\mathbf{M}_{\alpha}$ represents the second moment of the mass distribution relative to the particle center of mass. Similarly, the tensor $\mathbf{P}_{\alpha}$ is the first moment of the momentum distribution within the particle volume. To offer a clearer physical interpretation of the moment of momentum tensor, we decompose $\mathbf{P}_{\alpha}$ into three distinct components

$$ \mathbf{P}_{\alpha} = \mathbf{S}_{\alpha} + \mathbf{T}_{\alpha} + \mathbf{P}_{\alpha} \delta, \qquad (5.2) $$

where $\mathbf{S}_{\alpha} = \frac{1}{2}(\mathbf{P}_{\alpha} + \mathbf{P}_{\alpha}^{\dagger} - \frac{2}{3}(\delta : \mathbf{P}_{\alpha})\delta)$ is the symmetric traceless part, representing the deformation or stretching of momentum, $\mathbf{T}_{\alpha} = \frac{1}{2}(\mathbf{P}_{\alpha} - \mathbf{P}_{\alpha}^{\dagger})$ is the antisymmetric part, associated with the angular momentum, and $\mathbf{P}_{\alpha} = \frac{1}{3}\delta : \mathbf{P}_{\alpha}$ is a scalar that quantifies the momentum of expansion or compression of the droplet. The angular momentum of the particle can also be expressed as the pseudo-vector $\mu_{\alpha} = \int_{\Omega_{\alpha}} \rho_{d} \mathbf{r} \times \mathbf{u}_{d}^{0} dV$, and is related to $\mathbf{T}_{\alpha}$ via $(\mu_{\alpha})_{ij} = \epsilon_{ijk} (\mathbf{T}_{\alpha})_{jk}$ where $\epsilon_{ijk}$ is the Levi-Civita symbol. To illustrate the physical significance of these tensors, we present in Fig. 2 three representative inner velocity fields, each with its corresponding moment of momentum tensor. Injecting, $f_{d}^{0} = \rho_{d}$ in the second-order moment equation (derived in Appendix C) we obtain,

$$ \frac{d}{dt} \mathbf{M}_{\alpha} = 2(\mathbf{S}_{\alpha} + \mathbf{P}_{\alpha} \delta), \qquad (5.3) $$

which is the second-order moment of mass conservation equation assuming that the fluid within the drop is divergence free. From Eq. (5.3) we deduce that the evolution of the distribution of mass of a particle is solely determined by the of momentum $\mathbf{S}_{\alpha}$ and $\mathbf{P}_{\alpha}$. This indicates that angular momentum does not influence the evolution of the second moment of mass, a consequence of the symmetry of the tensor $\mathbf{M}_{\alpha}$, which must be preserved after differentiation with respect to time. However, this does not imply that the particle angular velocity $(\omega_{\alpha})$ does not appear in this equation. For instance, in the case of rigid body motion where $\mathbf{w}_{d}^{0} = \omega_{\alpha} \times \mathbf{r}$, we obtain the relation $2(\mathbf{S}_{\alpha} + \mathbf{P}_{\alpha} \delta)_{ij} = \epsilon_{iab}(\omega_{\alpha})_{a}(\mathbf{M}_{\alpha})_{bj} + \epsilon_{jab}(\omega_{\alpha})_{a}(\mathbf{M}_{\alpha})_{bi}$.

Now that we have described the kinematics of the particle shape, let us proceed to derive an equation for the moment of momentum. This equation is derived by injecting $Q_{\alpha}^{(1)} = \mathbf{P}_{\alpha}$ in Eq. (3.23), it reads,

$$ \frac{d}{dt} \mathbf{P}_{\alpha} - \int_{\Omega_{\alpha}} \rho_{d} \mathbf{w}_{d}^{0} \mathbf{w}_{d}^{0} d\Omega = - \int_{\Omega_{\alpha}} \sigma_{d}^{0} d\Omega - \oint_{\Gamma_{\alpha}} \gamma(\delta - \mathbf{n}\mathbf{n})d\Gamma + \oint_{\Gamma_{\alpha}} \mathbf{r} \sigma_{f}^{0} \cdot \mathbf{n} d\Gamma. \qquad (5.4) $$
<figure><img src="image_5.png" /><figcaption>**Fig. 2.** Graphical representation of the inner kinematics of an arbitrary particle under three scenarios. The arrows represent the velocity field inside the particle, $\mathbf{w}_d^0$, with the corresponding value of the moment of momentum tensor indicated below. The operator $|...|$ refers to the norm of the tensors. According to the inner velocity field: (Case 1) The particle experiences a mean deformation, resulting in non-zero stretching of momentum along the principal axis of deformation; (Case 2) The particle is rotating, leading to a non-zero angular momentum vector in the direction of rotation; (Case 3) The particle undergoes compression, resulting in a negative trace of the moment of momentum.</figcaption></figure>

In the following sections the normal vector noted $\mathbf{n}$ always refer to $\mathbf{n}_d$. The conservation equation of the angular momentum $\mu_\alpha$ is obtained by taking the double contracted product of Eq. (5.4) with $\epsilon$, which directly gives

$$ \frac{d}{dt} \mu_{\alpha} = \oint_{\Gamma_{\alpha}} \mathbf{r} \times \boldsymbol{\sigma}_{f}^{0} \cdot \mathbf{n} \, d\Gamma \qquad (5.5) $$

Note that every term on the right-hand side of Eq. (5.4) vanished due to their symmetric nature apart from the skew-symmetric part of the hydrodynamic stress, which is the hydrodynamic torque applied on the particle $\alpha$. In particular, the surface tension terms do not appear in the angular momentum balance since the tensor $\delta - \mathbf{n}\mathbf{n}$ is symmetric, which is consistent with the findings of Hesla et al. (1993). As a consequence, the surface tension does not affect the angular momentum regardless of the particle shape.

Taking the symmetric part of Eq. (5.4), and substrating the trace yields,

$$ \frac{d}{dt} \mathbf{S}_{\alpha} - \rho_{d} \int_{\Omega_{\alpha}} \left( \mathbf{w}_{d}^{0} \mathbf{w}_{d}^{0} - \frac{1}{3} (\mathbf{w}_{d}^{0} \cdot \mathbf{w}_{d}^{0}) \delta \right) d\Omega = -2\mu_{d} \int_{\Omega_{\alpha}} \mathbf{e}_{d}^{0} d\Omega - \oint_{\Gamma_{\alpha}} \gamma \left( \frac{1}{3} \delta - \mathbf{n}\mathbf{n} \right) d\Gamma + \frac{1}{2} \oint_{\Gamma_{\alpha}} \left( \mathbf{r}\boldsymbol{\sigma}_{f}^{0} + \boldsymbol{\sigma}_{f}^{0}\mathbf{r} - \frac{2}{3}(\boldsymbol{\sigma}_{f}^{0} \cdot \mathbf{r})\delta \right) \cdot \mathbf{n} \, d\Gamma, \qquad (5.6) $$

where we have introduced the rate of strain tensor for phase $k$, defined as $\mathbf{e}_k^0 = \frac{1}{2}(\nabla \mathbf{u}_k^0 + \mathbf{\nabla} \mathbf{u}_k^0)$. Here $\mathbf{\nabla}$ represents the transpose operator. On the left-hand side of Eq. (5.6), two inertial contributions can be identified: the time derivative of $\mathbf{S}_\alpha$, and the stress arising from the product of internal velocity field within the droplet. These inertial effects are counterbalanced by the terms appearing on the right-hand side of the equation, which include: the volume integral of the particle viscous stress; the surface tension moment; and the first moment of the hydrodynamic force. While surface tension effects do not influence the linear or angular momentum equations directly, they do impact the moment of momentum $\mathbf{P}_\alpha$, specifically its symmetric part $\mathbf{S}_\alpha$. Consequently, surface tension influences the hydrodynamic behavior of a particle exclusively through its effect on $\mathbf{S}_\alpha$, which is related to the shape of a particle represented by $\mathbf{M}_\alpha$, via Eq. (5.3). Inserting Eq. (5.6) in Eq. (5.3) yields,

$$ \frac{1}{2} \frac{d^2(\mathbf{M}_\alpha - \frac{1}{3}(\mathbf{M} : \delta)\delta)}{dt^2} = \rho_d \int_{\Omega_\alpha} \left( \mathbf{w}_d^0 \mathbf{w}_d^0 - \frac{1}{3} (\mathbf{w}_d^0 \cdot \mathbf{w}_d^0) \right) d\Omega - 2\mu_d \int_{\Omega_\alpha} \mathbf{e}_d^0 d\Omega - \oint_{\Gamma_\alpha} \gamma \left( \frac{1}{3} \delta - \mathbf{n}\mathbf{n} \right) d\Gamma + \frac{1}{2} \oint_{\Gamma_\alpha} \left( \mathbf{r}\boldsymbol{\sigma}_f^0 + \boldsymbol{\sigma}_f^0 \mathbf{r} - \frac{2}{3}(\boldsymbol{\sigma}_f^0 \cdot \mathbf{r})\delta \right) \cdot \mathbf{n} \, d\Gamma. \qquad (5.7) $$

Eq. (5.7) is a second-order differential equation governing the second-order mass moment that characterizes the droplet shape. In other words, it should be interpreted as an equation describing the droplet shape, represented by the traceless part of the tensor $\mathbf{M}_\alpha$. Accordingly, the right-hand side of Eq. (5.7) contains terms that promote deformation-such as the product of the internal velocity field and the traceless component of the first moment as well as terms that oppose deformation, including droplet internal viscous resistance and surface tension. Eq. (5.7) is particularly useful to compute the unknown internal stress within solid particles, in terms of surface integral, i.e. the first moment of the hydrodynamic force. This relation plays a key role in expressing the bulk stress of a suspension and ultimately leads to the effective viscosity, once a closed-form expression for the average first moment is obtained (Batchelor, 1970). In the inertial regime, and for solid particles, the tensors $\mathbf{M}_\alpha$ and the internal velocity field $\mathbf{w}_d^0$ are fully prescribed by the particle kinematics. As a result, in Eq. (5.7), $\mathbf{M}_\alpha$ and $\mathbf{w}_d^0$ appear not as unknowns but rather as known source terms. For spherical particles specifically, the inertial correction to the first force moment has been derived by Hwang and Shen (1989) and Lhuillier (1996).

## 5.2. Averaged momentum and mass conservation equations

Let us now focus on the averaged momentum equation of the continuous phase. By applying Eq. (2.17) with $f_k = \rho_f$ and $\rho_f \mathbf{u}_f$ we obtain the mass and momentum equation for the continuous phase under the two fluid formulation,

$$ (\partial_t + \mathbf{u}_f \cdot \nabla) \phi_f = -\phi_f \nabla \cdot \mathbf{u}_f \qquad (5.8) $$

$$ \phi_f \rho_f (\partial_t + \mathbf{u}_f \cdot \nabla) \mathbf{u}_f = \phi_f \rho_f \mathbf{g} + \nabla \cdot \left[ \phi_f \boldsymbol{\sigma}_f - \left\langle \chi_f \rho_f \mathbf{u}_f' \mathbf{u}_f' \right\rangle \right] - \left\langle \delta_f \boldsymbol{\sigma}_f^0 \cdot \mathbf{n} \right\rangle. \qquad (5.9) $$

To obtain the hybrid formulation one can expand the final term on the right-hand side of Eq. (5.9) using the Taylor expansion provided in Eq. (4.6). However, before doing so, it is important to examine the term $\phi_f \boldsymbol{\sigma}_f$, as it contains a non-closed contribution originating from the averaging of the strain rate tensor. Properly identifying and isolating this contribution is a necessary step prior to performing the Taylor expansion.

We begin by noting that $\phi_f \boldsymbol{\sigma}_f$ can be expressed in terms of the averaged fluid pressure $\rho_f$, and continuous phase averaged velocity ($\mathbf{u}_f$), or bulk-averaged velocity ($\mathbf{u} = \phi_f \mathbf{u}_f + \phi_d \mathbf{u}_d$). Expressing $\boldsymbol{\sigma}_f$ as a function of $\mathbf{u}_f$ instead of $\mathbf{u}$ (or vice versa) leads to two distinct forms of the momentum
equation, each associated with different formulations of the closure terms. We also review several alternative choices found in the literature (see also Jackson, 2000; Pähtz et al., 2025 for an overview focused on solid particles), emphasizing that while all these formulations are mathematically equivalent, they give rise to different closure problems. We argue that one particular formulation is probably the most suitable, especially in light of the closure terms available in the literature for non-dilute flows.

First we introduce what we refer to as the mean Newtonian stress, based either on the continuous phase averaged velocity $\mathbf{u}_f$ or on the bulk velocity $\mathbf{u}$, namely,

$$ \Sigma_f = -p_f \delta + 2\mu_f \mathbf{E}_f \qquad (5.10) $$

$$ \Sigma = -p_f \delta + 2\mu_f \mathbf{E} \qquad (5.11) $$

where $\mathbf{E}_f = [\nabla \mathbf{u}_f + (\nabla \mathbf{u}_f)^\dagger]/2$ and $\mathbf{E} = [\nabla \mathbf{u} + (\nabla \mathbf{u})^\dagger]/2$ represent the mean strain rate tensor based either on $\mathbf{u}_f$ or on $\mathbf{u}$. Additionally the ensemble averaged stress $\phi_f \sigma_f$ may be written as

$$ \phi_f \sigma_f = -\phi_f p_f \delta + 2\mu_f \phi_f \mathbf{e}_f \qquad (5.12) $$

where $\phi_f \mathbf{e}_f = \langle \chi_f (\nabla \mathbf{u}_f^0 + \nabla \mathbf{u}_f^0)^\dagger \rangle$. Inserting Eqs. (5.10) and (5.11) in Eq. (5.12) yields

$$ \phi_f \sigma_f = \phi_f \Sigma_f - \mu_f \langle \delta_f (\mathbf{u}'_f \mathbf{n} + \mathbf{n} \mathbf{u}'_f) \rangle \qquad (5.13) $$

$$ \phi_f \sigma_f = \phi_f \Sigma - \langle 2\mu_f \chi_d \mathbf{e}_d^* \rangle \qquad (5.14) $$

where we have used the relations

$$ \langle \chi_f \nabla \mathbf{u}_f^0 \rangle = \phi_f \nabla \mathbf{u}_f + \langle \delta_f \mathbf{n}_f \mathbf{u}_f' \rangle \qquad (5.15) $$

and

$$ \phi_f \mathbf{e}_f = \phi_f \mathbf{E} - \langle \chi_d \mathbf{e}_d^* \rangle, \qquad (5.16) $$

to derive Eqs. (5.13) and (5.14), respectively. In Eq. (5.14) we have introduced $\mathbf{e}_d^* = \mathbf{e}_d^0 - \mathbf{E} = \nabla(\mathbf{u}_d^0 - \mathbf{u}) + \nabla(\mathbf{u}_d^0 - \mathbf{u})$ which is the droplet internal shear rate relative to the ‘bulk’ shear rate $\mathbf{E}$. It is worth noting that Eq. (5.15) remains valid even when the interfacial rate of strain is nonzero. In contrast, Eq. (5.16) holds under the assumption that the bulk rate of strain satisfies $\mathbf{E} = \phi_f \mathbf{e}_f + \phi_d \mathbf{e}_d$, i.e. that $\langle \delta_f \mathbf{e}_f^0 \rangle = 0$. This condition is supposed to be met here because we did not consider interfacial viscosity (Nadim, 1996). The averaged momentum equation employing the stress formulation Eq. (5.13) read,

$$ \phi_f \rho_f (\partial_t + \mathbf{u}_f \cdot \nabla) \mathbf{u}_f = \phi_f (\nabla \cdot \Sigma_f + \rho_f \mathbf{g}) - \nabla \cdot \left[ \langle \chi_f \rho_f \mathbf{u}_f' \mathbf{u}_f' \rangle + \langle \delta_f \mu_f (\mathbf{u}_f' \mathbf{n} + \mathbf{n} \mathbf{u}_f') \rangle \right] - \langle \delta_f \sigma_f^{(1)} \cdot \mathbf{n} \rangle, \qquad (5.17) $$

where, $\sigma_f^{(1)} = \sigma_f^0 - \Sigma_f$ which also reads

$$ \sigma_f^{(1)} = -p_f' \delta + \mu_f [\nabla \mathbf{u}_f' + \nabla \mathbf{u}_f'] \qquad (5.18) $$

Likewise, using Eq. (5.14) one may also derive another form of the momentum equation, namely,

$$ \phi_f \rho_f (\partial_t + \mathbf{u}_f \cdot \nabla) \mathbf{u}_f = \phi_f (\nabla \cdot \Sigma + \rho_f \mathbf{g}) - \nabla \cdot \left[ \langle \chi_f \rho_f \mathbf{u}_f' \mathbf{u}_f' \rangle + \langle 2\mu_f \chi_d \mathbf{e}_d^* \rangle \right] - \langle \delta_f \sigma_f^{(2)} \cdot \mathbf{n} \rangle, \qquad (5.19) $$

where $\sigma_f^{(2)} = \sigma_f^0 - \Sigma$. It can be also expressed as,

$$ \sigma_f^{(2)} = -p_f' \delta + \mu_f [\nabla \mathbf{u}_f^* + \nabla \mathbf{u}_f^*] \qquad (5.20) $$

where $\mathbf{u}_f^* = \mathbf{u}_f^0 - \mathbf{u}$. In both cases (Eqs. (5.17) and (5.19)) one has to compute the local stress relative to the averaged continuous phase motion or bulk phase motion. Note the difference between the last two terms in formulations Eqs. (5.17) and (5.19), which is given by $\phi_f \nabla \cdot (\Sigma_f - \Sigma)$. Hence, the transition from one formulation to the other is straightforward and only requires the addition or subtraction of this term. However, it is important to recognize its physical significance. This difference can be expressed explicitly as,

$$ \phi_f \nabla \cdot (\Sigma_f - \Sigma) = \mu_f \phi_f \nabla \cdot [\nabla(\phi_d(\mathbf{u}_f - \mathbf{u}_d)) + \nabla(\phi_d(\mathbf{u}_f - \mathbf{u}_d))] \qquad (5.21) $$

where the decomposition

$$ \mathbf{u} = \phi_f \mathbf{u}_f + \phi_d \mathbf{u}_d, \qquad (5.22) $$

has been used in the derivation. Eq. (5.21) reveals that the difference between the two formulations introduces a non-Newtonian stress term. This non-Newtonian stress shares similarities with the second-order force moment closure (Jackson, 1997; Zhang and Prosperetti, 1997), and it is typically related to intrinsic convection in sedimentation processes (Lhuillier, 2022). Therefore, when addressing the closure problem, it is essential to treat each stress formulation carefully, as the inclusion or exclusion of this additional term can significantly impact the resulting model.

Although the stress decomposition and momentum formulations proposed herein are mathematically equivalent, each presents distinct advantages and limitations, which we examine in detail below. First, we may observe that a simplification arises for Eq. (5.14) in the case of solid particles, for which $\mathbf{e}_d^0 = 0$. As a result $\phi_f \sigma_f = -\phi_f p_f \delta + 2\mu_f \mathbf{E}$ (Joseph et al., 1990; Jackson, 2000). Owing to this simplification, the formulation based on $\Sigma$ and the associated closure relation (Eq. (5.14)) makes it a good candidate to be used in practice for solid particles. Second, employing Eq. (5.17) with Eq. (5.18) may appear more convenient as the governing equations are directly formulated in terms of the fluid velocity field $\mathbf{u}_f$, thereby circumventing the need to explicitly include an expression for the bulk velocity ($\mathbf{u}$). Third, there exists a conceptual reason for preferring the bulk stress $\Sigma$ when deriving closure relations for non-dilute suspensions. Specifically, because $\sigma_f^{(2)}$ depends on the disturbance velocity $\mathbf{u}_f^* = \mathbf{u}_f^0 - \mathbf{u}$, the averaged velocity $\mathbf{u}$ naturally serves as the “far-field” or “undisturbed” velocity boundary condition in the conditionally averaged Navier-Stokes equations (Hinch, 1977; Fintzi, 2025). Notably, most (if not all) existing theoretical solutions addressing the disturbance
field generated by a particle embedded in a Newtonian solvent or in an effective medium — representing other particles contribution — use a divergence-free background velocity field as limiting boundary condition (Hinch, 1977; Kim and Russel, 1985). Since $\nabla \cdot \mathbf{u} = 0$, it follows that the velocity field used in these theoretical study corresponds to $\mathbf{u}$. Consequently, the disturbance stress derived in those studies is $\sigma_f^{(2)} = \sigma_f^0 - \Sigma$. In opposition, the stress decomposition based on $\Sigma_f$, involves the fields $\mathbf{u}_f$ which is not divergence free. Hence, this ‘far field’ or ‘undisturbed’ velocity boundary condition far from the test particle does not correspond to the usual boundary condition assumed in most of the theoretical problems. In conclusion, it is important to recognize that the vector field $\mathbf{u}$ corresponds precisely to the ‘background velocity’ typically employed in theoretical derivations. Consequently, the closure terms obtained in such analyses — such as the hydrodynamic forces, the stresslet, and higher-order moments — are formulated in terms of the difference $\sigma_f^0 - \Sigma$. Of course, it is always possible to express $\Sigma_f$ in terms of $\Sigma$ (see Eq. (5.21)), meaning that the choice of formulation remains free. Nevertheless, we must be cautious when claiming that a specific expression of the momentum exchange term is exactly equivalent to a closure term reported in the literature especially for non-dilute suspension.

Furthermore, given that various stress decomposition have already been introduced in the literature, we propose to examine their respective advantages and limitations. In the seminal work by Zhang and Prosperetti (1997), as well as in numerous subsequent studies, the disturbance stress is defined as $\sigma_f' = \sigma_f^0 - \sigma_f$, a formulation that may be regarded as the “intuitive” definition. However, it is important to note that, by employing Eq. (5.13), one obtains

$$ \sigma_f' = \sigma_f^{(1)} + \frac{\mu_f}{\phi_f} \langle \delta_r(\mathbf{u}_f' \mathbf{n} + \mathbf{n}\mathbf{u}_f') \rangle. \qquad (5.23) $$

The last term represents a contribution from the stress induced by the dispersed phase, as illustrated below. To better understand this term, we consider the case of a suspension composed of rigid solid particles, for which the strain rate tensor of the dispersed phase vanishes, i.e., $\mathbf{e}_d = 0$. By subtracting Eq. (5.13) from Eq. (5.14), we directly obtain

$$ \langle \delta_r(\mathbf{n}\mathbf{u}_f' + \mathbf{u}_f' \mathbf{n}) \rangle = \mathbf{E}_f - \mathbf{E} - \phi_d \mathbf{E}_f = (\mathbf{u}_f - \mathbf{u}_d)\nabla\phi_d + \nabla\phi_d(\mathbf{u}_f - \mathbf{u}_d) - \phi[\nabla\mathbf{u}_d + (\nabla\mathbf{u}_d)^\dagger]. \qquad (5.24) $$

where we have used Eq. (5.22). Therefore, at least in the case of solid particles, this term becomes non-zero whenever there are significant gradients in volume fraction and mean particle velocity. The same observations hold if one considers expression Eq. (5.14) in the analysis above. Since $\sigma_f$ already includes closure contributions associated with the dispersed phase, directly subtracting the full expression of $\sigma_f$ from $\sigma_f^0$ can lead to inconsistencies unless the closure problem is handled with sufficient care. Note that the last term of Eq. (5.23) is proportional to the interface area concentration, hence it is of $O(\phi_d)$. Thus, when integrating $\sigma_f' \cdot \mathbf{n}$ on the interface of the droplets, the contribution from that term to the mean momentum exchange term is of $O(\phi_d^2)$.

Another commonly adopted approach in the literature is $\sigma^{(4)} = \sigma_f^0 - p_f \delta$ (Simonin, 1996; Lhuillier and Nadim, 2009; Morel, 2015; Guazzelli and Pouliquen, 2018). This definition leads to the expression

$$ \sigma_f^{(4)} = -p_f' \delta + \mu_f (\nabla \mathbf{u}_f^0 + \nabla^\dagger \mathbf{u}_f^0), \qquad (5.25) $$

which implies that the closure is based on the absolute local fluid velocity $\mathbf{u}_f^0$, rather than the fluctuating component $\mathbf{u}_f'$. Without going into the details, note that numerous closures are based on the reciprocal theorem formulation (Kim and Karrila, 2013; Stone et al., 2001; Raja et al., 2010). The closures (drag forces, stresslet etc...) provided by the reciprocal theorem are by construction expressed in terms of the disturbance fields ($p_f'$, $\mathbf{u}_f^*$), because they must decay to zero far from the test particle. For example, the Faxen contribution to the drag force in the case of a spherical solid particle of radius $a$ is given by $\mathbf{f} = \pi a^3 \mu_f \nabla^2 \mathbf{u}$. This formulation is obtained by considering the contribution from the disturbance stress, $\sigma^{(2)} = -p_f' \delta + \mu_f (\nabla \mathbf{u}_f^* + \nabla^\dagger \mathbf{u}_f^*)$, which vanish far from the test particle. If one uses the formulation based on $\sigma_f^{(4)} = -p_f' \delta + \mu_f (\nabla \mathbf{u}_f^0 + \nabla^\dagger \mathbf{u}_f^0)$, the Faxen contribution to the drag force becomes $\mathbf{f} = \nu_p \nabla \cdot \Sigma + \pi a^3 \mu_f \nabla^2 \mathbf{u} = \pi a^3 \mu_f (1 + \frac{4}{3}) \nabla^2 \mathbf{u}$ where the second term is the contribution from the mean stress. Although both formulations are equally valid, the second one appears to be less commonly used and could potentially lead to misinterpretations or inconsistencies if not carefully handled.

In light of these considerations, we adopt the formulation based on $\sigma^{(2)}$-i.e., Eq. (5.19)-for the remainder of this work, as it is less prone to errors when considering the closure problem. For simplicity, we will denote $\sigma^{(2)}$ as $\sigma^*$ throughout the rest of the paper. The last two terms on the right-hand side of Eq. (5.19) can be further expanded into a Taylor series using Eq. (4.5). Doing so leads us to the hybrid formulation of the continuous phase momentum equation,

$$ \phi_f \rho_f (\partial_t + \mathbf{u}_f \cdot \nabla) \mathbf{u}_f = \phi_f (\nabla \cdot \Sigma + \rho_f \mathbf{g}) + \nabla \cdot \sigma_f^{\text{eff}} - \left\langle \delta_p \oint_{\Gamma_a} \sigma_f^* \cdot \mathbf{n} \, d\Gamma \right\rangle, \qquad (5.26) $$

where,

$$ \sigma^{\text{eff}} = - \left\langle \chi_f \rho_f \mathbf{u}_f' \mathbf{u}_f' \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} \sigma_f^* \cdot \mathbf{n} \, d\Gamma - \delta_p \int_{\Omega_a} 2\mu_f e_d^* d\Omega \right\rangle - \nabla \cdot \left\langle \delta_p \frac{1}{2} \oint_{\Gamma_a} \mathbf{r} \mathbf{r} \sigma_f^* \cdot \mathbf{n} \, d\Gamma - \delta_p \int_{\Omega_a} 2\mu_f \mathbf{r} e_d^* d\Omega \right\rangle + \nabla \cdot \mathbf{v}(\cdot). \qquad (5.27) $$

Under this form, the left-hand side of Eq. (5.26) represents the total derivative of $\mathbf{u}_f$, while on the right-hand side we find: (1) the mean Newtonian stress contribution based on the bulk velocity $\Sigma$, (2) the mean buoyancy force, (3) the Reynolds stress term, and (4) the moments of momentum exchange terms, which are computed based on $\sigma_f^*$. The (...) refers to the higher order moments.

The averaged mass of droplets ($m_p$), the averaged center of mass velocity ($\mathbf{u}_p$), the averaged second moment of mass ($\mathbf{M}_p$), and the averaged first moment of momentum ($\mathbf{P}_p$), obey conservation laws that are given according to Eqs. (4.13), (4.14) and (C.7) (for conservation laws at the local scale, refer to the previous section). They read,

$$ (\partial_t + \mathbf{u}_p \cdot \nabla) n_p = -n_p \nabla \cdot \mathbf{u}_p \qquad (5.28) $$

$$ n_p (\partial_t + \mathbf{u}_p \cdot \nabla) \mathbf{M}_p + \nabla \cdot \langle \delta_p \mathbf{u}_p' \mathbf{M}_p \rangle = n_p 2 (\mathbf{S}_p + \mathbf{P}_p \delta) \qquad (5.29) $$

$$ m_p n_p (\partial_t + \mathbf{u}_p \cdot \nabla) \mathbf{u}_p + \nabla \cdot \langle \delta_p m_p \mathbf{u}_p' \mathbf{u}_p' \rangle = m_p n_p \mathbf{g} + \left\langle \delta_p \oint_{\Gamma_a} \sigma_f^* \cdot \mathbf{n} \, d\Gamma \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} \Sigma \cdot \mathbf{n} \, d\Gamma \right\rangle \qquad (5.30) $$
<table><caption>**Table 2**

Definition of physical quantities and dimensionless parameters. Note that the choice of velocity scales depends on the specific problem under consideration. For example, in the case of sedimenting droplets in a quiescent fluid, the characteristic velocity is typically $U \sim |\mathbf{u}_p|$.</caption><tr><td>Velocity scale</td><td>U</td></tr><tr><td>Macroscopic length scale</td><td>L</td></tr><tr><td>Droplets radius</td><td>a</td></tr><tr><td>Reynolds number</td><td>Re = ρ<sub>f</sub> aU / μ<sub>f</sub></td></tr><tr><td>Capillary number</td><td>Ca = μ<sub>f</sub> U / γ</td></tr><tr><td>Marangoni number</td><td>Ma = a|∇γ| / μ<sub>f</sub> U</td></tr><tr><td>Viscosity ratio</td><td>λ = μ<sub>d</sub> / μ<sub>f</sub></td></tr><tr><td>Density ratio</td><td>ζ = ρ<sub>d</sub> / ρ<sub>f</sub></td></tr></table>

$$n_p(\partial_t + \mathbf{u}_p \cdot \nabla)\mu_p + \nabla \cdot \langle \delta_p \mathbf{u}_a' \mu_a \rangle = \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} \times (\boldsymbol{\sigma}_f^* \cdot \mathbf{n}) d\Gamma \right\rangle \quad (5.31)$$

$$n_p(\partial_t + \mathbf{u}_p \cdot \nabla)S_p + \nabla \cdot \langle \delta_p \mathbf{u}_a' S_a \rangle = \rho_d \left\langle \delta_p \int_{\Omega_a} \mathbf{w}_d^0 \mathbf{w}_d^0 - \frac{1}{3} (\mathbf{w}_d^0 \cdot \mathbf{w}_d^0) \delta d\Omega \right\rangle - \left\langle \delta_p \int_{\Omega_a} 2\mu_d \mathbf{e}_d^* d\Omega \right\rangle \\ + \left\langle \delta_p \oint_{\Gamma_a} \frac{1}{2} (\mathbf{r}\boldsymbol{\sigma}_f^* + \gamma \mathbf{r}\boldsymbol{\sigma}_f^* - \frac{2}{3}(\boldsymbol{\sigma}_f^* \cdot \mathbf{r})\delta) \cdot \mathbf{n} d\Gamma \right\rangle \\ - \left\langle \delta_p \oint_{\Gamma_a} \gamma \left(\frac{1}{3}\delta - \mathbf{n}\mathbf{n}\right) d\Gamma \right\rangle + (1-\lambda) \left\langle \delta_p \int_{\Omega_a} 2\mu_f \mathbf{E} d\Omega \right\rangle + \frac{1}{2} \left\langle \delta_p \int_{\Omega_a} \mathbf{r}(\nabla \cdot \boldsymbol{\Sigma}) + (\nabla \cdot \boldsymbol{\Sigma})\mathbf{r} d\Omega \right\rangle, \quad (5.32)$$

The presented system of equations extends the averaged equations developed for non-spherical solid particles (Curtiss, 1956). It is worth noting that, in the case of solid particles, Eq. (5.32) becomes redundant, as the tensor $S_p$ can be expressed in terms of $M_p$, $\mu_p$ or the mean angular velocity, and the correlations of their fluctuations. In this case, Eq. (5.32) must be solved to determine the unknown internal stress of the solid particles.

The set of Eqs. (5.8), (5.26), (5.28)–(5.32) is completed by the following relations

$$\phi_f + \phi_d = \phi_f + \phi + \frac{1}{2} \nabla \nabla : (M_p n_p) + \dots = 1, \quad (5.33)$$

$$\mathbf{u} = \mathbf{u}_f \phi_f + \phi \mathbf{u}_p - \frac{1}{\rho_d} \nabla \cdot (\mathbf{P}_p n_p) + \dots \quad (5.34)$$

where we have introduced the notation $\phi = n_p v_p$. Eq. (5.34) is obtained by inserting expansion Eq. (4.4) in Eq. (5.22).

### 5.3. Symmetry of the effective stress tensor

We remark that the second moment and higher-order moments in Eq. (5.27) appear under two divergence operators in Eq. (5.26). Hence, if we note $\Sigma_{ijk}$ the third rank tensor that represent these moments, then only the vector $\partial_k \partial_j \Sigma_{ijk}$ is of physical significance in the momentum balance (Eq. (5.26)). Thus, one can demonstrate that (Lhuillier, 1996)

$$\partial_j \partial_k \Sigma_{ijk} = \partial_j \partial_k \Sigma_{i(jk)} = \partial_j \partial_k [\Sigma_{i(jk)} + \Sigma_{j(ik)} - \Sigma_{k(ij)}], \quad (5.35)$$

where $\Sigma_{i(jk)} = \frac{1}{2}[\Sigma_{ijk} + \Sigma_{ikj}]$ represents the symmetric part of $\Sigma_{ijk}$ over the index $jk$, as indicated by the parenthesis (and so on for the other tensor). This expression is allowed because $\partial_j \partial_k (\Sigma_{ijk} - \Sigma_{ikj}) = 0$ and $\partial_j \partial_k (\Sigma_{j(ik)} - \Sigma_{k(ij)}) = 0$. This manipulation highlights the fact that the effective stress due to the second order moments remains symmetric over the indices $ij$, in all circumstances. Hence, as already demonstrated by Lhuillier (1996) only the hydrodynamic torque can induce skew-symmetric stresses in the effective stress of the averaged momentum equation.

### 6. Closure for dilute suspensions of droplets in viscous dominated flows

We now consider the closures for a dilute, monodisperse suspension of spherical droplets with radius $a$ in Stokes flow. Our analysis recovers the results of Zhang and Prosperetti (1997, Appendix B). We extend these results by considering the influence of Marangoni effects on the suspension dynamics. Additionally, we demonstrate how higher-order moments can characterize droplet shapes and relate them to the averaged equations governing the continuous phase. We also address the closure of the various covariance terms arising in the averaged equations.

The closure terms in the equations above are expressed in terms of $p'_f$ and $\mathbf{u}_f^*$. Our objective is to determine the disturbance pressure and velocity fields generated by a spherical droplet immersed in an arbitrary flow. Specifically, the pair $(\mathbf{u}_f^*, p'_f)$ represents the solution of the single-particle conditionally averaged Navier-Stokes equations (Hinch, 1977; Zhang and Prosperetti, 1994a; Fintzi, 2025). When deriving averaged equations accurate to $O(\phi)$ we can neglect of droplet-droplet interactions (Hinch, 1977; Zhang and Prosperetti, 1994a). Furthermore, we assume negligible inertia in the closure problem by omitting all terms of $O(Re)$ in the closure equations, and correspondingly, all terms of $O(Re\phi)$ in the averaged equations (a summary of the dimensionless parameters is provided in Table 2). Although surface tension gradients are included, we assume spherical droplet shapes. This allows us to neglect all terms of $O(Ca)$ in the closure problem. As a result, we focus on the motion of an isolated spherical droplet in an arbitrary Stokes flow, subject to tangential stress discontinuities at the interface. The analytical solution to this problem is well established and available in various sources, including Subramanian (1985), Nadim and Stone (1991), Pozrikidis et al. (1992), Leal (2007), Raja et al. (2010), Pozrikidis (2011) and Kim and Karrila (2013), and forms the basis for evaluating the closure terms. For clarity, we detail the closure problem and its solution in Appendix F.

In the following, $L$ denotes the characteristic length scale over which the averaged quantities may vary (Jackson, 1997). Each averaged property may possess a different variation length scale. For instance, consider $\nabla u$, which typically varies at the length scale of the process, whereas $\nabla \gamma$ may
vary over a smaller length scale. For example, if surface tension gradient is generated by surfactant transport, $\gamma$ typically varies over a length scale of the drop size. If $\gamma$ varies due to a non-constant temperature field, the length scale of variation will correspond to that of the temperature field, which is related to the process or macro scale. To simplify the problem, we assume that $\nabla\gamma$, and $\nabla u$ typically vary on the order of ~ L. We then consider a proper separation of scales such that $a \ll L$ and therefore all contributions proportional to $O(a^2/L^2)$ become negligible in the averaged equations (Jackson, 1997; Zhang and Prosperetti, 1997).

## 6.1. Hydrodynamic stresses closures

In the first place we focus on the momentum exchange terms. We may directly compute the following expressions from the singularity solutions and find,<fnref n="3" />

$$ \left\langle \delta_p \oint_{\Gamma_a} \boldsymbol{\sigma}_f^* \cdot \mathbf{n} \, d\Gamma \right\rangle = \phi \frac{\mu_f}{a^2} \frac{3(2+3\lambda)}{2(1+\lambda)} \mathbf{u}_r + \phi \mu_f \frac{3\lambda}{4(\lambda+1)} \nabla^2 \mathbf{u} + \phi \frac{1}{a} \frac{1}{\lambda+1} \nabla\gamma + \phi a \frac{1}{10(\lambda+1)} \nabla^2(\nabla\gamma), \quad (6.1) $$

$$ \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r}\boldsymbol{\sigma}_f^* \cdot \mathbf{n} \, d\Gamma \right\rangle = \mu_f \phi \frac{3(5\lambda + 2)}{5(\lambda + 1)} \mathbf{E} + \phi a \frac{9}{25(\lambda + 1)} (\nabla\nabla\gamma - \delta\nabla^2\gamma/3), \quad (6.2) $$

$$ \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r}\boldsymbol{\sigma}_f^* \cdot \mathbf{n} \, d\Gamma \right\rangle = \mu_f \phi \frac{3}{5(\lambda+1)} (\mathbf{u}_r \delta + \delta^+ \mathbf{u}_r \delta) + \mu_f \phi \frac{3(5\lambda+2)}{10(\lambda+1)} \delta \mathbf{u}_r - \phi a \frac{2}{5(\lambda+1)} (\nabla\gamma \delta + \delta^+ \nabla\gamma \delta) + \phi a \frac{3}{5(\lambda+1)} \delta \nabla\gamma, \quad (6.3) $$

$$ \left\langle \delta_p \oint_{\Gamma_a} 2\mu_f \mathbf{e}_d^* d\Gamma \right\rangle = -\mu_f \phi \frac{2(5\lambda + 2)}{5(\lambda + 1)} \mathbf{E} - \phi a \frac{6}{25(\lambda + 1)} (\nabla\nabla\gamma - \delta\nabla^2\gamma/3), \quad (6.4) $$

$$ \left\langle \delta_p \oint_{\Gamma_a} 2\mu_f \mathbf{r}\boldsymbol{\mathbf{e}}_d^* d\Gamma \right\rangle = \phi \frac{3\pi}{10(\lambda+1)} (\delta\mathbf{u}_r + \delta\mathbf{u}_r^\dagger) - \phi \frac{\pi}{5(\lambda+1)} \mathbf{u}_r \delta - \phi a \frac{1}{5(\lambda+1)} (\delta\nabla\gamma + \delta\nabla\gamma^\dagger) + \phi a \frac{2}{15(\lambda+1)} \nabla\gamma\delta, \quad (6.5) $$

where we have introduced the relative velocity $\mathbf{u}_r = \mathbf{u} - \mathbf{u}_p$, and the viscosity ratio $\lambda = \mu_d/\mu_f$. Most of the terms proportional to $\mathbf{E}$ and $\mathbf{u}_r$ in this expression are already well known and will not be discussed here. In Eq. (6.1) the third term represents the force induced by $\nabla\gamma$, which typically represents the thermocapillary migration of droplets if one relates $\nabla\gamma$ to a temperature gradient. The last term of Eq. (6.1) corresponds to a “Faxon-like” contribution of the Marangoni forces and involves the third derivative $\gamma$. It is worth noting that the first Marangoni moments of hydrodynamic forces scale as $\propto \phi\delta\nabla\gamma$ (see Eq. (6.2)), while the second moments scale as $\propto \phi\delta\nabla\gamma$ (Eq. (6.3)). These tensorial form could have been anticipated based on symmetry considerations. Note that we have neglected the terms proportional to $\nabla^2\mathbf{E}$ in the first moment expression, and to $\nabla\nabla\mathbf{u}$ or $\nabla\nabla\nabla\gamma$ in the second moment expression, as they are of order $O(a^2/L^2)$ or smaller. The complete expressions for the second moments, including these higher-order contributions, are provided in Appendix F.

The above closure provides the contribution from the disturbance fields, however in the dispersed phase Eqs. (5.30) and (5.32) one need the contribution from the total momentum exchange including the contribution of the mean stress $\boldsymbol{\Sigma}$. These terms can be obtained following the procedure outlined in Zhang and Prosperetti (1997) and Morel (2015), and reads,

$$ \left\langle \delta_p \oint_{\Gamma_a} \boldsymbol{\Sigma} \cdot \mathbf{n} \, d\Gamma \right\rangle = \phi \nabla \cdot \boldsymbol{\Sigma} + O(a^2/L^2), \quad (6.6) $$

$$ \left\langle \delta_p \int_{\Omega_a} \mathbf{E} d\Omega \right\rangle = \phi \mathbf{E} + O(a^2/L^2), \quad (6.7) $$

$$ \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r}(\nabla \cdot \boldsymbol{\Sigma}) d\Gamma \right\rangle = O(a^2/L^2). \quad (6.8) $$

Where we have used the approximation of $\boldsymbol{\Sigma}_f(\mathbf{x} + \mathbf{r}) = \boldsymbol{\Sigma}_f(\mathbf{x}) + \mathbf{r} \cdot \nabla \boldsymbol{\Sigma}_f|_{\mathbf{r}=0} + \dots$ and neglected the $O(a^2/L^2)$ terms.

## 6.2. Velocity variance and covariance closures

The disturbance velocity field $\mathbf{u}'_f$ is proportional to $\propto \mathbf{u}_r$, $\nabla\mathbf{E}_f$ and $\nabla\nabla\mathbf{u}_f$ depending on the problem at hand. Additionally, the Reynolds stress tensor $\langle \chi_f \mathbf{u}'_f \mathbf{u}'_f \rangle$ is a symmetric second-order tensor. We deduce that the functional form of the Reynolds stress must be

$$ \langle \chi_f \rho_f \mathbf{u}'_f \mathbf{u}'_f \rangle = C_{uu}^1(\phi, \lambda) \rho_f \mathbf{u}_r \mathbf{u}_r + C_{uu}^2(\phi, \lambda) \rho_f (\mathbf{u}_r \cdot \mathbf{u}_r) \delta + a^2 C_{EE}^1(\phi, \lambda) \rho_f \mathbf{E} \cdot \mathbf{E} + a^2 C_{EE}^2(\phi, \lambda) \rho_f (\mathbf{E} : \mathbf{E}) \delta + \dots \quad (6.9) $$

where the remaining terms indicated by the ... represent linear combination of terms proportional to $a^4\nabla\nabla\mathbf{u}_f : \nabla\nabla\mathbf{u}_f$. However, since these terms are proportional to $a^2$, they scale as $O(a^2/L^2)$ in the averaged equations. As a result, only the first two terms in Eq. (6.9) are significant in the momentum equation.<fnref n="4" /> The exact values of the coefficients $C_{EE}$ are provided in Raja et al. (2010). Because the Reynolds stress is defined as an average over the continuous-phase domain (denoted by $\chi_f$), the disturbance velocity fields $\mathbf{u}'_f$ cannot be directly integrated to determine the constants $C_{uu}^1$ and $C_{uu}^2$. Nonetheless, based on experimental measurements by Cartellier et al. (2009), particle-resolved simulations by Fintzi (2025), and theoretical results in triply periodic domains by Hill et al. (2001), it is reasonable to expect that these constants follow the scaling:

$$ C_{uu}^1, C_{uu}^2 \propto \phi^{2/3} \frac{(2+3\lambda)^2}{(\lambda+1)^2}. \quad (6.10) $$

<footnote n="3">³ We follow the convention of Happel and Brenner (2012) for the transpose operator. For an arbitrary order tensor $\mathbf{A}$, its transpose is denoted either as $\dagger A_{ijkl...} = A_{ijkl...}$ or as $(A..._{ijkl...})^\dagger = A..._{ijkl...}$.</footnote>

<footnote n="4">⁴ Note that since $\mathbf{u}'_f \propto \frac{\delta\nabla\gamma}{\mu_f}$ it follows that $\mathbf{u}'_f \mathbf{u}'_f \propto O(a^2/L^2)$. Thus, the contribution from $\nabla\gamma$ ends up being negligible under the current modeling assumptions.</footnote>
In the present study, we have neglected droplet-droplet interactions in the closure problem, and therefore expect $\langle \delta_p u'_a u'_a \rangle = 0$. Nonetheless, by invoking symmetry arguments and drawing on experimental observations reported in the literature (Guazzelli and Hinch, 2011), we conclude that:

$$ \langle \delta_p m_p u'_a u'_a \rangle = \rho_d C_{up}^1 (\phi, \lambda) \mathbf{u}_r \mathbf{u}_r + \rho_d C_{up}^2 (\phi, \lambda) \delta(\mathbf{u}_r \cdot \mathbf{u}_r) \quad (6.11) $$

where $C_{up}^1$ and $C_{up}^2$ are unknown constants which are $\propto \phi^{2/3}$ (Guazzelli and Hinch, 2011). This result indicates that, due to the long-range interactions between droplets, the velocity variance of the particles remains non-zero even at order $O(\phi)$. It is important to note that both $\langle \delta_p m_p u'_a u'_a \rangle$ and $\langle \chi_f \rho_f u'_f u'_f \rangle$ represent inertial contributions by construction. However, when expressed in dimensionless form, both terms scale as $O(\text{Re}\phi^{2/3})$.

Given that $O(\phi^{2/3} \text{Re}) \gg O(\text{Re}\phi)$ in the dilute limit, we conclude that these velocity variance terms must be retained to achieve closure at order $O(\text{Re}\phi)$.

The covariance terms appearing on the left-hand side of Eqs. (5.29) to (5.32), reflect the correlation between the shape of the droplet ($M_\alpha$), its angular momentum ($\mu_\alpha$), and its stretching of momentum ($S_\alpha$), with its center of mass velocity $\mathbf{u}_\alpha$. If any of these quantities is statistically independent of $\mathbf{u}_\alpha$, the corresponding covariance term vanishes. In the dilute Stokes regime, a purely translating spherical droplet remains undeformed due to the balance of normal stresses at its surface (Leal, 2007). Similarly, translation of a spherical particles does not induce hydrodynamic torque, and thus no angular momentum is generated. Therefore, in the Stokes regime, the quantities $M_\alpha$, $P_\alpha$ are uncorrelated with $\mathbf{u}_\alpha$, implying that the covariance terms $\langle \delta_p u'_a M'_a \rangle$, $\langle \delta_p u'_a \mu'_a \rangle$ and $\langle \delta_p u'_a S'_a \rangle$ equal zero. However, these conclusions no longer hold at finite inertia. In that case, the translational-rotational coupling (Rubinow and Keller, 1961) leads to nonzero force on the particle, and the droplet can deform as a result of its motion relative to the surrounding fluid (Taylor and Acrivos, 1964).

### 6.3. Closed form of the hybrid model

Note that to leading order in the expansion in $Ca$ and $Re$, the linear momentum equations are decoupled from the dispersed phase moments ($S_p$, $\mu_p$, and $M_p$). As a result, Eqs. (5.29), (5.31) and (5.32) are not required to compute the mean variables ($\mathbf{u}_p$, $\mathbf{u}$, $\rho_f$, $\phi$). Therefore, one can directly substitute expressions Eqs. (6.1)–(6.8) into Eqs. (5.26) and (5.30) to obtain a closed form of the hybrid model which read,

$$ \phi_f + \phi = 1 \quad (6.12) $$

$$ \phi_f \mathbf{u}_f + \phi \mathbf{u}_p = \mathbf{u} \quad (6.13) $$

$$ \nabla \cdot \mathbf{u} = 0 \quad (6.14) $$

$$ (\partial_t + \mathbf{u} \cdot \nabla) \phi = \nabla \cdot (\phi \mathbf{u}_r) \quad (6.15) $$

$$ \rho_d \phi (\partial_t + \mathbf{u}_p \cdot \nabla) \mathbf{u}_p = \phi (\nabla \cdot \mathbf{\Sigma} + \rho_d \mathbf{g}) + \nabla \cdot \boldsymbol{\sigma}_p^{\text{eff}} + \mathbf{F} \quad (6.16) $$

$$ \phi_f \rho_f (\partial_t + \mathbf{u}_f \cdot \nabla) \mathbf{u}_f = \phi_f (\nabla \cdot \mathbf{\Sigma} + \rho_f \mathbf{g}) + \nabla \cdot \boldsymbol{\sigma}_f^{\text{eff}} - \mathbf{F} \quad (6.17) $$

where the closures are given by,

$$ \mathbf{F} = \phi \frac{\mu_f}{a^2} \frac{3(2 + 3\lambda)}{2(1 + \lambda)} \mathbf{u}_r + \phi \mu_f \frac{3\lambda}{4(\lambda + 1)} \nabla^2 \mathbf{u} + \phi \frac{1}{a} \frac{1}{\lambda + 1} \nabla \gamma + \phi a \frac{1}{10(\lambda + 1)} \nabla^2 (\nabla \gamma) \quad (6.18) $$

$$ \begin{aligned} \boldsymbol{\sigma}_f^{\text{eff}} = & \mu_f \phi \frac{5\lambda + 2}{2(\lambda + 1)} \mathbf{E} - \mu_f \frac{3\lambda}{4(\lambda + 1)} [\nabla(\phi \mathbf{u}_r) + \nabla(\phi \mathbf{u}_r)^\dagger] + \mu_f \frac{3\lambda - 2}{4(\lambda + 1)} \nabla \cdot (\phi \mathbf{u}_r)\delta \\ & + a\phi \frac{3}{5(\lambda + 1)} (\nabla \nabla \gamma - \delta \nabla^2 \gamma / 3) - a \frac{1}{(\lambda + 1)} \nabla (\phi \nabla \gamma) - a \frac{5}{6(\lambda + 1)} \nabla \cdot (\phi \nabla \gamma) \\ & - \rho_f C_{uu}^1 (\phi, \lambda) \mathbf{u}_r \mathbf{u}_r - \rho_f C_{uu}^2 (\phi, \lambda) (\mathbf{u}_r \cdot \mathbf{u}_r)\delta \end{aligned} \quad (6.19) $$

$$ \boldsymbol{\sigma}_p^{\text{eff}} = - \rho_d C_{up}^1 (\phi, \lambda) \mathbf{u}_r \mathbf{u}_r - \rho_d C_{up}^2 (\phi, \lambda) (\mathbf{u}_r \cdot \mathbf{u}_r) \delta \quad (6.20) $$

This system is constituted of 6 unknowns ($\phi_f$, $\phi$, $\mathbf{u}$, $\mathbf{u}_f$, $\rho_f$, $\mathbf{u}_p$), and 6 equations (Eqs. (6.12) to (6.17)). Upon the precise knowledge of $C_{up}^1$, $C_{up}^2$, $C_{uu}^1$ and $C_{uu}^2$ the system can be considered closed and accurate to order $O(\phi)$.

First note that most of the terms related to the phase relative motion (i.e. those $\propto \mathbf{u}_r$ or $\mathbf{E}$) are already presented in Zhang and Prosperetti (1997, Appendix B). As already noted in several studies $\nabla\gamma$ contributes to the hydrodynamic forces on the droplets (Subramanian, 1985), hence also contributes as a source term in the equation of $\mathbf{u}_f$. What has been unnoticed before, however, is that terms of the form $\sim \nabla\nabla\gamma$ also appear to induce non-newtonian stresses within the suspension. Although the problem considered in this work differs from that typically encountered in industrial processes –here, a fixed surface tension gradient is imposed rather than solving a scalar transport equation that governs surface tension locally– it is still insightful to understand under which conditions surface tension gradients (i.e., Marangoni effects) become significant. First, it is worth noting that Marangoni effects become negligible when the viscosity ratio is high, meaning the dispersed phase (e.g., droplets) is much more viscous than the surrounding fluid. Consequently, one expects Marangoni flows to be more prominent for bubbles than for viscous droplets. Consider a single bubble rising in a quiescent fluid under the influence of gravity and an imposed uniform surface tension gradient, or equivalently, a constant concentration gradient of a surfactant species at infinity, with a linear dependence of surface tension on concentration. Balancing the Marangoni force with the buoyant force gives the scaling, $\nabla\gamma \sim (\rho_d - \rho)ga$. Assuming the surface tension gradient scales as $\nabla\gamma \sim \Delta\gamma/L$, where $\Delta\gamma$ is the typical variation in surface tension over a characteristic length $L$, we can estimate when Marangoni effects become significant. For instance, in an aqueous emulsion (with $\rho \approx 1000 \text{ kg/m}^3$, $\mu \approx 10^{-3} \text{ Pa} \cdot \text{s}$, and $\rho_d - \rho \approx 100 \text{ kg/m}^3$), and assuming $\Delta\gamma \approx 10^{-2} \text{ N/m}$ and $a \approx 10 \text{ µm}$, Marangoni and buoyancy forces become comparable when $L \approx 1 \text{ m}$. This suggests that even when surface tension varies over meter-scale distances, Marangoni forces can significantly influence the dynamics of very small bubbles or low viscosity ratio droplets. Next, we consider a droplet subjected to a shear flow with an imposed second-order spatial variation of surface tension. Balancing the effective shear stress due to the flow with the stress induced by Marangoni effects yields the scaling $\nabla^2\gamma \sim \frac{\mu\dot{\gamma}}{a}$, where $\dot{\gamma}$ is the shear rate. For a typical shear rate of $\dot{\gamma} \approx 1 \text{ s}^{-1}$ and droplet size $a \approx 10 \text{ µm}$, the characteristic length scale over which surface tension must vary for Marangoni stresses to match shear-induced stresses is approximately $L \sim 1 \text{ cm}$. Therefore, in processes where surface tension varies over centimeter-scale distances, Marangoni effects can generate effective stresses comparable in
magnitude to those from shear-induced viscosity. This highlights the importance of accounting for surface tension gradients in modeling interfacial flows, particularly in systems involving small bubbles or droplets.

For completeness the averaged mass and second moment of mass read as,

$$n_p m_p = \rho_d \phi \qquad (6.21)$$

$$n_p M_p = \rho_d \phi \frac{a^2}{5} \delta \qquad (6.22)$$

Substituting these closures into Eqs. (5.31) and (5.32) yields

$$S_p = 0 \qquad (6.23)$$

$$n_p(\partial_t + \mathbf{u}_p \cdot \nabla)\mu_p = 0, \qquad (6.24)$$

as no hydrodynamic torque is exerted on the droplet.

## 6.4. Small droplet deformations

The usual procedure to determine the droplet deformation is to use the normal stress balance (Nadim and Stone, 1991; Nadim, 1996). However, to highlight the significance of the first-moment of momentum equation, we propose an alternative method based on Eqs. (5.29), (5.31) and (5.32). As we now demonstrate, with the closure relations derived in the previous section, Eq. (5.32) can be used to compute the first-order correction to the droplet spherical shape. In other words, using closure terms accurate to O(Ca), this approach yields the leading order contribution of Ca to droplet deformation.

Let us describe points lying in the droplet $\alpha$ using the parametric equation in the local spherical reference frame $(r, \theta, \phi)$,

$$\mathbf{r}(r, \phi, \theta) = r[1 + C a f_{\alpha}(\phi, \theta)]\mathbf{e}, \qquad (6.25)$$

where $0 < r < a$ is the radial parameter, $\theta$ the polar angle, and $\phi$ the azimuthal angle. $f_{\alpha}(\theta, \phi)$ is the shape function of the droplet, and $\mathbf{e} = \cos \phi \sin \theta \mathbf{e}_x + \sin \phi \sin \theta \mathbf{e}_y + \cos \theta \mathbf{e}_z$ is the radial unit vector. Using the parametrization given by Eq. (6.25) one can eventually compute the volume $d\Omega$ and surface $d\Gamma$ element in terms of $r, \phi, \theta$ and the deformation function $f_{\alpha}(\theta, \phi)$. The result are, $d\Gamma = a^2(1+2C a f_{\alpha}(\theta, \phi)) \sin \theta d\theta d\phi$ and $d\Omega = (1+3C a f_{\alpha}(\theta, \phi))r^2 \sin \theta dr d\theta d\phi$ at the leading order in Ca. Following, Nadim (1996) and Nadim and Stone (1991) we then expand $f_{\alpha}(\phi, \theta)$ in a series of surface harmonics centered at the droplet center of mass, namely

$$f_{\alpha}(\mathbf{e}) = \sum_{n=2}^{\infty} S^{(n)} : \mathbf{H}_{\alpha}^{(n)}, \qquad (6.26)$$

with $S^{(n)} = \frac{(-1)^{n}r^{n+1}}{1 \dots (2n-1)} \nabla^{(n)}(\mathbf{r})$ the nth order surface spherical harmonic, and $\mathbf{H}_{\alpha}^{(n)}$ an nth order symmetric and traceless tensor to be determined. Inserting Eq. (6.25) in the definition of the second moment of mass yields,

$$n_p M_p = \rho_d \phi \frac{a^2}{5} (\delta + 2Ca\mathbf{H}_p^{(2)}), \qquad (6.27)$$

where we used the definition $n_p \mathbf{H}_p^{(n)} = \langle \delta_p \mathbf{H}_\alpha^{(n)} \rangle$. To obtain, closed-form expression for the tensors $\mathbf{H}_p^{(2)}$, we consider the first moment of the forces, as given by Eq. (5.32). At leading order in Ca, Re and $a^2/L^2$, we have $S_p = 0$ and the first and last term on the left-hand side of Eq. (5.32) are negligible. Then Eq. (5.32) yields,

$$\left\langle \delta_p \oint_{\Gamma_a} \gamma \left( \frac{1}{3} \delta - \mathbf{n}\mathbf{n} \right) d\Gamma \right\rangle = - \left\langle \delta_p \int_{\Omega_a} 2\mu_d \mathbf{e}_d^* d\Omega \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} \frac{1}{2} (\mathbf{r}\sigma_f^* + {}^\dagger \mathbf{r}\sigma_f^* - \frac{2}{3}(\sigma_f^* \cdot \mathbf{r})\delta) \cdot \mathbf{n} d\Gamma \right\rangle + (1-\lambda) \left\langle \delta_p \int_{\Omega_a} 2\mu_f \mathbf{E} d\Omega \right\rangle$$

All the terms on the right-hand side of this equation are provided by the closure introduced in Section 6.1. However to close this equation, we must express the quantity $\langle \delta_p \oint_{\Gamma_a} \gamma(\delta/3 - \mathbf{n}\mathbf{n}) d\Gamma \rangle$ as a function of $\mathbf{H}_p^{(n)}$. We first Taylor expand the surface tension about the point $\mathbf{x}$ as $\gamma|_{\mathbf{x}_\alpha} = \gamma|_{\mathbf{x}} + \mathbf{r} \cdot \nabla\gamma|_{\mathbf{x}} + \frac{1}{2}\mathbf{r}\mathbf{r} : \nabla\nabla\gamma|_{\mathbf{x}} + \dots$ to get

$$\left\langle \delta_p \oint_{\Gamma_a} \gamma \left( \frac{\delta}{3} - \mathbf{n}\mathbf{n} \right) d\Gamma \right\rangle = \gamma \left\langle \delta_p \oint_{\Gamma_a} \left( \frac{\delta}{3} - \mathbf{n}\mathbf{n} \right) d\Gamma \right\rangle + \nabla\gamma \cdot \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} \left( \frac{\delta}{3} - \mathbf{n}\mathbf{n} \right) d\Gamma \right\rangle + \dots \qquad (6.28)$$

Then, by injecting the approximation $\mathbf{n} = \mathbf{e} - Ca(\delta - \mathbf{e}\mathbf{e}) \cdot \nabla f_{\alpha} + O(Ca^2)$ into the expression of $f_{\alpha}$ given by Eq. (6.26) yields after some algebra

$$\left\langle \delta_p \oint_{\Gamma_a} \left( \frac{\delta}{3} - \mathbf{n}\mathbf{n} \right) d\Gamma \right\rangle = \frac{1}{a} \frac{8}{5} C a \phi \mathbf{H}_p^{(2)} \qquad (6.29)$$

$$\left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} \left( \frac{\delta}{3} - \mathbf{n}\mathbf{n} \right) d\Gamma \right\rangle = \frac{6}{7} C a \phi \mathbf{H}_p^{(3)} \qquad (6.30)$$

$$\left\langle \delta_p \oint_{\Gamma_a} \left( \frac{\delta}{3} - \mathbf{n}\mathbf{n} \right)_{ij} (\mathbf{r}\mathbf{r})_{kl} d\Gamma \right\rangle = a \frac{4}{15} C a \phi (\mathbf{H}_p^{(2)} \delta)_{ijkl}^p + a C a \phi \frac{16}{35} \mathbf{H}_{ijkl}^{(4)} \qquad (6.31)$$

In the last equation we have introduced the shorthand,

$$(\mathbf{H}_p^{(2)} \delta)_{ijkl}^p = \{ H_{ij} \delta_{kl} + \frac{2}{9} H_{kl} \delta_{ij} - \frac{1}{6} (H_{ik} \delta_{jl} + H_{jl} \delta_{ik} + H_{il} \delta_{jk} + H_{jk} \delta_{il}) \}. \qquad (6.32)$$

Injecting Eqs. (6.29) to (6.31) and the closures derived in Section 6.1 into Eq. (5.32) gives directly,

$$\mathbf{H}_p^{(2)} + \frac{30}{56} \frac{a}{\gamma} \nabla\gamma \cdot \mathbf{H}_p^{(3)} + \dots = \frac{19\lambda + 16}{8(\lambda + 1)} \left( \frac{a\mu_f}{\gamma Ca} \right) E + \frac{3}{40} \frac{(2\lambda + 3)}{(\lambda + 1)} \frac{a^2}{Ca\gamma} (\nabla\gamma - \nabla^2\gamma\delta/3) \qquad (6.33)$$
Discarding first the surface tension gradients in Eq. (6.33), we recover the classical results of Taylor (1932) (see also Rallison, 1984) for droplet deformation in a pure linear flow. As a result, the mean deformation tensor $H_p^{(2)}$ is found to scale with the mean rate-of-strain tensor $E$, non-dimensionalized by the shear rate scale $a\mu_f/(Ca\gamma) = a/U$.<fnref n="5" /> Gradients in surface tension introduce considerable complexity into Eq. (6.33). The first consequence of the surface tension gradient is the appearance of $H_p^{(3)}$, along with all other $H_p^{(n)}$, on the left-hand side of Eq. (6.33), thereby coupling the equation for $H_p^{(2)}$ with those of the other $H_p^{(n)}$ terms. The second consequence comes from the second term on the right-hand side of Eq. (6.33). This term represents the contribution from the hydrodynamic stress generated by $\nabla\nabla\gamma$. However, one may perform several simplifications by noting first that all the terms in $H_p^{(n)}$ with $n \ge 4$ are of order $a^2/L^2$ or smaller. The same comment applies to the second term on the right-hand side of Eq. (6.33), which is of order $a^2/L^2$. In the current modeling hypothesis one may neglect all the $O(a^2/L^2)$ terms, hence Eq. (6.33) reduces to

$$H_p^{(2)} + \frac{30}{56} \frac{a}{\gamma} \nabla\gamma \cdot H_p^{(3)} = \frac{19\lambda + 16}{8(\lambda + 1)} \left( \frac{a\mu_f}{\gamma Ca} \right) E. \qquad (6.34)$$

The effect of the surface tension gradient is smaller than the effect of strain by a factor of $a/L$, yet it remains non-negligible within the current assumptions framework. To compute the components of $H^{(3)}$, we rely on the second moment of the momentum equation (see Eq. (F.17)), along with the closure relations for both the second moment of hydrodynamic forces and the first moment of internal shear rate. Given the complexity of the derivation, full details are provided in Appendix F. Since buoyancy effects are not central to our discussion, they are neglected in this calculation. The resulting expression is,

$$\begin{aligned} Ca \frac{8}{7} \gamma (H_p^{(3)})_{ijk} + Ca a(\nabla\gamma)_l \frac{16}{35} (H_p^{(2)}\delta)_{ijkl}^{p} + \frac{64}{105} H_{ijkl}^{(4)} = & \\ -\mu_f a^2 \frac{5\lambda + 2}{21(\lambda + 1)} (\delta\nabla^2\mathbf{u})_{ijk} + \mu_f a^2 \frac{2(\lambda - 1)}{21(\lambda + 1)} [(\delta\nabla^2\mathbf{u})_{ikj} + (\delta\nabla^2\mathbf{u})_{jki}] & \\ + \mu_f \frac{11\lambda + 10}{21(\lambda + 1)} [(\nabla\nabla\mathbf{u})_{ijk} + (\nabla\nabla\mathbf{u})_{ikj} + (\nabla\nabla\mathbf{u})_{kji}] & + \mu_f \frac{2(3\lambda + 2)}{5(\lambda + 1)} [(\delta\mathbf{u}_r)_{ikj} + (\delta\mathbf{u}_r)_{jki} - \frac{2}{3}(\delta\mathbf{u}_r)_{jki}] & \\ + a^3 \frac{4(16\lambda + 19)}{735(\lambda + 1)} (\nabla\nabla\nabla\gamma)_{ijk} & + a^3 \frac{2(66\lambda + 109)}{3675(\lambda + 1)} [(\delta\nabla^2\nabla\gamma)_{ikj} + (\delta\nabla^2\nabla\gamma)_{ikj}] & \\ - a^3 \frac{8(73\lambda + 102)}{11025(\lambda + 1)} (\delta\nabla^2\nabla\gamma)_{ijk} & + a \frac{4(2\lambda + 3)}{15(\lambda + 1)} [(\delta\nabla\gamma)_{ikj} + (\delta\nabla\gamma)_{jki} - \frac{2}{3}(\delta\nabla\gamma)_{ijk}]. & \end{aligned} \qquad (6.35)$$

At first glance, it may be seen that $H^{(3)}$ can be explicitly computed in terms of $\nabla\gamma$, $u_r$, $\nabla\nabla\gamma$ and $\nabla\nabla\mathbf{u}$. However, it is important to note that, in Stokes flow, the relative translation ($u_r$) does not contribute to the deformation of the droplet (Taylor and Acrivos, 1964; Nadim and Stone, 1991). To clarify this point, we consider the steady-state momentum balance on a force-free droplet (Eq. (6.1)), that is

$$(2 + 3\lambda)u_r = -\frac{\lambda}{2}a^2\nabla^2\mathbf{u} - \frac{2a}{3\mu_f}\nabla\gamma - \frac{2a^3}{30\mu_f}\nabla^2(\nabla\gamma). \qquad (6.36)$$

Inserting this expression in Eq. (6.35) yields,

$$\begin{aligned} (H_p^{(3)})_{ijk} + \frac{2a}{5\gamma}(\nabla\gamma)_l [(H_p^{(2)}\delta)_{ijkl}^{p} + \frac{4}{3}H_{ijkl}^{(4)}] + \dots = & \\ \frac{a^2\mu_f}{Ca\gamma} \frac{11\lambda + 10}{24(\lambda + 1)} [(\nabla\nabla\mathbf{u})_{ijk} + (\nabla\nabla\mathbf{u})_{ikj} + (\nabla\nabla\mathbf{u})_{kji} - \frac{1}{5}((\delta\nabla^2\mathbf{u})_{ijk} + (\delta\nabla^2\mathbf{u})_{ikj} + (\delta\nabla^2\mathbf{u})_{jki})] & \\ + \frac{a^3}{Ca\gamma} \frac{(16\lambda + 19)}{210(\lambda + 1)} (\nabla\nabla\nabla\gamma)_{ijk} & + \frac{a^3}{Ca\gamma} \frac{(11\lambda + 10)}{350(\lambda + 1)} [(\delta\nabla^2\nabla\gamma)_{ikj} + (\delta\nabla^2\nabla\gamma)_{ikj}] & \\ - \frac{a^3}{Ca\gamma} \frac{146\lambda + 155}{3150(\lambda + 1)} (\delta\nabla^2\nabla\gamma)_{ijk} & + \frac{a}{Ca\gamma} \frac{8}{15} ((\delta\nabla\gamma)_{kji} + (\delta\nabla\gamma)_{kij} - \frac{2}{3}(\delta\nabla\gamma)_{ijk}) & \end{aligned} \qquad (6.37)$$

Examining the first line, one can identify the results of Nadim and Stone (1991) for the deformation of a force-free droplet embedded in a quadratic flow.<fnref n="6" /> The second line represents the contributions arising from deformation due to the first and third derivatives of the surface tension gradient. All terms contributing to $H_p^{(3)}$ are of order $O(a/L)$, implying that in Eq. (6.34) the term involving $H_p^{(3)}$ scales as $O(a^2/L^2)$ and is therefore negligible. Including terms of order $O(a^2/L^2)$ significantly increases the complexity of the problem for two main reasons: (1) source terms involving gradients of $\gamma$ are present on the right-hand sides of Eqs. (6.33) and (6.35); and (2) the equations governing each $H_p^{(n)}$ become coupled.

Once the averaged deformation is determined via the moment of momentum equations, it becomes possible to express the closure terms of the hybrid model (i.e., from Eqs. (6.1) to (6.5)) in terms of the tensors $H_p^{(n)}$, as illustrated, for example, in Haber and Hetsroni (1971). This leads to a hybrid model that is accurate to $O(Ca)$ and will be the focus of future work.

## 7. Conclusion

In this work, we provided a comprehensive presentation of the averaged equations and the derivation of closure relations for complex dispersed phases. The general derivation, developed from Sections 2.2 to 4.1, offers a detailed explanation of the structure of the averaged system under the so-called hybrid formulation. To illustrate our methodology, Sections 5 and 6 examine the case of a dilute, viscous-dominated droplet suspension, incorporating effects due to non-uniform surface tension. We showed that surface tension gradients give rise to a momentum source, driving the Marangoni drift of the dispersed phase, and also introduce non-Newtonian behavior in the momentum equations. Additionally, we derived a method to compute the averaged droplet shape using the moment of momentum equations. A natural extension of this work involves introducing transport equations to describe the surfactant concentration or temperature field, thereby linking surface tension gradients to the corresponding scalar fields.

<footnote n="5">⁵ While the contribution from the mean cubic flow $\nabla^2 E$, could in principle, be accounted for using Faxen relations applied to the first hydrodynamic moment, this term is $a^2/L^2$ smaller than this other and is therefore negligible in the present context.</footnote>

<footnote n="6">⁶ Note that $-15H_p^{(3)}$ corresponds to Nadim and Stone (1991)'s $H^{(3)}$, thereby recovering the same result.</footnote>
There has been limited effort to extend existing results to finite Reynolds number regimes that fully address the closure problem, despite the fact that most flows of interest exhibit significant inertial effects. Stone et al. (2001) showed that in a dilute suspension of solid spherical particles immersed in a general linear flow, finite inertia gives rise to normal stresses. This prediction was extended to the case of drops by Raja et al. (2010). However, both studies focus exclusively on neutrally buoyant particles, thereby neglecting the influence of relative motion on the suspension rheology. To begin exploring how inertia may influence both suspension rheology and droplet deformation, it is necessary to determine at least the first-order correction in $O(Re)$ to the closure terms in Eqs. (5.26) and (5.32). Based on symmetry arguments — similar to those applied in Eq. (6.9) — one may state that the functional form of the Stresslet due to pure relative motion is given by,

$$ \left\langle \delta_p \oint_{\Gamma_a} \mathbf{r} \sigma_f^* \cdot \mathbf{n} d\Gamma \right\rangle - \left\langle \delta_p \int_{\Omega_a} 2\mu_f e_f^* d\Omega \right\rangle = \mathrm{Re}\phi [C_s^1 \mathbf{u}_r \mathbf{u}_r + C_s^2 (\mathbf{u}_r \cdot \mathbf{u}_r) \delta] + O(Re^2, \phi^2), \quad (7.1) $$

at first order in Reynolds number. Here, $C_s^1$ and $C_s^2$ are functions of $\lambda$ that must be determined. We conclude that at finite Reynolds values, relative motions between droplets induce additional stresses in the suspension, scaling as $\sim \mathbf{u}_r \mathbf{u}_r$, in addition to the Reynolds stress contribution given in (Eq. (6.9)). These relative motions also lead to droplet deformation, consistent with the predictions of (Eq. (5.32)) and in agreement with the results of Taylor and Acrivos (1964). In an ongoing work, we explore how droplet translation at finite Reynolds numbers affects the Stresslet. In particular, we determine the exact values of the constants $C_s^1$ and $C_s^2$ appearing in Eq. (7.1) (see also Fintzi, 2025).

### CRediT authorship contribution statement

**Nicolas Fintzi:** Methodology, Writing – review & editing, Conceptualization, Writing – original draft. **Jean-Lou Pierson:** Supervision, Conceptualization, Writing – review & editing, Writing – original draft, Methodology.

### Declaration of competing interest

The authors declare the following financial interests/personal relationships which may be considered as potential competing interests: Reports a relationship with that includes: Has patent pending to. If there are other authors, they declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### Acknowledgments

It is a great pleasure to dedicate this paper to Daniel Lhuillier whose work on the “hybrid formalism” has greatly influenced the vision of both co-authors and the present work. As a scientist, Daniel Lhuillier serves as a mentor for young researchers. He constantly challenges us to question the significance and novelty of our findings and encourages us to explore papers that the authors may have overlooked. His constant guidance have also been a source of motivation to complete this work. The authors are also grateful for many in-depth discussions with Professor Stéphane Popinet as well.

## Appendix A. Topological equations for $\delta_r$

In this appendix, we derive some useful relation related to $\delta_r$. In the sense of distribution (Appel, 2007),

$$ \langle \nabla \delta_r, \varphi \rangle = - \langle \delta_r, \nabla \varphi \rangle. \quad (A.1) $$

Hence,

$$ \langle \nabla \delta_r, \varphi \rangle = - \int_{\Gamma} \nabla \varphi d\Gamma. \quad (A.2) $$

Since $\mathbf{n}$ is a unit vector, we have $\nabla \mathbf{n} \cdot \mathbf{n} = 0$ and $\nabla \varphi = \nabla (\varphi \mathbf{n}) \cdot \mathbf{n}$, which yields

$$ \langle \nabla \delta_r, \varphi \rangle = - \int_{\Gamma} \nabla (\varphi \mathbf{n}) \cdot \mathbf{n} d\Gamma. \quad (A.3) $$

By applying the result of multiplying a distribution by a function we get (Appel, 2007)

$$ \langle \nabla \delta_r, \varphi \rangle = \langle \nabla (\delta_r \mathbf{n}) \cdot \mathbf{n}, \varphi \rangle. \quad (A.4) $$

With a slight abuse of notation the preceding relation may be rewritten as

$$ \nabla \delta_r = \nabla (\mathbf{n} \delta_r) \cdot \mathbf{n} \quad (A.5) $$

Another relation may be obtained for the surface gradient operator of $\delta_r$. By definition,

$$ \nabla_{\parallel} \delta_r = \nabla \delta_r - \mathbf{n} \cdot \nabla \delta_r \quad (A.6) $$

Since, $\mathbf{n} \cdot \nabla \delta_r = \nabla (\delta_r \mathbf{n}) \cdot \mathbf{n} - \nabla \mathbf{n} \cdot \mathbf{n} \delta_r$ and using Eq. (A.5) we obtain

$$ \nabla_{\parallel} \delta_r = (\nabla \mathbf{n} \cdot \mathbf{n}) \delta_r \quad (A.7) $$
## Appendix B. Generalized form of the conservation equation on the interface

In this appendix we derive Eq. (2.10). Multiplying Eq. (2.8) by $\delta_R$ yields

$$ \delta_R [\partial_t f_R^0 + f_R^0(\mathbf{u}_R^0 \cdot \mathbf{n})(\nabla \cdot \mathbf{n}) + \nabla_\parallel \cdot (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] = \delta_R s_I^0 - \delta_R \llbracket f_k^0(\mathbf{u}_I^0 - \mathbf{u}_k^0) + \mathbf{\Phi}_k^0 \rrbracket \cdot \mathbf{n}_k. \quad (B.1) $$

We first focus on the first two terms on the left-hand side of Eq. (B.1). Since, $\delta_R \partial_t f_R^0 = \partial_t(f_R^0 \delta_R) - f_R^0 \partial_t \delta_R$ and using Eq. (2.5) yields

$$ \delta_R [\partial_t f_R^0 + f_R^0(\mathbf{u}_R^0 \cdot \mathbf{n})(\nabla \cdot \mathbf{n})] = \partial_t(f_R^0 \delta_R) + f_R^0 \nabla \cdot [(\mathbf{u}_I^0 \cdot \mathbf{n})\mathbf{n}\delta_R]. \quad (B.2) $$

Since the gradient of $f_R^0$ lies in the plane parallel to the interface, we have $\nabla f_R^0 \cdot \mathbf{n} = 0$. Consequently, the expression $f_R^0 \nabla \cdot [(\mathbf{u}_I^0 \cdot \mathbf{n})\mathbf{n}\delta_R]$ simplifies to $\nabla \cdot [(\mathbf{u}_I^0 \cdot \mathbf{n})\mathbf{n}f_R^0 \delta_R]$ which yields

$$ \delta_R [\partial_t f_R^0 + f_R^0(\mathbf{u}_R^0 \cdot \mathbf{n})(\nabla \cdot \mathbf{n})] = \partial_t(f_R^0 \delta_R) + \nabla \cdot [(\mathbf{u}_I^0 \cdot \mathbf{n})\mathbf{n}f_R^0 \delta_R]. \quad (B.3) $$

Next, we turn our attention to the third term on the left-hand side of Eq. (B.1). By applying the chain rule

$$ \begin{aligned} \delta_R \nabla_\parallel \cdot (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) = & \nabla \cdot [\delta_R(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] \\ & - \mathbf{n}(\mathbf{n} \cdot \nabla) \cdot [\delta_R(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] - (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) \cdot \nabla_\parallel \delta_R. \end{aligned} \quad (B.4) $$

Since $(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) \cdot \nabla_\parallel \delta_R = (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) \cdot \nabla_\parallel \delta_R$ and using Eq. (A.7) we may rewrite the last term of Eq. (B.4) as $(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) \cdot (\nabla \mathbf{n} \cdot \mathbf{n}) \delta_R$. This yields,

$$ \begin{aligned} \delta_R \nabla_\parallel \cdot (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) = & \nabla \cdot [\delta_R(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] \\ & - \mathbf{n}(\mathbf{n} \cdot \nabla) \cdot [\delta_R(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] - (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) \cdot (\nabla \mathbf{n} \cdot \mathbf{n}) \delta_R. \end{aligned} \quad (B.5) $$

By noticing that $f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel = (\delta - \mathbf{n}\mathbf{n}) \cdot (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)$ we can expand the second term on the right hand side of Eq. (B.5) which yields

$$ -\mathbf{n}(\mathbf{n} \cdot \nabla) \cdot [\delta_R(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] = (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) \cdot (\nabla \mathbf{n} \cdot \mathbf{n}) \delta_R \quad (B.6) $$

To derive the last equation, we used the following properties in indicial notation: $n_i n_j (\delta_{ik} - n_i n_k) = 0$ and $n_i n_j \delta_j (n_i n_k) = \partial_j (n_k) n_j$. Inserting Eq. (B.6) in Eq. (B.5) yields

$$ \delta_R \nabla_\parallel \cdot (f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel) = \nabla \cdot [\delta_R(f_R^0 \mathbf{u}_R^0 - \mathbf{\Phi}_R^0_\parallel)] \quad (B.7) $$

Inserting Eqs. (B.3) and (B.7) in Eq. (B.1) gives

$$ \partial_t(\delta_R f_R^0) + \nabla \cdot (\delta_R f_R^0 \mathbf{u}_I^0 - \delta_R \mathbf{\Phi}_R^0_\parallel) = \delta_R s_I^0 - \delta_R \llbracket f_k^0(\mathbf{u}_I^0 - \mathbf{u}_k^0) + \mathbf{\Phi}_k^0 \rrbracket \cdot \mathbf{n}_k, \quad (B.8) $$

which is the final form of the conservation equation on the interface.

## Appendix C. Arbitrary order moments equations

In this appendix, we extend the Lagrangian conservation laws to an arbitrary order moment equation. Let us first define the arbitrary moment of the Eulerian field $f_d^0$, in index notation it reads,

$$ (\mathbf{q}_a^{(n)})_{i_1 \dots i_n} = \int_{\Omega_a} \prod_{m=1}^n r_{im} f_d^0 d\Omega, $$

where we recall that $r_{im} = x_{im} - (x_a)_{im}$. Applying the Reynolds transport theorem (Eq. (3.10)) we show that:

$$ \frac{d}{dt} (\mathbf{q}_a^{(n)})_{i_1 \dots i_n} = \int_{\Omega_a} \left[ \partial_t \left( f_d^0 \prod_{m=1}^n r_{im} \right) + \nabla \cdot \left( \mathbf{u}_d^0 f_d^0 \prod_{m=1}^n r_{im} \right) \right] d\Omega + \oint_{\Gamma_a} \prod_{m=1}^n r_{im} f_d^0 (\mathbf{u}_R^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma. \quad (C.1) $$

Applying the product rule to the terms inside the first integral results in the following expression:

$$ \frac{d}{dt} (\mathbf{q}_a^{(n)})_{i_1 \dots i_n} = \int_{\Omega_a} f_d^0 \left[ \partial_t \left( \prod_{m=1}^n r_{im} \right) + (\mathbf{u}_d^0 \cdot \nabla) \left( \prod_{m=1}^n r_{im} \right) \right] d\Omega + \int_{\Omega_a} \prod_{m=1}^n r_{im} \left[ \partial_t f_d^0 + \nabla \cdot (\mathbf{u}_d^0 f_d^0) \right] d\Omega + \oint_{\Gamma_a} \prod_{m=1}^n r_{im} f_d^0 (\mathbf{u}_R^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma. \quad (C.2) $$

By applying the product rule, once again, to the terms within the first integral on, noting that $\partial_t r_{i_e} = (\mathbf{w}_d^0)_{i_e}$, and utilizing the conservation Eq. (2.7) on the second integral, leads us to the relation:

$$ \frac{d}{dt} (\mathbf{q}_a^{(n)})_{i_1 \dots i_n} = \sum_{e=1}^{n} \int_{\Omega_a} f_d^0 \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{im} (\mathbf{w}_d^0)_{i_e} d\Omega + \int_{\Omega_a} \prod_{m=1}^{n} r_{im} (\nabla \cdot \mathbf{\Phi}_d^0) d\Omega + \int_{\Gamma_a} \prod_{m=1}^{n} r_{im} f_d^0 (\mathbf{u}_R^0 - \mathbf{u}_d^0) \cdot \mathbf{n}_d d\Gamma. \quad (C.3) $$

The second term on the right-hand side of this equation can be reformulated as,

$$ \begin{aligned} \int_{\Omega_a} \prod_{m=1}^{n} r_{im} (\nabla \cdot \mathbf{\Phi}_d^0) d\Omega &= \int_{\Omega_a} \nabla \cdot \left( \prod_{m=1}^{n} r_{im} \mathbf{\Phi}_d^0 \right) d\Omega - \int_{\Omega_a} \mathbf{\Phi}_d^0 \cdot \nabla \left( \prod_{m=1}^{n} r_{im} \right) d\Omega \\ &= \oint_{\Gamma_a} \prod_{m=1}^{n} r_{im} (\mathbf{\Phi}_d^0 \cdot \mathbf{n}_d) d\Gamma - \sum_{e=1}^{n} \int_{\Omega_a} (\mathbf{\Phi}_d^0)_{i_e} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{im} d\Omega \end{aligned} $$
Including this relation into the former equation yields,

$$ \frac{d}{dt}(\mathbf{q}_{\alpha}^{(n)})_{i_1 \dots i_n} = \sum_{e=1}^{n} \oint_{\Omega_{\alpha}} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{i_m} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_e} d\Omega + \oint_{\Omega_{\alpha}} \prod_{m=1}^{n} r_{i_m} (f_d^0 \mathbf{u}_d^0 + f_d^0 (\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0)) \cdot \mathbf{n}_d d\Gamma. \quad (C.4) $$

This represents the final form of the Lagrangian conservation equation for the nth order moment of the quantity $f_d^0$ within the particle. When $f_d^0$ is a scalar quantity, this expression indicates that $\frac{d}{dt}\mathbf{q}_\alpha^{(n)}$ is entirely symmetric since it involves only sums of products of the $r_{i_n}$ times a scalar quantity. Regarding the seemingly non-symmetric terms involving $\mathbf{w}_d^0$ and $\mathbf{\Phi}_d^0$, it is important to note that, due to the summation operator's presence, these terms also become symmetric.

Regarding the surface property conservation equations, the derivation is similar and will not be displayed here. The result yields,

$$ \frac{d}{dt}(\mathbf{q}_{\alpha\Gamma}^{(n)})_{i_1 \dots i_n} = \sum_{e=1}^{n} \oint_{\Gamma_{\alpha}} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{i_m} (f_{\Gamma}^0 \mathbf{w}_{\Gamma}^0 - \mathbf{\Phi}_{\parallel\Gamma}^0)_{i_e} d\Gamma + \oint_{\Gamma_{\alpha}} \prod_{m=1}^{n} r_{i_m} (\mathbf{s}_{\Gamma}^0)_{k} d\Gamma + \oint_{\Gamma_{\alpha}} \prod_{m=1}^{n} r_{i_m} \|\mathbf{\Phi}_{k}^0 + f_{k}^0 (\mathbf{u}_{\Gamma}^0 - \mathbf{u}_{k}^0) \cdot \mathbf{n}_{d}\| \cdot \mathbf{n}_{k} d\Gamma. \quad (C.5) $$

where we have defined,

$$ (\mathbf{q}_{\alpha\Gamma}^{(n)})_{i_1 \dots i_n} = \oint_{\Gamma_{\alpha}} \prod_{m=1}^{n} r_{i_m} f_{\Gamma}^0 d\Gamma. $$

Summing Eqs. (C.4) and (C.5) one obtains the equation for the total nth order moment, $\mathbf{Q}_\alpha^{(n)} = \mathbf{q}_\alpha^{(n)} + \mathbf{q}_d^{(n)}$, namely,

$$ \begin{aligned} \frac{d}{dt}(\mathbf{Q}_{\alpha}^{(n)})_{i_1 \dots i_n} = & \sum_{e=1}^{n} \oint_{\Omega_{\alpha}} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{i_m} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_e} d\Omega + \oint_{\Omega_{\alpha}} \prod_{m=1}^{n} r_{i_m} (\mathbf{s}_d^0)_{k} d\Omega \\ & + \sum_{e=1}^{n} \oint_{\Gamma_{\alpha}} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{i_m} (f_{\Gamma}^0 \mathbf{w}_{\Gamma}^0 - \mathbf{\Phi}_{\parallel\Gamma}^0)_{i_e} d\Gamma + \oint_{\Gamma_{\alpha}} \prod_{m=1}^{n} r_{i_m} (\mathbf{s}_{\Gamma}^0)_{k} d\Gamma \\ & + \oint_{\Gamma_{\alpha}} \prod_{m=1}^{n} r_{i_m} ([\mathbf{\Phi}_f^0 + \mathbf{f}_f^0 (\mathbf{u}_f^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d)_k d\Gamma. \end{aligned} \quad (C.6) $$

By averaging Eq. (C.6), we obtain the particle averaged equation for $\langle \delta_p \mathbf{Q}_\alpha^{(n)} \rangle = n_p \mathbf{Q}_p^{(n)}$, it yields,

$$ \begin{aligned} & \partial_t \langle \delta_p (\mathbf{Q}_\alpha^{(n)})_{i_1 \dots i_n}^\alpha \rangle + \nabla \cdot \langle \delta_p \mathbf{u}_\alpha (\mathbf{Q}_\alpha^{(n)})_{i_1 \dots i_n}^\alpha \rangle = \sum_{e=1}^{n} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{i_m} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_e} d\Omega \right\rangle \\ & + \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} (\mathbf{s}_d^0)_k d\Omega \right\rangle + \sum_{e=1}^{n} \left\langle \delta_p \oint_{\Gamma_\alpha} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{i_m} (f_\Gamma^0 \mathbf{w}_\Gamma^0 - \mathbf{\Phi}_{\parallel\Gamma}^0)_{i_e} d\Gamma \right\rangle \\ & + \left\langle \delta_p \oint_{\Gamma_\alpha} \prod_{m=1}^{n} r_{i_m} (\mathbf{s}_\Gamma^0)_k d\Gamma \right\rangle + \left\langle \delta_p \oint_{\Gamma_\alpha} \prod_{m=1}^{n} r_{i_m} ([\mathbf{\Phi}_f^0 + \mathbf{f}_f^0 (\mathbf{u}_f^0 - \mathbf{u}_f^0)] \cdot \mathbf{n}_d)_k d\Gamma \right\rangle. \end{aligned} \quad (C.7) $$

## Appendix D. Expansion of phase and surface quantities as series of moments

We first provide a concise proof of the following relation

$$ \delta(\mathbf{x} - \mathbf{x}_{\alpha} - \mathbf{r}) = \delta(\mathbf{x} - \mathbf{x}_{\alpha}) - \mathbf{r} \cdot \nabla \delta(\mathbf{x} - \mathbf{x}_{\alpha}) + \frac{1}{2} \mathbf{r}\mathbf{r} : \nabla\nabla\delta(\mathbf{x} - \mathbf{x}_{\alpha}) - \dots \quad (D.1) $$

which is widely used in the literature (see for example Zhang et al. (2023)). We begin by applying the Dirac delta function $\delta(\mathbf{x} - \mathbf{x}_{\alpha} - \mathbf{r})$, to a test function $\varphi(\mathbf{x})$. This results in the following expression,

$$ \langle \delta(\mathbf{x} - \mathbf{x}_{\alpha} - \mathbf{r}), \varphi(\mathbf{x}) \rangle = \varphi(\mathbf{x}_{\alpha} + \mathbf{r}). $$

By applying the Taylor expansion to the test function $\varphi(\mathbf{x}_\alpha + \mathbf{r})$ around the point $\mathbf{r} = 0$, we obtain the following expression

$$ \langle \delta(\mathbf{x} - \mathbf{x}_{\alpha} - \mathbf{r}), \varphi(\mathbf{x}) \rangle = \varphi(\mathbf{x}_{\alpha}) + \mathbf{r} \cdot \nabla \varphi|_{\mathbf{x}_{\alpha}} + \frac{1}{2} \mathbf{r}\mathbf{r} : \nabla\nabla \varphi|_{\mathbf{x}_{\alpha}} + \dots \quad (D.2) $$

By the definition of the Dirac delta function and its derivative (Appel, 2007), we have $(\delta_\alpha, \varphi) = \varphi(\mathbf{x}_\alpha)$, $(\nabla\delta_\alpha, \varphi) = -\nabla\varphi|_{\mathbf{x}_\alpha}$, ... Substituting these definitions into Eq. (D.2), we obtain in the distributional sense Eq. (D.1).

Let us now turn our attention to the demonstration of Eq. (4.1). We can show that,

$$ \langle (f_d^0 \chi_\alpha)(\mathbf{x}), \varphi(\mathbf{x}) \rangle = \langle (f_d^0 \chi_\alpha)(\mathbf{x}_\alpha + \mathbf{r}), \varphi(\mathbf{x}_\alpha + \mathbf{r}) \rangle \quad (D.3) $$

$$ = \langle (f_d^0 \chi_\alpha)(\mathbf{x}_\alpha + \mathbf{r}), \langle \delta(\mathbf{x} - \mathbf{x}_\alpha - \mathbf{r}), \varphi(\mathbf{x}) \rangle \rangle. \quad (D.4) $$

The first equality follows directly from the change of variables $\mathbf{x} = \mathbf{x}_\alpha + \mathbf{r}$. The second equality holds by the definition of the Dirac delta functions and of the convolution product in the distributional sense (Appel, 2007). Substituting Eq. (D.1) into the previous equation gives,

$$ f_d^0 \chi_\alpha = \delta_\alpha \int_{\Omega_\alpha} f_d^0 d\Omega - \nabla \cdot \left( \delta_\alpha \int_{\Omega_\alpha} \mathbf{r} f_d^0 d\Omega \right) + \frac{1}{2} \nabla \cdot \nabla \cdot \left( \delta_\alpha \int_{\Omega_\alpha} \mathbf{r} \mathbf{r} f_d^0 d\Omega \right) - \dots \quad (D.5) $$

The previous proof can be easily extended to surface quantities by following the same steps as in Eqs. (D.3) and (D.4), leading to

$$ \langle f_{\Gamma}^{0} \delta_{f_{\alpha}}, \varphi(\mathbf{x}) \rangle = \langle (f_{\Gamma}^{0} \delta_{f_{\alpha}})(\mathbf{x}_{\alpha} + \mathbf{r}), \langle \delta(\mathbf{x} - \mathbf{x}_{\alpha} - \mathbf{r}), \varphi(\mathbf{x}) \rangle \rangle, \quad (D.6) $$
from which we obtain

$$f_{\Gamma}^{0} \delta_{\Gamma \alpha} = \delta_{\alpha} \oint_{\Gamma_{\alpha}} f_{\Gamma}^{0} d\Gamma - \nabla \cdot \left( \delta_{\alpha} \oint_{\Gamma_{\alpha}} \mathbf{r} f_{\Gamma}^{0} d\Gamma \right) + \frac{1}{2} \nabla \nabla : \left( \delta_{\alpha} \oint_{\Gamma_{\alpha}} \mathbf{r} \mathbf{r} f_{\Gamma}^{0} d\Gamma \right) - \dots \quad (D.7)$$

## Appendix E. Arbitrary order equivalence

In this appendix, we provide a general proof of Eq. (4.9). Let us begin by re-writing the phase averaged equation:

$$\partial_t \langle \chi_d f_d^0 \rangle = \nabla \cdot \left\langle \chi_d \Phi_d^0 - \chi_d f_d^0 \mathbf{u}_d^0 \right\rangle + \langle \chi_d s_d^0 \rangle + \langle \delta_{\Gamma} [\Phi_d^0 + f_d^0 (\mathbf{u}_{\Gamma}^0 - \mathbf{u}_d^0)] \cdot \mathbf{n}_d \rangle. \quad (E.1)$$

Our objective here is to demonstrate how Eq. (E.1) is related to the Lagrangian moments' equation, given by Eq. (C.4). The first step is to expand each term of Eq. (E.1) using the relation Eq. (4.6) which directly gives,

$$\begin{aligned} 0 = & -\partial_t \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} f_d^0 d\Omega \right\rangle \\ & + \nabla \cdot \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} (\Phi_d^0 - f_d^0 \mathbf{u}_d^0) d\Omega \right\rangle \\ & + \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} s_d^0 d\Omega \right\rangle \\ & + \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \oint_{\Gamma_\alpha} \prod_{m=1}^{n} r_{i_m} [\Phi_d^0 + f_d^0 (\mathbf{u}_\Gamma^0 - \mathbf{u}_d^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle. \end{aligned}$$

The third term can be reformulated using the decomposition: $\mathbf{u}_d^0 = \mathbf{u}_\alpha + \mathbf{w}_d^0$, which gives,

$$\begin{aligned} & \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} f_d^0 \mathbf{u}_d^0 d\Omega \right\rangle \\ & = \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \mathbf{u}_\alpha \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} f_d^0 d\Omega \right\rangle \\ & + \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} f_d^0 \mathbf{w}_d^0 d\Omega \right\rangle. \end{aligned} \quad (E.2)$$

Injecting this formulation in the former equation yields,

$$\partial_t \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} f_d^0 d\Omega \right\rangle \quad (E.3)$$

$$+ \nabla \cdot \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \mathbf{u}_\alpha \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} f_d^0 d\Omega \right\rangle \quad (E.4)$$

$$= \nabla \cdot \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} (\mathbf{\Phi}_d^0 - f_d^0 \mathbf{w}_d^0) d\Omega \right\rangle \quad (E.5)$$

$$+ \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} s_d^0 d\Omega \right\rangle \quad (E.6)$$

$$+ \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left\langle \delta_p \oint_{\Gamma_\alpha} \prod_{m=1}^{n} r_{i_m} [\mathbf{\Phi}_d^0 + f_d^0 (\mathbf{u}_I^0 - \mathbf{u}_d^0)] \cdot \mathbf{n}_d d\Gamma \right\rangle. \quad (E.7)$$

Then, note that the first term on the right-hand side can be re-written when evaluated at the order $n-1$, as follows,

$$\begin{aligned} \nabla \cdot \frac{(-1)^{(n-1)}}{(n-1)!} \prod_{m=1}^{(n-1)} \partial_{i_m} \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{(n-1)} r_{i_m} (\mathbf{\Phi}_d^0 - f_d^0 \mathbf{w}_d^0) d\Omega \right\rangle = \\ \frac{(-1)^n}{(n)!} \prod_{m=1}^{n} \partial_{i_m} n \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n-1} r_{i_m} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_m} d\Omega \right\rangle \end{aligned}$$

Upon using that expression in Eq. (E.7) we can factor out the gradient operators, i.e. the $\frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m}$, which gives,

$$\begin{aligned} 0 = & \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left[ -\partial_t \langle \delta_p (\mathbf{q}_\alpha^{(n)})_{i_1 \dots i_n} \rangle - \nabla \cdot \left\langle \delta_p \mathbf{u}_\alpha (\mathbf{q}_\alpha^{(n)})_{i_1 \dots i_n} \right\rangle \right. \\ & \left. + n \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n-1} r_{i_m} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0) d\Omega \right\rangle + \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} s_d^0 d\Omega \right\rangle \right. \\ & \left. + \left\langle \delta_p \int_{\Omega_\alpha} \prod_{m=1}^{n} r_{i_m} [\mathbf{\Phi}_d^0 + f_d^0 (\mathbf{u}_I^0 - \mathbf{u}_d^0)] \cdot \mathbf{n}_d d\Omega \right\rangle \right]. \end{aligned} \quad (E.8)$$
At this stage, one might immediately recognize Eq. (C.4) in the square bracket. However, note that the third term of Eq. (E.8) differs from the second term of Eq. (C.4). Indeed,

$$n \left\langle \delta_p \int_{\Omega_a} \prod_{m=1}^{n-1} r_{im} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_n} d\Omega \right\rangle \neq \sum_{e=1}^{n} \left\langle \int_{\Omega_a} \prod_{m \neq e}^{n} r_{im} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_e} d\Omega \right\rangle. \quad (E.9)$$

Even if the above inequality holds as it is, it must be understood that what matter in Eq. (E.8) is the gradient of that term. Additionally, the contraction of one of these terms with the gradient operator $\prod_{m=1}^{n} \partial_{i_m}$ makes the skew-symmetric part of these tensors (in the indices $i_1...i_n$) having a vanishing contribution to the expression. Thus, one might apply the operator $\prod_{m=1}^{n} \partial_{i_m}$ on each side of Eq. (E.9) and by noticing that the skew-symmetric part of the left-hand side of Eq. (E.9) vanish, this leads to,

$$\prod_{m=1}^{n} \partial_{i_m} \left[ n \left\langle \delta_p \int_{\Omega_a} \prod_{m=1}^{n-1} r_{im} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_n} d\Omega \right\rangle \right] = \prod_{m=1}^{n} \partial_{i_m} \left[ \sum_{e=1}^{n} \left\langle \int_{\Omega_a} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{im} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_e} d\Omega \right\rangle \right]. \quad (E.10)$$

Injecting this last equality in Eq. (E.8) and using Eq. (C.5) to reformulate the interfacial term directly gives,

$$0 = \frac{(-1)^n}{n!} \prod_{m=1}^{n} \partial_{i_m} \left[ -\partial_t \langle \delta_p (\mathbf{Q}_a^{(n)})_{i_1 \dots i_n}^{\alpha} \rangle - \nabla \cdot \langle \delta_p \mathbf{u}_a (\mathbf{Q}_a^{(n)})_{i_1 \dots i_n}^{\alpha} \rangle \right] \\ + \sum_{e=1}^{n} \left\langle \delta_p \int_{\Omega_a} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{im} (f_d^0 \mathbf{w}_d^0 - \mathbf{\Phi}_d^0)_{i_e} d\Omega \right\rangle + \left\langle \delta_p \int_{\Omega_a} \prod_{m=1}^{n} r_{im} (\mathbf{s}_d^0)_k d\Omega \right\rangle \\ + \sum_{e=1}^{n} \left\langle \delta_p \oint_{\Gamma_a} \prod_{\substack{m=1 \\ m \neq e}}^{n} r_{im} (f_\Gamma^0 \mathbf{w}_\Gamma^0 - \mathbf{\Phi}_{||\Gamma}^0)_{i_e} d\Gamma \right\rangle + \left\langle \delta_p \oint_{\Gamma_a} \prod_{m=1}^{n} r_{im} (\mathbf{s}_\Gamma^0)_k d\Gamma \right\rangle \\ + \left\langle \delta_p \oint_{\Gamma_a} \prod_{m=1}^{n} r_{im} (|\mathbf{\Phi}_f^0 + \mathbf{r}_f^0 (\mathbf{u}_\Gamma^0 - \mathbf{u}_f^0)|) \cdot \mathbf{n}_d)_k d\Gamma \right\rangle, \quad (E.11)$$

which proves Eq. (4.9).

## Appendix F. Singularity solution for droplets with surface tension gradient

At the leading order in droplet volume fraction, the closure problem is equivalent to that of an isolated droplet in an infinite medium (Hinch, 1977). Hence, we consider the problem of an isolated droplet, immersed in an arbitrary quadratic flow without the presence of fluid inertia. The disturbances pressure and velocity field are noted $\mathbf{u}_{out}$, $\mathbf{u}_{in}$, $p_{out}$ and $p_{in}$, for the velocity outside the droplet, the velocity inside the droplet, the pressure outside the droplet and the pressure inside the droplet, respectively. The corresponding “undisturbed field” are noted $\mathbf{u}$ and $p_f$, corresponding to the ensemble averaged mixture velocity and averaged pressure of the continuous phase in the main text.

Under these hypotheses, in dimensionless form, $\mathbf{u}_{out}$, $\mathbf{u}_{in}$, $p_f$ and $p_o$ are governed by the Stokes equations, namely,

$$\nabla \cdot \mathbf{u}_{in} = 0 \qquad \nabla \cdot \mathbf{u}_{out} = 0 \qquad (F.1)$$

$$\nabla^2 \mathbf{u}_{in} = \nabla p_{in} \qquad \nabla^2 \mathbf{u}_{out} = \nabla p_{out} \qquad (F.2)$$

At the surface of the droplet ($r=1$) the condition of continuity of velocity and continuity of the normal stresses read,

$$\mathbf{u}_{out} - \mathbf{u}_{in} = 0 \qquad (F.3)$$

$$(\mathbf{u}_{in} + \mathbf{u}_r) \cdot \mathbf{n} = 0 \qquad (F.4)$$

$$\mathbf{n} \cdot 2(\mathbf{e}_{out} - \lambda \mathbf{e}_{in} + \mathbf{E} - \lambda \mathbf{E}) \cdot (\mathbf{\delta} - \mathbf{n}\mathbf{n}) = \mathbf{b} \cdot (\mathbf{\delta} - \mathbf{n}\mathbf{n}) \qquad (F.5)$$

where $\mathbf{u}_r = \mathbf{u} - \mathbf{u}_p$ corresponds to the velocity of the mixture with respect to the droplet mean velocity. We recall that $\mathbf{e}_{in/out} = (\nabla \mathbf{u}_{in/out} + \nabla^\dagger \mathbf{u}_{in/out})/2$ and $\mathbf{E} = (\nabla \mathbf{u} + \nabla^\dagger \mathbf{u})/2$, which correspond to the rate of strain inside or around or droplet, and to the ensemble averaged rate of strain, respectively. In Eq. (F.5) the vector $\mathbf{b}$ correspond to the dimensionless tangential stress jump across the interface and reads,

$$\mathbf{b} = \frac{a \nabla \gamma}{\mu_f U} \qquad (F.6)$$

Far from the droplet (centered at $\mathbf{r}=0$) the velocity and pressure fields satisfy,

$$\lim_{r \to \infty} (\mathbf{u}_{out} \cdot \mathbf{p}_{out}) = 0. \qquad (F.7)$$

In general $\mathbf{b}$ as well as $\mathbf{u}$ are not constant vectors and vary across space as a function of $\mathbf{r}$. However, it can be assumed that these quantities are slowly varying functions of space at the scale of a droplet radius, hence we may introduce the simplifications

$$\mathbf{u}(\mathbf{n}) = \mathbf{u}_r|_{\mathbf{r}=0} + \mathbf{r} \cdot \nabla \mathbf{u}|_{\mathbf{r}=0} + \frac{1}{2} \mathbf{r} \mathbf{r} : \nabla \nabla \mathbf{u}|_{\mathbf{r}=0} + \dots$$

$$\mathbf{E}(\mathbf{n}) = \mathbf{E}|_{\mathbf{r}=0} + \mathbf{r} \cdot \nabla \mathbf{E}|_{\mathbf{r}=0} + \frac{1}{2} \mathbf{r} \mathbf{r} : \nabla \nabla \mathbf{E}|_{\mathbf{r}=0} + \dots$$

$$\mathbf{b}(\mathbf{n}) = \mathbf{b}|_{\mathbf{r}=0} + \mathbf{r} \cdot \nabla \mathbf{b}|_{\mathbf{r}=0} + \frac{1}{2} \mathbf{r} \mathbf{r} : \nabla \nabla \mathbf{b}|_{\mathbf{r}=0} + \dots$$
Note that because **b** ~ $\nabla\gamma$ the tensor: $\nabla\mathbf{b}$ and $\nabla\nabla\mathbf{b}$ are symmetric tensors. Injecting, these expressions into Eqs. (F.3) to (F.5) one deduce that the velocity and pressure fields are linearly related to $\mathbf{u}|_{\mathbf{r}=0}$, $\mathbf{b}|_{\mathbf{r}=0}$, and to the gradients of these vectors.

According to the linearity of the Stokes equations we deduce that $\mathbf{u}_{l/o}$ and $p_{l/o}$ must be linear combination of spherical harmonics proportional to $\mathbf{u}|_{\mathbf{x}=0}$, $\nabla\gamma|_{\mathbf{x}=0}$, and their derivatives. Hence, we can write that the disturbance velocity and pressure fields are given by,

$$ \begin{pmatrix} \mathbf{u}_o \\ p_o \\ \mathbf{u}_i \\ p_i \end{pmatrix} = \begin{pmatrix} \mathbf{U}_o^{(1)} + \mathbf{U}_o^{(2)} \cdot \nabla + \mathbf{U}_o^{(3)} : \nabla\nabla & \mathbf{U}_o^{(b-1)} + \mathbf{U}_o^{(b-2)} \cdot \nabla + \mathbf{U}_o^{(b-3)} : \nabla\nabla \\ \mathbf{P}_o^{(1)} + \mathbf{P}_o^{(2)} \cdot \nabla + \mathbf{P}_o^{(3)} : \nabla\nabla & \mathbf{P}_o^{(b-1)} + \mathbf{P}_o^{(b-2)} \cdot \nabla + \mathbf{P}_o^{(b-3)} : \nabla\nabla \\ \mathbf{U}_i^{(1)} + \mathbf{U}_i^{(2)} \cdot \nabla + \mathbf{U}_i^{(3)} : \nabla\nabla & \mathbf{U}_i^{(b-1)} + \mathbf{U}_i^{(b-2)} \cdot \nabla + \mathbf{U}_i^{(b-3)} : \nabla\nabla \\ \mathbf{P}_i^{(1)} + \mathbf{P}_i^{(2)} \cdot \nabla + \mathbf{P}_i^{(3)} : \nabla\nabla & \mathbf{P}_i^{(b-1)} + \mathbf{P}_i^{(b-2)} \cdot \nabla + \mathbf{P}_i^{(b-3)} : \nabla\nabla \end{pmatrix} \cdot \begin{pmatrix} \mathbf{u}_r \\ \mathbf{b} \end{pmatrix} $$

The tensors $\mathbf{U}_{o/i}^{(n)}$ (resp. $\mathbf{P}_{o/i}^{(n)}$) are $(n+1)$th (resp. $n$th) order tensors, solely function of $\mathbf{r}$ and the viscosity ratio $\lambda$. The tensors $\mathbf{U}_{o/i}^{(1)}$ and $\mathbf{P}_{o/i}^{(1)}$ represent the functional form of the famous Hadamard-Ribczynski solution (Pozrikidis et al., 1992; Kim and Karrila, 2013). The dependency of $\mathbf{u}_{l/o}$ and $p_{l/o}$ with the mean gradient velocity $\nabla\mathbf{u}$ is given by the tensor $\mathbf{U}_{o/i}^{(2)}$ and $\mathbf{P}_{o/i}^{(2)}$. The disturbance velocity and pressure fields of a drop immersed in pure linear flow can be found in Rallison (1978), Leal (2007) and Raja et al. (2010). Finally, the disturbance field of a droplet immersed in a pure quadratic flow is defined through the tensor $\mathbf{U}_{o/i}^{(3)}$ and $\mathbf{P}_{o/i}^{(3)}$. The definition of the former tensors may be found in Nadim and Stone (1991).

The tensors $\mathbf{U}_{o/i}^{(b-n)}$ (resp. $\mathbf{P}_{o/i}^{(b-n)}$) are $(n+1)$th (resp. $n$th) order tensors, solely function of $\mathbf{r}$ and the viscosity ratio $\lambda$ as well. The disturbance fields induced due to a constant, linear, or quadratic stress jump at the interface of the droplet is less common in the literature. For instance one might refer to Subramanian (1985) and Leal (2007) to compute the disturbance fields due to a constant surface tension gradient, hence defining the expression of $\mathbf{U}_{o/i}^{(b-1)}$ and $\mathbf{P}_{o/i}^{(b-1)}$. Then to express the disturbance fields in terms of $\nabla\mathbf{b}$ or higher order derivatives one might refer to the solution given in Raja et al. (2010, Appendix C). Without loss of generality we consider only the symmetric traceless part of $(\nabla\mathbf{b})_{jk}$ in the following. Based on the classic method of spherical harmonic (Leal, 2007, chapter 8) we could derive the expression for the disturbance fields due to an imposed $\nabla\nabla\mathbf{b}$, it reads,

$$ \begin{aligned} (\mathbf{U}_o^{(b-3)})_{ijk} = & -\frac{7r^4 + 10r^2 - 9}{1260(\lambda + 1)r^5} (\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj}) - \frac{7r^4 - 60r^2 + 45}{1260(\lambda + 1)r^7} (x_i x_j \delta_{kl} + x_i x_k \delta_{jl} + x_i x_l \delta_{kj}) \\ & + \frac{r^2 - 3}{84(\lambda + 1)r^7} (\delta_{ij}x_k x_l + \delta_{ik}x_j x_l + \delta_{il}x_k x_l) - \frac{5r^2 - 7}{28(\lambda + 1)r^9} x_i x_j x_k x_l \end{aligned} $$

$$ (\mathbf{P}_o^{(b-3)})_{jkl} = -\frac{7r^2 - 45}{630(\lambda + 1)r^5} (\delta_{jk}x_l + \delta_{jl}x_k + \delta_{lk}x_j) - \frac{5}{14(\lambda + 1)r^7} x_j x_k x_l $$

$$ \begin{aligned} (\mathbf{U}_i^{(b-3)})_{ijk} = & \frac{9r^4 - 20r^2 + 7}{630(\lambda + 1)} (\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj}) + \frac{9r^2 - 5}{630(\lambda + 1)} (x_i x_j \delta_{kl} + x_i x_k \delta_{jl} + x_i x_l \delta_{kj}) \\ & - \frac{3r^2 - 2}{42(\lambda + 1)} (\delta_{ij}x_k x_l + \delta_{ik}x_j x_l + \delta_{il}x_k x_l) + \frac{1}{14(\lambda + 1)} x_i x_j x_k x_l \end{aligned} $$

$$ (\mathbf{P}_i^{(b-3)})_{jkl} = \frac{54r^2 - 35}{315(\lambda + 1)r} (\delta_{jk}x_l + \delta_{jl}x_k + \delta_{lk}x_j) - \frac{6}{7(\lambda + 1)} x_j x_k x_l $$

It is implied that each of these tensors must be contracted with the full symmetric tensor $(\nabla\nabla\mathbf{b})_{ijk} \sim (\nabla\nabla\nabla\gamma)_{ijk}$ to obtain the velocity or pressure fields.

Hence, based on these expressions and the solutions provided in the literature, one can compute the force traction at the surface of the droplet $(\sigma_{out} \cdot \mathbf{n})$, as well as the internal shear rate $(\mathbf{e}_{in})$, and obtains the following expressions,

$$ \int_{\Gamma_a} \sigma_{out} \cdot \mathbf{n} d\Gamma = 2\pi \frac{2+3\lambda}{1+\lambda} \mathbf{u}_r + \pi \frac{\lambda}{\lambda+1} \nabla^2 \mathbf{u} + \frac{4\pi}{3} \frac{1}{\lambda+1} \mathbf{b} + \frac{2\pi}{15(\lambda+1)} \nabla^2 \mathbf{b} \quad (F.8) $$

$$ \int_{\Gamma_a} \mathbf{x} \sigma_{out} \cdot \mathbf{n} d\Gamma = \frac{2\pi(5\lambda + 2)}{5(\lambda + 1)} [\nabla \mathbf{u} + (\nabla \mathbf{u})^{\dagger}] + \frac{12\pi}{25(\lambda + 1)} (\nabla \mathbf{b} - \frac{1}{3} \delta \nabla \cdot \mathbf{b}) \quad (F.9) $$

$$ \int_{\Gamma_a} \mathbf{x}\mathbf{x}\sigma_{out} \cdot \mathbf{n} d\Gamma = \frac{4\pi}{5(\lambda+1)} (\mathbf{u}_r \delta + \delta^{\dagger} \mathbf{u}_r \delta) + \frac{2\pi(5\lambda + 2)}{5(\lambda + 1)} \delta\mathbf{u}_r - \frac{8\pi}{15(\lambda + 1)} (\mathbf{b}\delta + \delta^{\dagger}\mathbf{b}\delta) + \frac{4\pi}{5(\lambda + 1)} \delta\mathbf{b} \quad (F.10) $$

and,

$$ \int_{\Omega_a} 2\mathbf{e}_{in} d\Omega = -\frac{4\pi}{15} \frac{5\lambda + 2}{\lambda + 1} [\nabla \mathbf{u} + (\nabla \mathbf{u})^{\dagger}] - \frac{8\pi}{25} \frac{1}{\lambda + 1} (\nabla \mathbf{b} - \frac{1}{3} \delta \nabla \cdot \mathbf{b}) \quad (F.11) $$

$$ \int_{\Omega_a} 2\mathbf{r}\mathbf{e}_{in} d\Omega = \frac{2\pi}{5(\lambda+1)} (\delta\mathbf{u}_r + \delta\mathbf{u}_r^{\dagger}) - \frac{4\pi}{15(\lambda+1)} \mathbf{u}_r \delta - \frac{4\pi}{15(\lambda+1)} \frac{1}{\lambda+1} (\delta\mathbf{b} + \delta\mathbf{b}^{\dagger}) + \frac{8\pi}{45(\lambda+1)} \mathbf{b}\delta. \quad (F.12) $$

In these expressions we did not include the dependence on $\propto \nabla\nabla\mathbf{u}$ and $\propto \nabla\nabla\mathbf{b}$ in Eqs. (F.10) and (F.12) because they turn out to be of $O(a^2/L^2)$, hence negligible in the averaged equation. Nevertheless, they will be needed to compute the deformation of the droplets, hence we display here the higher order terms, namely

$$ \begin{aligned} \int_{\Gamma_a} (\mathbf{x}\mathbf{x})_{jk} (\sigma_{out} \cdot \mathbf{n})_l d\Gamma = & \frac{7\lambda^2 + 190\lambda + 88}{105(\lambda^2 + 5\lambda + 4)} [(\nabla\nabla\mathbf{u})_{ijk} + (\nabla\nabla\mathbf{u})_{ikj}] + \frac{119\lambda^2 + 190\lambda - 24}{105(\lambda^2 + 5\lambda + 4)} (\nabla\nabla\mathbf{u})_{jki} \\ & - \frac{7\lambda^2 + 80\lambda - 72}{105(\lambda^2 + 5\lambda + 4)} [(\delta\nabla^2\mathbf{u})_{ijk} + (\delta\nabla^2\mathbf{u})_{ikj}] + \frac{2(13\lambda - 4)}{21(\lambda^2 + 5\lambda + 4)} (\delta\nabla^2\mathbf{u})_{jki} \\ & \frac{76}{735(\lambda + 1)} (\nabla\nabla\mathbf{b})_{ijk} + \frac{218}{3675(\lambda + 1)} (\delta\nabla^2\mathbf{b})_{jki} - \frac{272}{3675(\lambda + 1)} [(\delta\nabla^2\mathbf{b})_{ijk} + (\delta\nabla^2\mathbf{b})_{ikj}], \end{aligned} \quad (F.13) $$

$$ \int_{\Gamma_a} (\mathbf{x})_k (\mathbf{e}_{in})_{ij} d\Gamma = - \frac{7\lambda^2 + 20\lambda + 3}{105(\lambda^2 + 5\lambda + 4)} [(\nabla\nabla\mathbf{u})_{kij} + (\nabla\nabla\mathbf{u})_{kji}] - \frac{8(2\lambda + 1)}{21(\lambda^2 + 5\lambda + 4)} (\nabla\nabla\mathbf{u})_{ijk} $$
$$ \begin{aligned} & -\frac{\lambda - 10}{21(\lambda^2 + 5\lambda + 4)}[(\delta\nabla^2\mathbf{u})_{jki} + (\delta\nabla^2\mathbf{u})_{ikj}] + \frac{2(3\lambda - 2)}{21(\lambda^2 + 5\lambda + 4)}(\delta\nabla^2\mathbf{u})_{ijk} \\ & -\frac{32}{735(\lambda + 1)}(\nabla\nabla\mathbf{b})_{ijk} + \frac{292}{11025(\lambda + 1)}(\delta\nabla^2\mathbf{b})_{ijk} - \frac{22}{1225(\lambda + 1)}[(\delta\nabla^2\mathbf{b})_{jki} + (\delta\nabla^2\mathbf{b})_{ikj}]. \end{aligned} \quad (F.14) $$

### Deformation of the droplets

In the body of the text, we use the second moment of the momentum equation to determine the third mode of deformation $H_p^{(3)}$ and show that it is negligible. Here, we provide additional information on the derivation of the second moment of the momentum equation. Based on the general formula Eq. (C.6) one obtain by neglecting the inertial effects as well as the body forces,

$$ \int_{\Omega_a} r_j(\sigma_d^0)_{ik} + r_k(\sigma_d^0)_{ji} d\Omega + \oint_{\Gamma_a} r_j(\sigma_R^0)_{ik} + r_k(\sigma_R^0)_{ji} d\Gamma = \oint_{\Gamma_a} r_k r_j(\sigma_f^0 \cdot \mathbf{n})_i d\Gamma \quad (F.15) $$

Using the approach outlined by Lhuillier (1996), we can rewrite this equation by summing and subtracting the permutations: $B_{ijk} + B_{jik} - B_{kij}$, where $B_{ijk}$ corresponds to Eq. (F.15). This yields:

$$ \int_{\Omega_a} 2r_k(\sigma_d^0)_{ij} d\Omega + \oint_{\Gamma_a} 2\gamma r_k(\delta - \mathbf{n}\mathbf{n})_{ij} d\Gamma = + \oint_{\Gamma_a} r_k r_j(\sigma_f^0 \cdot \mathbf{n})_i + r_k r_i(\sigma_f^0 \cdot \mathbf{n})_j - r_j r_i(\sigma_f^0 \cdot \mathbf{n})_k d\Gamma \quad (F.16) $$

Taking the traceless part of this equation on the index $ij$, and substituting $\sigma_f^0$ by $\sigma_f^* + \Sigma$ and $e_d^0$ by $e_d^* + E$ one obtains,

$$ \begin{aligned} \oint_{\Gamma_a} 2\gamma r_k(\delta/3 - \mathbf{n}\mathbf{n})_{ij} d\Gamma = &+ \oint_{\Gamma_a} r_k r_j(\sigma_f^* \cdot \mathbf{n})_i + r_k r_i(\sigma_f^* \cdot \mathbf{n})_j - r_j r_i(\sigma_f^* \cdot \mathbf{n})_k d\Gamma \\ &- \frac{1}{3}\delta_{ij} \oint_{\Gamma_a} r_k r_l(\sigma_f^* \cdot \mathbf{n})_l + r_k r_l(\sigma_f^* \cdot \mathbf{n})_l - r_l r_i(\sigma_f^* \cdot \mathbf{n})_k d\Gamma \\ &+ \mu_f(1-\lambda) \int_{\Omega_a} 4\mathbf{E}_{ij} r_k d\Omega - \mu_f \lambda \int_{\Omega_a} 2r_k(\mathbf{e}_d^*)_{ij} d\Omega, \end{aligned} \quad (F.17) $$

where we have noticed that $\nabla \cdot \Sigma = -\rho_f g$ at $O(1)$ in $\phi$, and because we neglected the body forces in this example $\nabla \cdot \Sigma = 0$. Injecting Eqs. (F.13) and (F.14) in Eq. (F.17) one arrive at Eq. (6.35).

### Data availability

No data was used for the research described in the article.

### References

Advani, S.G., Tucker, C.L., 1987. The use of tensors to describe and predict fiber orientation in short fiber composites. J. Rheol. 31 (8), 751-784.

Appel, W., 2007. Mathematics for Physics and Physicists. Princeton University Press.

Aris, R., 2012. Vectors, Tensors and the Basic Equations of Fluid Mechanics. Courier Corporation.

Batchelor, G., 1970. The stress system in a suspension of force-free particles. J. Fluid Mech. 41 (3), 545-570.

Bothe, D., 2022. Sharp-interface continuum thermodynamics of multicomponent fluid systems with interfacial mass. Internat. J. Engrg. Sci. 179, 103731.

Buyevich, Y.A., Shchelchkova, I., 1979. Flow of dense suspensions. Prog. Aerosp. Sci. 18, 121-150.

Cartellier, A., Andreotti, M., Sechet, P., 2009. Induced agitation in homogeneous bubbly flows at moderate particle Reynolds number. Phys. Rev. E 80 (6), 065301.

Chapman, S., Cowling, T.G., 1990. The Mathematical Theory of Non-Uniform Gases: an Account of the Kinetic Theory of Viscosity, Thermal Conduction and Diffusion in Gases. Cambridge University Press.

Curtiss, C., 1956. Kinetic theory of nonspherical molecules. J. Chem. Phys. 24 (2), 225-241.

Drew, D.A., 1983. Mathematical modeling of two-phase flow. Annu. Rev. Fluid Mech. 15 (1), 261-291.

Drew, D.A., 1990. Evolution of geometric statistics. SIAM J. Appl. Math. 50 (3), 649-666.

Fintzi, N., 2025. Statistical Modeling of Disperse Two-Phase Flows: Application to Buoyancy-Driven Droplets Suspensions (Ph.D. thesis). Sorbonne Universite.

Guazzelli, É., Hinch, J., 2011. Fluctuations and instability in sedimentation. Annu. Rev. Fluid Mech. 43 (1), 97-116.

Guazzelli, E., Morris, J.F., 2011. A Physical Introduction to Suspension Dynamics. Vol. 45, Cambridge University Press.

Guazzelli, É., Pouliquen, O., 2018. Rheology of dense granular suspensions. J. Fluid Mech. 852, P1.

Haber, S., Hetronsi, G., 1971. The dynamics of a deformable drop suspended in an unbounded Stokes flow. J. Fluid Mech. 49 (2), 257-277.

Happel, J., Brenner, H., 2012. Low Reynolds Number Hydrodynamics: with Special Applications to Particulate Media. Vol. 1, Springer Science & Business Media.

Hesla, T.I., Huang, A.Y., Joseph, D.D., 1993. A note on the net force and moment on a drop due to surface forces. J. Colloid Interface Sci. 158 (1), 255-257.

Hill, R.J., Koch, D.L., Ladd, A.J., 2001. The first effects of fluid inertia on flows in ordered and random arrays of spheres. J. Fluid Mech. 448, 213-241.

Hinch, E., 1977. An averaged-equation approach to particle interactions in a fluid suspension. J. Fluid Mech. 83 (4), 695-720.

Hu, X., Ilgun, A.D., Passalacqua, A., Fox, R.O., Bertola, F., Milosevic, M., Visscher, F., 2021. CFD simulations of stirred-tank reactors for gas-liquid and gas-liquid-solid systems using OpenFOAM®. Int. J. Chem. React. Eng. 19 (2), 193-207.

Hwang, G., Shen, H., 1989. Modeling the solid phase stress in a fluid-solid mixture. Int. J. Multiph. Flow 15 (2), 257-268.

Ishii, M., Hibiki, T., 2010. Thermo-Fluid Dynamics of Two-Phase Flow. Springer Science & Business Media.

Jackson, R., 1997. Locally averaged equations of motion for a mixture of identical spherical particles and a Newtonian fluid. Chem. Eng. Sci. 52 (15), 2457-2469.

Jackson, R., 2000. The Dynamics of Fluidized Particles. Cambridge University Press.

Joseph, D., Lundgren, T., Jackson, R., Saville, D., 1990. Ensemble averaged and mixture theory equations for incompressible fluid-particle suspensions. Int. J. Multiph. Flow 16 (1), 35-42.

Junqua-Moullet, A., 2003. Détermination Expérimentale et Modélisation de la Concentration D'aire Interfaciale en Écoulement Stratifié Horizontal (Ph.D. thesis). Grenoble INPG.

Kataoka, I., 1986. Local instant formulation of two-phase flow. Int. J. Multiph. Flow 12 (5), 745-758.

Kim, S., Karrila, S.J., 2013. Microhydrodynamics: Principles and Selected Applications. Courier Corporation.

Kim, S., Russel, W.B., 1985. Modelling of porous media by renormalization of the Stokes equations. J. Fluid Mech. 154, 269-286.

Leal, L.G., 2007. Advanced Transport Phenomena: Fluid Mechanics and Convective Transport Processes. Vol. 7, Cambridge University Press.

Lhuillier, D., 1992. Ensemble averaging in slightly non-uniform suspensions. Eur. J. Mech. B Fluids 11 (6), 649-661.

Lhuillier, D., 1996. Contribution of the particulate phase to the overall stress of a suspension. Comptes Rendus de L'Academie Des Sciences-Serie IIb-Mecanique Phys. Chim. Astron. 323 (1), 3-8.

Lhuillier, D., 2022. Revisiting the two-phase modelling of non-Brownian Stokesian suspensions. Unpublished note.

Lhuillier, D., Chang, C.-H., Theofanous, T.G., 2013. On the quest for a hyperbolic effective-field model of disperse flows. J. Fluid Mech. 731, 184-194.
Lhuillier, D., Morel, C., Delhaye, J.-M., 2000. Bilan d'aire interfaciale dans un mélange diphasique: approche locale vs approche particulaire. Comptes Rendus de L'Académie Des Sciences-Series IIB-Mechanics-Physics-Astronomy 328 (2), 143-149.

Lhuillier, D., Nadim, A., 1998. Fluid dynamics of particulate suspensions: selected topics. In: Continuum Models and Discrete Systems - Proceedings of the 9th International Symposium. World Scientific, pp. 180-197.

Lhuillier, D., Nadim, A., 2009. Rheology of suspensions of mass-polarized particles in a gravitational field. Chem. Eng. Commun. 197 (1), 76-91.

Lhuillier, D., Nozières, P., 1992. Volume averaging of slightly non-homogeneous suspensions. Phys. A 181 (3-4), 427-440.

Lhuillier, D., Theofanous, T.G., Liou, M.-S., 2010. Multiphase flows: compressible multi-hydrodynamics. In: Handbook of Nuclear Engineering.

Marle, C., 1982. On macroscopic equations governing multiphase flow with diffusion and chemical reactions in porous media. Internat. J. Engrg. Sci. 20 (5), 643-662.

Morel, C., 2015. Mathematical Modeling of Disperse Two-Phase Flows. Springer.

Nadim, A., 1996. A concise introduction to surface rheology with application to dilute emulsions of viscous drops. Chem. Eng. Commun. 148 (1), 391-407.

Nadim, A., Stone, H.A., 1991. The motion of small particles and droplets in quadratic flows. Stud. Appl. Math. 85 (1), 53-73.

Nott, P.R., Guazzelli, E., Pouliquen, O., 2011. The suspension balance model revisited. Phys. Fluids 23 (4), 043304.

Orlando, G., Barbante, P.F., Bonaventura, L., 2023. On the evolution equations of interfacial variables in two-phase flows. Int. J. Multiph. Flow 169, 104617.

Pähtz, T., Chen, Y., Zhu, R., Tholen, K., He, Z., 2025. General average balance equations for a mixture of rigid particles and a continuous fluid. arXiv preprint arXiv:2308.09661.

Panicker, N., Passalacqua, A., Fox, R.O., 2018. On the hyperbolicity of the two-fluid model for gas-liquid bubbly flows. Appl. Math. Model. 57, 432-447.

Pozrikidis, C., 2011. Introduction to Theoretical and Computational Fluid Dynamics. Oxford University Press.

Pozrikidis, C., et al., 1992. Boundary Integral and Singularity Methods for Linearized Viscous Flow. Cambridge University Press.

R. Gatignol, R., Prud'homme, R., 2001. Mechanical and Thermodynamical Modeling of Fluid Interfaces. Vol. 58, World Scientific.

Raja, R.V., Subramanian, G., Koch, D.L., 2010. Inertial effects on the rheology of a dilute emulsion. J. Fluid Mech. 646, 255-296.

Rallison, J., 1978. Note on the Faxén relations for a particle in Stokes flow. J. Fluid Mech. 88 (3), 529-533.

Rallison, J., 1984. The deformation of small viscous drops and bubbles in shear flows. Annu. Rev. Fluid Mech. 16 (1), 45-66.

Randolph, A., 2012. Theory of Particulate Processes: Analysis and Techniques of Continuous Crystallization. Elsevier.

Rao, K.K., Nott, P.R., Sundaresan, S., 2008. An Introduction to Granular Flow. Vol. 490, Cambridge University Press Cambridge.

Rubinow, S.I., Keller, J.B., 1961. The transverse force on a spinning sphere moving in a viscous fluid. J. Fluid Mech. 11 (3), 447-459.

Simonin, O., 1996. Continuum modelling of dispersed two-phase flows. Lect. Series-Van Kareman Inst. Fluid Dyn. 2, K1-K47.

Stone, H., Brady, J., Lovalenti, P., 2001. Inertial effects on the rheology of suspensions and on the motion of individual particles. Preprint.

Subramanian, R.S., 1985. The Stokes force on a droplet in an unbounded fluid medium due to capillary effects. J. Fluid Mech. 153, 389-400.

Taylor, G.I., 1932. The viscosity of a fluid containing small drops of another fluid. Proc. R. Soc. Lond. Ser. A Containing Pap. A Math. Phys. Character 138 (834), 41-48.

Taylor, T., Acrivos, A., 1964. On the deformation and drag of a falling viscous drop at low Reynolds number. J. Fluid Mech. 18 (3), 466-476.

Teigen, K.E., Li, X., Lowengrub, J., Wang, F., Voigt, A., 2009. A diffuse-interface approach for modeling transport, diffusion and adsorption/desorption of material quantities on a deformable interface. Commun. Math. Sci. 4 (7), 1009.

Tryggvason, G., Scardovelli, R., Zaleski, S., 2011. Direct Numerical Simulations of Gas-Liquid Multiphase Flows. Cambridge University Press.

Wegener, M., Paul, N., Kraume, M., 2014. Fluid dynamics and mass transfer at single droplets in liquid/liquid systems. Int. J. Heat Mass Transfer 71, 475-495.

Young, N., 1959. The motion of bubbles in a vertical temperature gradient. J. Fluid Mech. 6 (3), 350-356.

Zaepffel, D., 2011. Modélisation des Écoulements Bouillants à Bulles Polydispersées (Ph.D. thesis). Grenoble.

Zhang, D.Z., 2021. Ensemble average and nearest particle statistics in disperse multiphase flows. J. Fluid Mech. 910.

Zhang, D., Prosperetti, A., 1994a. Averaged equations for inviscid disperse two-phase flow. J. Fluid Mech. 267, 185-219.

Zhang, D., Prosperetti, A., 1994b. Ensemble phase-averaged equations for bubbly flows. Phys. Fluids 6 (9), 2956-2970.

Zhang, D., Prosperetti, A., 1997. Momentum and energy equations for disperse two-phase flows and their closure for dilute suspensions. Int. J. Multiph. Flow 23 (3), 425-453.

Zhang, D.Z., Wang, M., Balachandar, S., 2023. Evolution of the age-included nearest pair distribution in disperse multiphase flows. Phys. Fluids 35 (6).