Drag on a spherical particle at the air-liquid interface

We investigate the flow of viscous interfaces carrying an insoluble surface active material, using numerical methods to shed light on the complex interplay between Marangoni stresses, compressibility, and surface shear and dilatational viscosities. We find quantitative relations between the drag on a particle and interfacial properties as they are required in microrheology, i.e., going beyond the asymptotic limits. To this end, we move a spherical particle probe at constant tangential velocity, symmetrically immersed at either the incompressible or compressible interface, in the presence and absence of surfactants, for a wide range of system parameters. A full three-dimensional finite element calculation is used to reveal the intimate coupling between the bulk and interfacial flows and the subtle effects of the different physical effects on the mixed-type velocity field that affects the drag coefficient, both in the bulk and at the interface. For an inviscid interface, the directed motion of the particle leads to a gradient in the concentration of the surface active species, which in turn drives a Marangoni flow in the opposite direction, giving rise to a force exerted on the particle. We show that the drag coefficient at incompressible interfaces is independent of the origin of the incompressibility (dilatational viscosity, Marangoni effects or a combination of both) and that its higher value can not only be related to the Marangoni effects, as suggested earlier. In confined flows, we show how the interface shear viscosity suppresses the vortex at the interface, generates a uniform flow, and consequently increases the interface compressibility and the Marangoni force on the particle. We mention available experimental data and provide analytical approximations for the drag coefficient that can be used to extract surface viscosities. VC 2021 Author(s). All article content, except where otherwise noted, is licensed under a Creative Commons Attribution (CC BY) license (http:// creativecommons.org/licenses/by/4.0/). https://doi.org/10.1063/5.0050936


I. INTRODUCTION
The rheology of complex fluid-fluid interfaces is a relevant and rich physical problem. 1,2 For example, the seminal theoretical work of Saffman and Delbr€ uck 3 on the Brownian motion in biological membranes revealed how the momentum diffusion at the interface and in the bulk are strongly coupled and, for non-compressible interfaces, the ratio between bulk and surface viscosities determines the range of the hydrodynamic effects at the interface. However, when Brownian motion is used to more generally deduce interfacial rheological properties beyond lipid membranes, values obtained from microrheology have been found to be very different from those with macroscopic interfacial shear rheometers. 4 Several factors have been identified to contribute to the drag of particles close to and at an interface, and possible sources of discrepancies have been identified, with some being related to dissipation modes of the surface or to electrokinetic effects 5 and others focused on the non-correct analysis of the hydrodynamic conditions of the measurements, 6 although for direct shear rheological techniques such hydrodynamic analysis has led to clear operating windows. 7,8 One aspect which merits further investigation is the role of interfacial compressibility and the possible influence of Marangoni flows on the interfacial and bulk velocity fields. First, there is clear experimental evidence from studying needles of different aspect ratios that compression might be of relevance in interpreting microrheological experiments. 9 Second, the analysis of different types of interfacial rheometers, which use mixed interfacial flow fields for measuring the surface rheological properties in mixed shear-dilation/compression 10 or shear-extension-dilation/compression, 11 has already revealed a complex and subtle interplay between the Marangoni stresses, bulk flow effects, and the interfacial stresses stemming from the desired rheological material functions. In the present work, we investigate to what extent this analysis holds true for the flow field caused by a particle dragged through the interface (the disturbance velocity field), and resolving the effects of the various coupled contributions to the forces and concentration fields, as such a translating probe produces a mixed-type flow field. 12 In particular, we are going to investigate the role of interfacial compressibility and the interplay between the different interfacial deformation modes (shear and local dilatation/ compression).
The drag of an object at an interface has already received attention previously, and typically, two effects have been studied, i.e., the effect of coupling of the interfacial momentum diffusion to that in the bulk and the effect of Marangoni contributions. The original work of Saffman and Delbr€ uck 3 examined a cylindrical probe translating at a membrane in the limit of large membrane viscosity and an incompressible interfacial layer. In this limit, a dependence of drag on the particle size was found to be essentially absent (only a weak logarithmic effect had been predicted), and the range of the hydrodynamic forces was identified as the Saffman-Delbr€ uck length. Later, Hughes et al. 13 obtained the exact solution to the equations of motion for any combination of membrane and adjacent liquid viscosities. In both these works, the domain of the bulk flow was considered infinite. The influence of the bulk layer thickness was studied by Evans and Sackmann. 14 They imposed a proportionality between the velocity and the shear stress exerted by the subphase on the membrane. This work has been generalized by Stone and Ajdari 15 to any value of subphase depth using numerical calculations. Barentin et al. 16,17 used the same geometry as Evans and Sackmann 14 and studied both cases of an incompressible Langmuir (insoluble) monolayer and a Gibbs (soluble) monolayer. They exploited the lubrication approximation for shallow films and obtained analytical expressions for the shear stress and the drag force. In all studies mentioned so far, the particle was cylindrical and non-protruding. For a protruding spherical particle trapped within a very thin, Newtonian fluid-fluid interface, the drag has been calculated only numerically. While Danov et al. 18,19 considered a compressible interface characterized by a Boussinesq-Scriven model with shear and dilational viscosities, Fischer 20 argued that due to the short time scales on which a sustained gradient in surface tension gives rise to Marangoni contributions to the forces, the monolayer can under certain conditions be considered incompressible, as the presence of a surfactant strongly suppresses any motion at the surface that compresses or expands the interface. Manikantan and Squires 21 provided an exhaustive, excellent discussion of these phenomena. In the study of monolayers of fatty acids and phospholipids, the surface shear viscosity calculated following Danov's approach was shown to underestimate the experimental value. 22 Fischer et al. 23 then solved the Stokes flow equations for a sphere immersed in a viscous interface (monolayer) upon incorporating the Marangoni effect into the treatment, by solving the equations for an incompressible interface. Fischer showed that the Marangoni effect contributes a significant part of the total drag at the limit of vanishing surface compressibility. However, not all systems will adhere to the conditions required for a vanishing surface compressibility, and how this will come into play is ill-understood, in particular when surface viscosities are significant as well. A recent approach by Elfring et al. 12 to characterize the latter used the lubrication approximation for a very shallow subphase and considered a non-protruding disk-shaped probe at a viscous, compressible interface that generates Marangoni flows. Elfring et al. assumed the presence of a soluble surfactant that equilibrates on a finite timescale and did not consider transient effects. While the lubrication approximation allows for analytical solution through a perturbation expansion, using a shallow subphase is not ideal for probing rheological properties of the interface, because in this regime, bulk stresses dominate over surface stresses. 12 The present work can be considered as an extension of the work of Elfring et al. to a protruding particle and a deep subphase. We study the transient Marangoni flow and the interface compressibility of a Newtonian interface endowed with both shear and dilatational viscosity. The interface is covered by an insoluble surface active component (surfactant, protein, fatty acid, lipid, or polymeric substance) and a large spherical particle, which is set to translate at the interface. The interfacial contribution to the drag on the particle is a combination of the Marangoni effect and the interfacial stresses arising from the interfacial strain fields and relevant interfacial viscosities. This interfacial disturbance velocity field depends on both the balance between the Marangoni flows and the shear and local dilational/compressional components and how the momentum sink of the bulk influences this. To the best of our knowledge, this is the first study on a protruding three-dimensional particle at compressible interfaces with Marangoni flows. Dimova et al. 24 also considered concentration variations with a spherical particle, but did not couple the interface to a bulk fluid and did not determine the effects of Marangoni flow.
We begin by recalling the hydrodynamic equations and formulate the boundary and initial conditions for the problem at hand in Sec. II. We then demonstrate how we calculate the drag on the particle and introduce dimensionless quantities to come up with reduced equations and a number of dimensionless numbers characterizing all parameters of the particle, the interface, and the bulk phases. The numerical implementation is described in Sec. III. In Sec. IV, we present results for various regimes, the surfactant dominated regime (Sec. IV A), a shear-free interface without (Sec. IV B), and with Marangoni flow (Sec. IV C). In Sec. IV D, we focus on the flow field of regime of finite shear viscosity, as well as a reduced parameter space that captures the situation encountered in real surfactant systems. The effects of interface viscosites on the drag coefficient are discussed in Sec. IV E. Conclusions are provided in Sec. V, with a caveat regarding the use of these drag flows for determining shear viscosities or rheological properties.

II. HYDRODYNAMIC EQUATIONS
Consider the rather general situation of a spherical particle symmetrically immersed at a planar interface, with its equator residing in the interface plane, and translating with a constant velocity U ¼ Ue x (Fig. 1). The two fluids in domains X 1 and X 2 are considered incompressible and Newtonian, i.e., modeled by where r denotes position, p i is the stress tensor in domain i, p is the pressure field, I is the unity tensor, s i ¼ g i ½$u þ ð$uÞ T is the viscous stress tensor, u is the fluid velocity, and g i is the viscosity of fluid i 2 f1; 2g. All fields appearing in our equations, such as u and s i , are functions of spatial coordinate r and time t, and we skip these arguments for better readability. A similar decomposition can also be written for the interface stress tensor 1,25 where s s is the extra surface stress tensor, I s ¼ I À nn is the surface or tangential projection tensor 25,26 with surface normal n, and c is the surface tension of the interface. For an isothermal system, the surface tension is solely a function of the surfactant concentration C, i.e., c ¼ cðCÞ. Its functional form remains to be specified (Sec. IV). The extra surface stress tensor is described by the Boussinesq-Scriven constitutive law 27,28 where g s and j s are the shear and dilatational viscosities of the interface, respectively. The surface viscosities can be made functions of the surface pressure P s ðCÞ ¼ c 0 À cðCÞ, where c 0 is the surface tension of a clean interface (C ¼ 0). It introduces additional complexities, such as irreversible particle motion. 29 However, for simplicity in the current analysis, they are here assumed to be system parameters, and an explicit dependence on C is not included. 28,30 In most practical cases, one needs some surface active species at the interface to get significant extra (viscous) stresses, 31,32 but surface viscosities for clean interfaces have also been reported. 33 The surface rate of deformation tensor D s is defined as where u s is the velocity u evaluated at the interface, and $ s ¼ I s Á $ is the surface gradient operator. 26 Note that this is a pseudo-definition of $ s , as $ does not exist for a surface field. Given the small length scale of the particle involved, the dynamics of the incompressible bulk fluids can be modeled with the creeping-flow (or Stokes) equations, which reduce the momentum balance to $ Á p i ¼ 0. With the help of Eq. (1), this becomes where the vanishing divergence of u expresses the balance of mass assuming a constant density within the domains. The velocity and pressure fields thus receive their time-dependency through C. The evolution of the surfactant concentration C is governed by the unsteady surface convection-diffusion (SCD) equation 26,[34][35][36] where D s is the surface diffusivity of the surfactant and where we assumed a planar interface. Equations (5) and (6) are the governing equations for u and C as a function of time and r. They require proper boundary and initial conditions to be solved. In this study, it is assumed that the interface remains planar, which is a valid assumption for the typical particle sizes and interfacial tensions used in experiments. 37 Moreover, we solve the hydrodynamic problem only for the translational motion of the spherical colloids along the interface, i.e., under the assumption of contact line pinning, which prevents the particle from rotating. It has been shown experimentally that the presence of surface roughness leads to pinning of the contact line. 38 The pinned contact line is also used for modeling single 37 and pair 39 particles at the interfaces. Due to pinning of the contact line, an extremely small deformation of the interface is needed to balance the torque acting on the particle, which justifies our approach of a non-rotating particle at a planar interface.

A. Boundary and initial conditions
For the surfactant concentration C governed by Eq. (6), we use the following initial and boundary conditions: where S I is the fluid-fluid interface, @S p is the circular particleinterface boundary, @S b is the square box boundary, and n p is the unit vector normal to the particle surface and tangential to the interface. We thus assume an initially undisturbed, uniform C profile at the fluid-fluid interface, which remains at the initial value at all times at the box boundary. The third boundary condition means no surfactant flux from the particle interface boundary. The boundary conditions for the velocity field u are composed of boundary conditions on the surface S p of the particle, at the fluid-fluid interface S I , and far away from the particle on the simulation box surface S b . In this study, we assume that there is no-slip at the particlefluid interface as well as at the interface between two Newtonian fluids, FIG. 1. Schematic representation of the system under study. A spherical particle of radius R translates with constant velocity U ¼ Ue x at the flat interface (y ¼ 0) between two Newtonian fluids, in domains X 1 and X 2 . The normal unit vectors at the particle surface S p and fluid-fluid interface are denoted by n p and n ¼ e y , respectively. The fluid-fluid interface between the two domains X 1 and X 2 we denote by S I , the circular intersection of the particle and the fluid-fluid interface is @S p . The interface boundary has contributions from @S p and the square intersection @S b (large square) between S I and the simulation box boundary enclosing X 1 and X 2 . For an unbounded domain, the interfacial area is infinitely large, but for the numerical investigation and to study the effect of subphase depth we consider a large but finite, cubic box of size L ¼ L x ¼ L y ¼ L z (a snapshot of the real system will be provided in Sec. III). The particle is symmetrically placed at the fluid-fluid interface and initially located at the box center (if not otherwise mentioned).

ARTICLE
scitation.org/journal/phf while more generally, there can be relative motion between surfaces in contact. 25,40,41 However, the slip boundary condition at the interface is usually only employed for macromolecular, non-Newtonian fluids, and its necessity for Newtonian fluids remains unclear. 40 We assume that the particle translates with velocity U ¼ _ r p , where r p is the location of the center-point of the particle, and the particle does not rotate. 42 Applying no-slip on the particle surface S p and vanishing velocity on the simulation box surface S b yields Boundary conditions for u at the fluid-fluid interface are defined by jump equations. In the absence of mass transfer (phase change) across the interface, there are the mass jump balance 25 and the momentum jump balance n Á ðs 2 À s 1 Þ À ðp 2 À p 1 Þn ¼ À$ s Á s s À $ s c þ cð$ s Á nÞn; (10) where n ¼ e y is interface surface normal pointing to region X 2 ( Fig. 1), u i and p i denote the velocity and pressure of fluid i at r 2 S I , and $ s Á n is the interface mean curvature. The subscripts 1 and 2 indicate the value of quantity in the bulk phases extrapolated to the interface. The momentum jump balance can be split into normal, n Á ðs 2 À s 1 Þ Á n À ðp 2 À p 1 Þ ¼ Àð$ s Á s s Þ Á n þ cð$ s Á nÞ; (11) and tangential components, where t is a unit tangent vector residing in the x-z-plane (Fig. 1). In our setup, the interface remains at y ¼ 0, coinciding with the y-coordinate of center of the sphere, the normal component of u s therefore vanishes, and the momentum jump balance in normal direction is replaced by u s Á n ¼ 0. For the tangential component of u s , we assume a no-slip boundary conditions at the interface, which implies continuity of the tangential component of the velocity across the interface, Finally, we note that the tangential momentum jump balance can be written in in terms of the surface pressure 25 The Gibbs-Marangoni modulus K p ¼ C@P s =@C allows one to relate the gradient in surface pressure to the gradient in concentration, 21,26 and to then rewrite the tangent component of the stress jump equation in the following final form: B. Drag on the particle The general relation between the drag on the particle F at the interface translating with a constant velocity U is complex and depends on the particle geometry, bulk and interface rheological properties, the interfacial equation of state, and the transport properties of the surfactant. The drag coefficient on a spherical particle with radius R in an unbounded Newtonian fluid translating with a constant velocity is the Stokes' drag coefficient f ¼ jFj=jUj ¼ 6pgR. 15 For a spherical particle symmetrically embedded in a clean, inviscid interface between a Newtonian and an inviscid fluid, the drag is exactly half the Stokes' drag, i.e., f ¼ 3pgR. 43 The drag force on the spherical particle embedded at the interface is the sum of interface and bulk contributions, where @S p is the perimeter of the particle at the interface and n p is the unit normal vector to the surface of the particle (Fig. 1). The first integral on the right hand side of (17) is the bulk force F b . The second integral is the interface force F i and can be decomposed into two con- The first integral on the right hand side of Eq. (18) is the elastic or Marangoni part of the drag force F M due to the non-uniform distribution of the surfactant in the contact line. 24 The second integral is the interface viscous force F s due to the extra interface stress tensor. Hence, the total drag force (17) on the particle is Because our particle is non-rotating, a torque will act on the particle, which can be evaluated via suitable adapted integrals involving the relative position x ¼ r À r p , where r p denotes the position of the particle's center. As explained earlier, it is assumed that this torque is balanced by an extremely small deformation of the interface.

C. Dimensionless quantities and equations
For the problem at hand, we can use the constant velocity U ¼ jUj, the radius R of the spherical particle, the viscosity g 1 of fluid 1, and the initial surfactant concentration C 0 to introduce reduced units, and to come up with a number of dimensionless parameters. Dimensionless variables, marked by asterisk (only here), are therefore defined uniquely in terms of their dimensional counterparts, such as From now on, and for the rest of this work, variables are meant to be dimensionless, and we omit the asterisk for brevity. Using them in the above sections, Eq. (5) is replaced by Physics of Fluids ARTICLE scitation.org/journal/phf where k ¼ g 2 =g 1 is the ratio between the two viscosities. The normal and tangential parts of the jump balance, Eqs. (11) and (16) remain formally unaltered, i.e., with s 1 ¼ 2D and s 2 ¼ 2kD. To close the system of equations, and to fully specify the dimensionless quantities, an interfacial equation of state must be chosen, which will be done in Sec. IV. Because we are going to assume (i) a linearized equation of state there, K p will be replaced by Ma C, where Ma is a dimensionless Marangoni number, (ii) a flat interface, the last term in the first line of Eq. (22) will disappear, and (iii) an air-liquid interface, k ¼ 0.
The dimensionless extra surface stress tensor (3) is written in terms of the Bousinessq number Bq 1 ¼ g s =Rg 1 , the ratio of surface shear viscosity to bulk viscosity, and the ratio Bq 2 ¼ j s =Rg 1 between surface dilation viscosity to bulk viscosity. Equation (6) becomes introducing the surface P eclet number Pe s ¼ UR=D s , a ratio of diffusion time, R 2 =D s , to the characteristic time for particle motion, R/U. The force expression (17) remains unchanged, and the integral over S p is split into parts that arise from the hemispheres in the two fluids with their corresponding p 1 and p 2 . The constitutive Eq. (1) as well as the mass jump balance (9) and boundary condition (13) remains unchanged in their dimensionless form, while the initial and boundary conditions (7) and (8) for the surfactant concentration and velocity field become Cðr; t ¼ 0Þ ¼ 1 for r 2 S I ; Cðr; tÞ ¼ 1 for r 2 @S b ; n Á ½$ s Cðr; tÞ ¼ 0 for r 2 @S p ; u ¼ e x for r 2 S p , and

III. NUMERICAL METHOD
The problem is composed of two sub-problems: the Stokes equation in the fluid domain and the SCD equation at the interface. These two sub-problems are coupled and must be solved together. We employed the finite element method (FEM) in a domain of finite size. To avoid tracking the moving particle, the equations are solved in a frame that is moving with the particle, i.e., the particle is kept stationary at x ¼ ðx; y; zÞ ¼ 0, and the walls move with velocity ÀUe x . Note that all our results are presented in a frame where the particle is moving and the walls are stationary. The Stokes equation (21) is thus solved for with a vanishing velocity on the particle surface @S p , the dimensionless boundary condition u ¼ Àe x on the outer surfaces of the simulation box, together with the stress boundary condition (22) at the fluid-fluid interface, using constitutive equation (23). The solution to the Stokes problem yields the velocity and pressure fields u i ðxÞ and p i ðxÞ, hence s i ðxÞ. To include the effect of the surfactant, in this work two approaches are employed: an implicit approach and an explicit approach. In the implicit approach, the SCD is solved together with the Stokes equation in one system. The non-linearity in the SCD is handled by a Picard iteration. In the explicit approach, the velocity at the fluid-fluid interface is used to solve the SCD Eq. (24) with boundary and initial conditions (7). After solving this equation, the new C is used to update interface tension c using the interfacial equation of state, to be specified in Eq. (25). The explicit approach typically works well for problems at low Ma, i.e., for low interfacial elasticity. However, at larger values of Ma, the time step must be chosen extremely small for stable simulations. For those cases, the implicit approach is used, which is more expensive per time step, but which allows for an arbitrarily large time step. Steady-state results are obtained by running the simulation until the dimensionless time t ¼ 100, which was found to be sufficiently long (see supplementary material Sec. S4).
We used an in-house FEM code to solve the full problem, i.e., the complete set of equations. A validation and convergence study is provided in the supplementary material Secs. S1(ii) and S1(iii). We make use of the open-source mesh generator Gmsh, 44 which offers great control over the local element size in the domain (Fig. 2). Due to symmetry in the xy-plane, only half of the domain is meshed, and appropriate symmetry conditions are employed. In order to diminish boundary effects, a large simulation box is used: L ¼ 3200 if not otherwise mentioned, which was shown to be sufficiently large (see supplementary material Sec. S2). The spherical particle is symmetrically immersed at the interface and translates with a constant velocity in the x-direction. The interface is located at y ¼ 0, and it divides the simulation box into two domains: domain X 1 in y < 0 is filled by a viscous liquid and X 2 in y > 0 is filled by an inviscid fluid. The Stokes equation is solved in X 1 . Details of the numerical implementation are discussed by Dietrich et al. 37 and Carrozza et al. 45 We tested mesh convergence and used our mesh M3 as a compromise between quality and efficiency if not mentioned otherwise [supplementary material Sec. S1(iii)]. To solve Eq. (24), second-order finite difference schemes are used to approximate the time derivatives.

IV. RESULTS AND DISCUSSION
In all the results presented in this work, we assume a linear interfacial equation of state: c ¼ c 0 À Ck B T (dimensional form) where c 0 is the surface tension of a clean interface, and thus P s ¼ K p ¼ Ck B T (dimensional form). Making the equation of state dimensionless using Eq. (20), we obtain Physics of Fluids ARTICLE scitation.org/journal/phf giving rise to the dimensionless capillary number Ca ¼ g 1 U=c 0 (ratio between bulk stress and interface tension force) and Marangoni number Ma ¼ C 0 k B T=g 1 U, 46 the ratio between the force tending to deform the interface and surface elasticity, which tends to restore the original shape of the film. Note that from this we obtain P s ¼ K p ¼ Ma C in our dimensionless form, which shows that the Marangoni contribution to the force is proportional to Ma. Assuming a linear interfacial equation of state is a strong assumption, employed in most of previous works, and to be discussed later below. The presented equations and methods capture the case of curved fluid-fluid interfaces. However, for this study, phase 2 is considered inviscid, therefore g 2 ¼ 0; s 2 ¼ 0, and the dimensionless number k vanishes. Moreover, the interface is assumed to remain straight, i.e., the capillary number Ca ( 1 and the mean curvature vanishes, which is a valid assumption for the typical particle sizes used in experiments. 37 As the last term in Eq. (22) vanishes for a flat interface, results do not depend on the choice of Ca. Due to this assumed linear relationship between c and C, we are left with dimensionless parameters Bq 1 ; Bq 2 , Ma, and Pe s , which completely describe the problem.
We begin by investigating the situation in the absence of any extra surface stresses, Bq 1 ¼ Bq 2 ¼ 0, in Sec. IV A. Since one of our main objectives is to distinguish between the Marangoni and dilatational (Bq 2 ) effects on the interface compressibility, we are going to introduce them separately here in Sec. IV B and to then study the interface when both mechanisms are present (Sec. IV C). It is worth noting at the outset that the case of vanishing Bq 2 ¼ 0 is an artificial case that has no analogue in the world of real interfaces. It is however important to delineate the contributions from shear and dilation to the measurable quantities and to allow for a comparison with previous results by Elfring et al. 12 The ratio which describes the viscous resistance to dilation, has been used 9,11 instead of Bq 2 to present and discuss results, as this ratio controls which deformation mode is "cheaper," i.e., shear or dilation/compression. An analogy can be made to linear elastic materials, where the Poisson ratio can be expressed in terms of the shear and bulk (compression) modulus. Whenever we fix Bq 1 and vary Bq 2 or vice versa, our graphs describe the effect of H as well. The interface will be considered to exhibit a non-vanishing shear viscosity in Sec. IV D, while in Sec. IV E we focus on the influences of both interface viscosities on the drag coefficient. To guide a reader, the structure of the manuscript is further visualized in Table I.

A. Surfactant dominant regime
While the drag coefficient on the particle had been studied before by Avramov et al. 47 in the absence of extra surface stresses, we here present the most relevant results (including interface compressibility aspects) for this special case, as they serve as a benchmark for the results in Secs. IV B-IV E, and to validate our implementation.
When a particle translates at the interface covered with a surfactant, it pushes the surfactant, resulting in an accumulation of surfactant in the front of the particle, and a depletion at its rear. The surfactant concentration around a translating particle, for Pe s ¼ 1, is shown in Fig. 3(a). This difference in the surfactant concentration causes a Marangoni flow; a flow from high to low concentration. 48 Surfactant concentration profiles in the x-direction along a line at z ¼ 0, passing through the center of the spherical particle at different Ma are shown in Fig. 3(b). To quantify the surfactant variation, we Physics of Fluids ARTICLE scitation.org/journal/phf calculate the difference in the maximum and minimum concentrations, DC. The inset in Fig. 3(b) shows DC at different Ma indicating that the relative surfactant variation is high for low Ma and decreases with increasing Ma. The Marangoni flow is also highly affected by the Pe s (supplementary material Fig. S11). At very small Pe s ( 1, the time required for recovering homogeneity of the surfactant concentration due to fast diffusion of the surfactants relative to the particle motion is very short; therefore, uniformity is obtained instantaneously, t ( 1. At Pe s < 0:02, there are negligible surfactant variations at the interface and at Pe s ¼ 0:001 the surfactant concentration remains uniform. Beyond Pe s > 1, the variation does not change much upon further increasing Pe s .
The variation in surfactant concentration, in turn, creates a gradient in $c or $P s , which acts against the surface compression in front of, and dilation at the rear of the particle. 49,50 Analogous to the bulk compressibility, the compressibility of monolayers can be defined as C ¼ ÀCdA=dP s j T , where A $ C À1 is the area per surfactant molecule at S I .
The interface symmetrically compresses at the front of the particle and dilates at the rear. The maximum dilation at the rear of the interface, denoted as ð$ s Á u s Þ max , serves as an indicator for the maximum interface compressibility. Barentin et al. 16 argued that the monolayer can be considered incompressible if j$ s Á u s j ( 1, as one can write the incompressibility condition in the bulk close to a planar interface as Provided that @u y =@y ¼ 0, the compression $ s Á u near the interface vanishes. To further demonstrate the effect of coupling to the subphase on the compressibility of the clean interface, we varied the depth of the subphase. For a very shallow subphase there is a "layered" structure in the subphase (supplementary material Fig. S7), the normal velocity basically vanishes, and the interface becomes less compressible [supplementary material Fig. S8(a)]. Moreover, as the subphase depth is decreased, the bulk stresses, and thus the drag coefficient, diverge [supplementary material Fig. S8(b)]. These effects can be logically attributed to an increased importance of the subphase for decreasing subphase depth, inherently limiting the sensitivity of the probe to properties of the interface, as also concluded by Elfring et al. 12 For a deep subphase, two different mechanisms can make an interface effectively incompressible: Marangoni effects and a viscous resistance to compressional flow. Figure 4(a) shows the effects of Pe s on the interface compressibility at Ma ¼ 10. Maximum interface dilation ð$ s Á u s Þ max for Ma ¼ 10 and Ma ¼ 20 as a function of Pe s are shown in Fig. 4(b). These results show that at the limits of small and high Pe s , interface compression is independent of Ma. At Pe s % Oð1Þ, an increase in Ma from 10 to 20 reduces the interface compressibility. These results show that when there is an increase in Ma in order to keep the interface compression constant, we should decrease Pe s . Following the reasoning of Chisholm et al., we can use the following simple argument to justify these results. 51 Exploiting the interfacial equation of state, we can rewrite the SCD in the following form: where the steady state is considered for simplicity. This equation shows that at high Ma the convection term is less important. The diffusion term at the right hand side also vanishes when ðMa Pe s Þ À1 ( 1, which shows when there is an increase in Ma, in order to recover the interface compression Pe s should decrease. An increase in Ma also increases the Marangoni contribution to the force on the particle (in the regime where the linear equation of state holds), which consequently increases the drag coefficient, see Figs. 4(c) and 4(d). This is a known result first observed by Plateau, 52 by comparing the damping of a magnetic needle immersed in a liquid with one placed at the liquid surface. He found that the resistance against deformation was higher at the surface than in the bulk and attributed this to a non-vanishing extra surface shear viscosity. Marangoni realized that trace contaminations are almost always present at liquid interfaces and that a non-uniform distribution of surface-active material, together with an interfacial equation of state, is responsible for Plateau's observation. 53 Interestingly, 100 years after the works of Plateau and The difference DC between the maximum concentration at the front of the particle, and minimum C at the rear of the particle is shown in the inset. All simulation results are presented for dimensionless quantities defined in Sec. II C.

ARTICLE scitation.org/journal/phf
Marangoni the role of such impurities and Marangoni contribution in coalescence is still being resolved. 54 It should be pointed out that such effects will be less strong when the surfactant concentration is higher, and the interfacial equation of state is no longer linear, but a weaker dependency of the surface tension with increasing concentration is typically observed.
To investigate the effects of Ma, we study the interface at different Ma at constant Pe s ¼ 1. The interface compression and ð$ s Á u s Þ max are shown in supplementary material Fig. S12. Results show that the interface compressibility decreases significantly by adding a small amount of surfactant to the interface. If Ma > 50, its dependency on the concentration decreases, and there is a very slight change in the interface compression. As we have already pointed out, the Marangoni flow imposes a force on the particle and we expect the drag coefficient f to increase. Drag coefficients change from a diffusion dominant regime at Pe s ¼ 0:05 to a convection dominant regime at Pe s ¼ 5 [ Fig. 5(a)]. The Marangoni contribution F M to the force amplitude as a function of Ma at different Pe s is shown in Fig. 5(b). The F M increases with Ma, and it is higher at high Pe s . Moreover, Fig. 5(a) shows how the value of the drag coefficient converges to 11.5 for large Pe s and Ma. This value is consistent with that for an incompressible interface, as will be shown in Sec. IV B.
To explicitly show that surface dilation depends on both Pe s and Ma, we show the maximum interface dilation, normalized by the value for a clean, inviscid interface, as a function of both Pe s and Ma in Fig. 6. In this plot, the values on the axes are similar to the results for a clean, inviscid interface [supplementary material Sec. S1(ii)], since at small Ma the effect of Marangoni stresses becomes negligible, and at small Pe s diffusion dominates over the convection in the surface convection-diffusion (SCD) equation. Hence, in the latter case, the diffusion timescale for recovering a uniform distribution of the surfactant at the interface is much shorter than the convection timescale. In the limit of Pe s ! 0, the surfactant concentration remains uniform, and regaining the uniform distribution is instantaneous. For values of Ma and Pe s that are both non-zero, the contourlines in this plot can be described well by a function of the form Ma Pe s ¼ constant, clearly showing that is the product between Pe s and Ma that plays a crucial role in surface dilation.

B. Purely dilatational interfaces
For a purely expanding or compressing interface or for an interface with Bq 2 ) Bq 1 , the interface stress tensor (2) simplifies to Using the above relation in the tangential momentum jump at the interface, Eq. (22), we have revealing the dilatational properties of the surface to be the surface tension c and the surface dilatational viscosity Bq 2 . The expression in the parenthesis on the right hand side is also known as dynamic surface tension and can be expressed as and is thus composed of both thermodynamic and viscous parts. Identifying the reversible and irreversible contributions to the dynamic surface tension is a difficult task, although it can be achieved by realizing that the thermodynamic part will be influenced by mass transport considerations (and hence geometry) and the viscous part only depends on the local strain rate. For the limiting case of Pe s ! 0, the surfactant remains uniform, the interface tension is constant, and the related term in Eq. (29) does not play any role for the dynamics at the interface. If we assume that the left hand side is a finite constant value, when j s is getting larger, the interface compression $ s Á u s is getting smaller. Figure 7 shows our simulation results for the interface compression in the presence of the spherical particle that moves at constant speed in x-direction. The maximum dilation ð$ s Á u s Þ max is shown in the inset. We find that increasing Bq 2 helps to make the interface more incompressible, while for Bq 2 > 1000 the interface dilation does not change anymore significantly. The drag coefficient increases (by about one third) with increasing Bq 2 (Fig. 8) and is expected to reach a constant value at Bq 2 ! 1, when the interface is in an incompressible state, as shown in Fig. 8, f levels off at 11.5. As can be seen in Fig. 7, at this stage the interface is incompressible. These findings are consistent with the results for large Pe s and Ma, as shown in the previous section, but the

ARTICLE
scitation.org/journal/phf origin of the surface incompressibility is different in both cases: here, rheological surface stresses cause a resistance against surface compression, whereas in Sec. IV A surface incompressiblity arises due to Marangoni stresses, which are thermodynamic in nature.

C. Marangoni effects at dilatational interfaces
When Pe s % Oð1Þ or higher, both Marangoni flow and dilational viscosity can be responsible for the interface incompressibility. As we have shown in Eq. (29), it is not easy to distinguish these two phenomena. The relative importance of these two depends on Bq 2 and Ma. Figure 9 shows the surfactant concentration profiles and DC at different Bq 2 . At high Bq 2 , the Marangoni flow diminishes due to the suppression of the surfactant gradient. These results indicate that the Marangoni flow is important at low Bq 2 when the interface is still compressible. Maximum dilations ð$ s Á u s Þ max and DC at different Bq 2 ¼ 0; 1; 10; 100 as a function of surfactant concentration are shown in Fig. 10(a). At low surfactant concentration, the effect of Bq 2 in both quantities is significant. Calculating the magnitude of the Marangoni contribution to the force also shows that F M decreases with increasing Bq 2 [ Fig. 11

ARTICLE
scitation.org/journal/phf indicates that the Marangoni effects are still present. We can consider the interface with Bq 2 ¼ 100 as being "almost incompressible," and adding the surfactant makes the interface incompressible. Comparison [ Fig. 11(b)] of the drag coefficient at the interface with Bq 2 ¼ 0 and Bq 2 ¼ 1 shows negligible changes at different Ma.
The sensitivity of f on Ma diminishes at high Bq 2 , where the interface is incompressible, and Marangoni effects absent [ Fig. 11(b)].

D. Interplay between surface viscosities, Marangoni flow, and (in)compressibility of the interface
So far, we have assumed that the interface has a vanishing shear viscosity. However, the presence of an interfacial shear viscosity will also lead to an energy dissipation when shear deformations occur. The detailed nature of the interface flow field hence depends on the relative contributions of other parameters such as interface compressibility and Marangoni flow at the interface. Yet, in microrheological experiments often only the shear viscosity is calculated from the observed drag forces or diffusivities and assumptions about incompressibility are typically made. It is thus worthwhile to investigate the interplay between these effects and the interfacial shear viscosity.
A compressible interface exhibits the possibility for the occurrence of Marangoni flow at the limit of Bq 2 ¼ 0. For an interface with Pe s ! 0, the result shows that at Bq 1 < 10, interface dilation increases with increasing Bq 1 which means shear deformations cost more and more energy, the interface appears even more compressible, see Fig.  12. At Bq 1 > 10, a further increase in surface shear viscosity slightly decreases the interface compressibility, most likely as, due to even a small shear deformations at these higher viscosities, the dissipation of momentum now is more localized close to the particle. In general, the interface shear viscosity does not have strong effects on the interface compressibility, as the effects are mainly indirect.
For a monolayer with a finite P eclet number Pe s ¼ 0:5 (and Ma ¼ 10), for which now Marangoni effects are present, the interface dilation is also shown in Fig. 12. The behavior is similar to the case of Pe s ! 0, where the surfactant remains uniform. Results show that with increasing the interface shear viscosity, the interface compressibility increases.
An illustration of an increase in the interface compressibility with Bq 1 is suppression of the backflow in a confined flow. When a particle moves at the compressible interface, due to the confinement, it generates a backflow, where a part of the fluid in front of the particle moves in the opposite direction to compensate the motion of the particle. Figure 13 shows the interface velocity fields for low Bq 1 ¼ 0:1 and high Bq 1 ¼ 100 at Pe s ! 0 and Bq 2 ¼ 0. The simulation domain length is L ¼ 100. At low interface shear viscosity, there is a vortex at the interface. We can see that increasing the interface shear viscosity suppresses the vortex and causes the interface to move more uniformly. The center of the vortex, located at ð0; 0; Z b Þ, can be an indicator of the interface compressibility. With decreasing compressibility, the backflow increases and Z b is getting closer to the particle surface and vice versa.
Profiles of the x-component of the surface velocity along the z-direction, from the particle surface, at different Bq 1 are shown in

ARTICLE
scitation.org/journal/phf Fig. 14, with the corresponding Z b vs Bq 1 in the inset. These results demonstrate that upon increasing the interface shear viscosity, Z b moves further away from the particle, which indicates an increase in the interface compressibility with increasing z.
The uniform motion of the interface at high Bq 1 can also be observed in an unbounded flow when walls are at infinity. Figure 15 shows velocity field, interface dilatation, and shear component of surface stress tensor for an interface with finite compressbility,

ARTICLE
scitation.org/journal/phf Bq 2 ¼ 0:1, at three different interface shear viscosity Bq 1 2 f0:1; 1; 5g and limit of Pe s ! 0. The shear component of the surface velocity gradient is necessary to fully characterize the velocity gradient. The velocity fields show clearly the uniform motion at high Bq 1 . This uniform motion minimizes surface shear and increases the cost for dissipation of energy. This effect can be seen from the third column in Fig. 15. At low Bq 1 , the shear effect is very important at close distance from the surface of the particle, where there is a large gradient in shear component of the stress tensor. At high shear viscosity, the gradient decreases and shear effects happen at far distance from the particle. Suppression of the vortex and uniform motion also affects the Marangoni flow. Figure 16 shows the velocity field u s and the surfactant concentration at Ma ¼ 10 and finite Pe s ¼ 0:1. The interface has the same shear and dilatational viscosity as Fig. 15. To highlight the

ARTICLE
scitation.org/journal/phf qualitative role of the Marangoni flow, its contribution to the surface velocity field, Du s ¼ u s À u 0 s is also shown in Fig. 16, where u 0 s denotes the velocity field of the interface with uniform surfactant, at Pe s ! 0. At low Bq 1 , the Marangoni flow occurs near the particle surface. Suppression of the vortex at higher shear viscosity enforces the surfactant to find another path to flow. Therefore, the flow is less focused and occurs even at far distance from the particle (see the second column). The results also show that the maximum flow is getting further from the particle with increasing Bq 1 . This is likely due to moving the vortex by increasing the interface compressibility. This difficulty for the surfactant to flow creates a higher gradient in front and back of the particle (third column), which consequently increases the Marangoni force acting on the particle. The Marangoni contribution F M to the force and the interface viscous drag F s amplitudes are compared in Fig. 17. At low Bq 1 , F M increases significantly with Bq 1 , and it reaches a constant value with increasing Bq 1 . At high Bq 1 , the F M does not make a significant contribution to the total drag on the particle. An increase in F M with increasing surface shear viscosity was also reported by Elfring et al. 12 for the motion of a disk at the interface. This can again be attributed to the enhancement of dilatation/compressional deformations as shear deformations become more costly.
When considering shear and dilational viscosities, a wide range of values has been reported (Table II), but for simple surface active components, such as small molecular weight surfactants, the ratio of shear and dilatational surface viscosities has been suggested to be of comparable magnitude, H ¼ Bq 2 =Bq 1 ¼ Oð1Þ. 16,18 Depending on the chemical characteristics of the system at hand and the mobility of the surfactants (small molecular weight vs polymeric), the P eclet number might either be small or large compared to unity, we are going to discuss both types upon choosing Pe s 2 f0:1; 2:0g. All qualitative effects are captured by a finite Bq 1 ¼ 1 and moderate Ma ¼ 10, which is the choice we are making for the present purpose.
The surface flow field and its characteristics, as well as the surfactant concentration profiles, are shown for three different H 2 f0:1; 1; 10g and Pe s ¼ 0:1 in Fig. 18 (corresponding data for Pe s ¼ 2:0 are available in supplementary material Fig. S13). The different components of the velocity field in Fig. 18 again reveal the interplay between the shear and dilational/compression. This figure demonstrates how the relative contributions of compression/dilation and the shear deformations depend on the relative values of H and the magnitudes of the Marangoni and bulk to interfacial stresses, congruent with work on macroscopic rheometers, which employ mixed flows, such as the double wall sinusoidal ring 10 and the interfacial Cambridge extensional rheometer. 11,57 E. Influences of the surface viscosities on the drag coefficient of a spherical particle Although recent available techniques make the experimental realization of surface microreology relatively straightforward, especially in the particle tracking techniques, one has to rely on a hydrodynamic model of the monolayer to obtain interface quantities such as the interface shear viscosity. Particle tracking techniques are a branch of, so-called, microrheological techniques, which have advantages over macrorheology including smaller sample sizes, access to higher frequency ranges, the ability to measure small-scale material heterogeneities, and greater sensitivity. 9 In particle tracking, the position of a Brownian probe is recorded as a function of time, and the diffusivity D of the probe is calculated from a mean squared displacement (MSD) which in two dimensions is related to diffusivity by MSD ¼ hDr 2 ðsÞi ¼ 8Ds; (31) where r is the distance between a pair of particles, s is the lag time, and and hi denotes an ensemble average. In general, D is related to the particle's drag coefficient f by the Einstein-Sutherland relation 58,59 For finite Bq 2 the curves remain qualitatively unchanged, inline with the expectations from results presented in Fig. 11(a).  Fig. 19. For the larger Pe s , both the drag coefficient and the flow fields are basically unaffected by H over the experimentally relevant range. The surface flow field is only weakly compressed and dilated, while the surfactant concentration is distorted by about 10%. At the smaller Pe s the situation is reversed (Fig. 18). While the surfactants remain homogeneously spread at the interface, the surface velocity field exhibits larger shear gradients and weaker compression/dilation effects with increasing H.
According to these results, the effect of Pe s may be best detected from the C variation, while H affects the z-component of the surface velocity field, evaluated along the x ¼ z-diagonal. The above discussion pertains to small molecule surfactants; 12,60 for other types of structure interfaces, the Boussinesq-Scriven model is probably not adequate and strong viscoelastic effects are present, and dilatational elasticity becomes important. 9 However overall the trends and ideas observed in the simulations can be used to assess under which conditions the drag force is going to be determined by shear viscosities.
While Fig. 19 represents the drag coefficient for an interface in general cases when all dimentionless parameters are present, models

ARTICLE
scitation.org/journal/phf for the drag coefficient of protruding particle as a function of interface shear viscosity have been reported for the limiting cases of compressible and incompressible interfaces. Danov et al. 18 solved numerically the problem of a translating particle at a compressible Newtonian interface with constant surface tension. They reported the particle diffusion coefficient / 1=f at different contact angles as a function of Bq 1 and Bq 2 . We simulated the interface under similar conditions. The drag coefficient f and magnitude of bulk F b and interface viscous F s drags are shown in Fig. 20. Comparing our results with Danov's study for a particle with contact angle 90 shows that these two results are in a very good agreement [ Fig. 20(b)], given that the numerical methods differ and that we are operating in a differently bounded domain. Bonales et al. 61 have used Danov's theory to calculate the shear viscosity of two polymer Langmuir films. They compared the results with the data obtained by canal viscosimetry. Maestro et al. 4 used this theory at the glass transition and compared the results with particle tracking techniques. Sickert et al. 22 used Danov's theory in the study of monolayers of fatty acids and phospholipids. They have shown that the surface shear viscosity calculated from Danov's theory was lower than experimental data. A critical assessment of these findings can be found in Samaniuk and Vermant. 9 The incompressibility of interfaces is a widely accepted assumption and has been used in many pioneering works in the problem of particle motion at the interfaces and membranes in both nonprotruding particles 3,13,15 and protruding particle. 23,62 Fischer et al. 23 solved the problem of translation of a particle with radius R in an incompressible viscous monolayer with the surface shear viscosity g s , between two viscous phases with viscosities g 1 and g 2 . They solved the incompressible surface Stokes equations where F s is an external surface force parallel to the monolayer and Cj s represents the jump in C across interface. They present solutions to their model for drag coefficients on spherical particles at interfaces for the limiting cases where Bq 1 ( 1 or Bq 1 ) 1. They have obtained the following result for the translational drag coefficient f as a series expansion of Boussinesq number Bq ¼ g s =ðg 1 þ g 2 ÞR for 0 < Bq ( 1. For our setup Bq ¼ Bq 1 and thus Fitted expressions for f from numerical results give the following formulas for f 0 and f 1 : 23 where d is the signed distance from the apex of the sphere to the plane of the interface, which defines the contact angle. Equation (35) 23 For a translational drag on a half immersed sphere in a non viscous monolayer, they find f % f 0 % 11:7, which is higher than at a free surface (f ¼ 3p). This is due to surface incompressibility, which they attribute to Marangoni effects. The value of f 0 % 11:7 obtained by Fischer et al. (the above fit formula with unrevised a 0 ¼ 32 instead gives f 0 ¼ 11:08) for an incompressible and inviscid interface is slightly higher than the value of f % 11:5 we found for an incompressible interface due to Marangoni flow in Sec. IV A or dilational viscosity in Sec. IV B or combinations of both effects in Sec. IV C. These results are important and show that the drag coefficient of a particle at incompressible interfaces is independent of the origin of the incompressibility (dilatational viscosity, Marangoni effects or a combination of both) and that its higher value can not only be related to the Marangoni effects, as suggested by Fischer et al. The conditions required for surface incompressibility were not investigated in detail by Fischer et al. 23 Figure 11(a) shows that at Bq 2 ¼ 100 the Marangoni effects are still present but they do not have significant effects on the drag coefficient [ Fig. 11(b)]. These results do not support their claim that the Marangoni effects become very noticeable in the limit of vanishing surface compressibility. 23 We compared our results with Fischer's study for Bq 1 ( 1 for an incompressible interface with Bq 2 ¼ 10 6 using a linear fit to the data. At Bq 1 ! 0, the drag coefficient f ¼ f 0 % 11:52 is smaller than the value reported by Fischer et al., 23 and for the slope f 1 we obtain f 1 % 6:85. Drag coefficients of the particle for intermediate Bq 1 (1 Bq 1 10) and high values (10 < Bq 1 10 4 ) at an incompressible interface (with Bq 2 ¼ 10 6 ) are shown in Fig. 21 regime. A comparison between the drag coefficient at small and intermediate Bq 1 clearly shows a noticeable difference between these possibly experimentally relevant two regimes (Table II). A similar analysis for a compressible interface (Bq 2 ¼ 0:1; Pe s ! 0) is captured by the analytical expressions f ¼ 11:234 þ 1:931Bq 1 and ln f ¼ 1:07 þ 0:857 ln ðBq 1 Þ for the intermediate and high Bq 1 , respectively, in the same range of Bq 1 used for the incompressible fluid in Fig. 21.
Fischer's model has been used for the regime of Bq 1 ( 1 to calculate the interface shear viscosity from D, with R, g 1 and T at hand. 4,9,22,30,63,64 However, if one utilizes it within the 0:1 < Bq 1 < 10 regime, there is at least an order-of-magnitude agreement with estimates of g s calculated from an extrapolation of the two limiting solutions of the Fischer model. 9 The model presented in this work for drag coefficient at 1 Bq 1 10-in dimensional form f ¼ g 1 Rð14:830 þ 3:677 Bq 1 Þcan be used for such an intermediate regime.

V. CONCLUSIONS
This work offers a comprehensive study of the effects of interfacial properties on the drag of a spherical probe partially immersed into a viscous bulk fluid. We calculated the force on a spherical probe particle translating at constant velocity and centered at the flat interface, while the interface was initially covered uniformly by a surface active component. The translating particle compresses the interface at its front and dilates at its rear. We demonstrated how the quantitative occurrence of this phenomenon is affected by the properties of the surface active component, namely the concentration, the interfacial equation of state and the P eclet number, and the interface shear and dilatational viscosities. The presence of a surface active material was shown to reduce the interface compressibility and this effect strongly depends on the surface P eclet number. We quantified how a nonvanishing interface dilatational viscosity helps to reduce the local interface compressibility, while shear viscosity has only a smaller, opposite, and obviously indirect effect. A low value for the shear viscosity is found to provide an easy way to deform locally and reduces the divergence of the surface velocity field.
For the special case of inviscid interfaces carrying surfactants, the product between Ma and Pe s plays a crucial role in determining the dilation of the interface, and we show that for MaPe s > 15, the interface dilation is only 20% as compared to a clean, inviscid interface. Increasing the product between Ma Pe s further can thus make the interface effectively incompressible.
When interfacial viscosities are included as well, the Marangoni contribution F M to the forces experienced by the particle is found to be affected by the magnitude of both interfacial shear and dilatational viscosities. The particle induces a flow field that compresses/dilates the interface, which leads to surfactant concentration variations, and thus to variations in the interfacial tension along the interface. The flow field depends on the rheological properties of the interface, here to be taken purely viscous, and the interface P eclet number. A non-uniform distribution of surfactants leads to Marangoni flows, but these flows are coupled to the surface viscosities of the interface, resulting in a complex interplay of effects. Markedly, our simulation results showed that interface viscosities have an opposite effect on the Marangoni force F M acting on the particle; dilatational viscosity decreases but shear viscosity slightly increases F M .
The present work demonstrated that an interface can be incompressible by dilatational viscosity, Marangoni effects, or combinations of both. The drag coefficient of the particle at incompressible interfaces for these three cases, within the limit of vanishing interface shear viscosity, exhibited the same value of 11.5. We also calculated the particle's drag coefficient over a wide range of dimensionless interface shear viscosities (Bq 1 ) for compressible and incompressible interfaces and provided approximate analytical expressions for them.
The results presented in this work validate conclusions from the literature that were limited to a narrow asymptotic regime and only valid for non-protruding particles in shallow subphases. 12 The results of the present study should therefore help interpreting shear and dilatational rheological, as well as particle tracking experiments, where particles are generally protruding and the subphase should be deep in order to effectively probe the interface. Moreover, our numerical method allowed us to explore the parameter space far from the asymptotic regimes.
Here, we focused on a spherical probe symmetrically immersed at the interface, with a contact angle equal to 90 . An extension to non-spherical probes to study the effects of the particle aspect ratio, orientation, and contact angle is straightforward and subject of a subsequent study. We here assumed that the surface tension decreases linearly with increasing surfactant concentration. This assumption is limited to dilute condutions 21,65 and should be released within the otherwise unchanged framework if one aims for a quantitative comparison with experimental data, provided the nonlinear features of the surface tension have been fully resolved beforehand, which would extend the linearized equation used here. Nonequilibrium molecular dynamics and Monte Carlo simulations along the lines indicated in recent works 56,66-68 for complex interfaces may help to provide a route to calculate the missing ingredients, the material functions c, and surface viscosities as a function of Ma for atomistically detailed or coarsegrained models of amphiphilic surfactants.

SUPPLEMENTARY MATERIAL
See the supplementary material for the validation of the FEM program, mesh convergence, effects of subphase depths, transient behavior, the coefficient a 0 , and additional data, table, and figures referred to within the manuscript.

DEDICATION
This manuscript is dedicated to the memory of Dr. Robert "Bob" Byron Bird.