1. Introduction
The fluid-structure interaction between flows and boundaries is a central situation in continuum mechanics, encountered at many length and velocity scales. A classical example is lubrication, where the addition of a liquid film – a lubricant – between two contacting objects allows for a drastic reduction of the friction between them. Such a process occurs in a large variety of contexts with hard materials, such as roller bearings, pistons and gears in industry (Dowson & Higginson Reference Dowson and Higginson2014), or faults (Brodsky & Kanamori Reference Brodsky and Kanamori2001) and landslides (Campbell Reference Campbell1989) in geological settings. At large velocity, or moderate loading, the liquid film is continuous with no direct contact between the solids. When the solids are deformable, the friction force can be described using elastohydrodynamic (EHD) models within the soft-lubrication approximation (Dowson & Higginson Reference Dowson and Higginson2014).
The previous EHD coupling is also widely encountered in soft condensed matter, but at very different pressure and velocity scales (Karan, Chakraborty & Chakraborty Reference Karan, Chakraborty and Chakraborty2018). Examples encompass the remarkable frictional properties of eyelids (Jones et al. Reference Jones, Fulford, Please, McElwain and Collins2008) and cartilaginous joints (Mow, Holmes & Lai Reference Mow, Holmes and Lai1984; Jahn, Seror & Klein Reference Jahn, Seror and Klein2016), as well as biomimetic gels (Gong Reference Gong2006) and rubbers (Sekimoto & Leibler Reference Sekimoto and Leibler1993; Moyle et al. Reference Moyle, Wu, Khripin, Bremond, Hui and Jagota2020; Wu et al. Reference Wu, Moyle, Jagota and Hui2020; Hui et al. Reference Hui, Wu, Jagota and Khripin2021). Of interest as well are the collisions and rebounds of spheres in viscous environments (Davis, Serayssol & Hinch Reference Davis, Serayssol and Hinch1986; Gondret, Lance & Petit Reference Gondret, Lance and Petit2002; Tan, Wang & Frechette Reference Tan, Wang and Frechette2019), the rheological properties of soft suspensions and pastes (Sekimoto & Leibler Reference Sekimoto and Leibler1993; Meeker, Bonnecaze & Cloitre Reference Meeker, Bonnecaze and Cloitre2004), and the self-similar properties of the contact (Snoeijer, Eggers & Venner Reference Snoeijer, Eggers and Venner2013).
In the last decade, EHD interactions have been of great interest in the materials science community with the emergence of contactless rheological methods to measure the mechanical properties of confined liquids and soft surfaces (Chan, Klaseboer & Manica Reference Chan, Klaseboer and Manica2009; Vakarelski et al. Reference Vakarelski, Manica, Tang, O'Shea, Stevens, Grieser, Dagastine and Chan2010; Leroy & Charlaix Reference Leroy and Charlaix2011; Leroy et al. Reference Leroy, Steinberger, Cottin-Bizonne, Restagno, Léger and Charlaix2012; Villey et al. Reference Villey, Martinot, Cottin-Bizonne, Phaner-Goutorbe, Léger, Restagno and Charlaix2013; Wang, Dhong & Frechette Reference Wang, Dhong and Frechette2015; Guan et al. Reference Guan, Charlaix, Qi and Tong2017; Wang et al. Reference Wang, Pilkington, Dhong and Frechette2017a; Wang, Tan & Frechette Reference Wang, Tan and Frechette2017b; Wang et al. Reference Wang, Zeng, Alem, Zhang, Charlaix and Maali2018; Lainé et al. Reference Lainé, Jubin, Canale, Bocquet, Siria, Donaldson and Niguès2019; Bertin et al. Reference Bertin, Zhang, Boisgard, Grauby-Heywang, Raphael, Salez and Maali2021). Typically, in such experimental systems, a spherical colloidal probe is immersed in a fluid and driven to oscillate, with a nanometric amplitude, near a surface of interest. The force exerted on the probe is measured by an atomic force microscope, a surface force apparatus or a tuning-fork microscope, and depends on the properties of both the fluid and the solid boundary.
Generally, an object that moves in a confined fluid environment experiences an enhanced drag force with respect to the bulk Stokes law, as a result of the boundary-induced flow modification (Happel & Brenner Reference Happel and Brenner2012). Furthermore, near a soft wall, the hydrodynamic interactions are modified by the deformation of the boundary that they generate, yielding a nonlinear coupling. Perturbation methods, assuming a small deformation of the interface, have been employed in order to calculate the soft-lubrication interactions exerted on a free infinite cylinder immersed in a viscous fluid and near a thin compressible elastic material (Salez & Mahadevan Reference Salez and Mahadevan2015). In particular, interesting inertial-like features have been predicted despite the low-
$Re$-number aspect of the flow.
Perhaps the most emblematic example of soft-lubrication interaction is the non-inertial lift force predicted for a particle sliding near a soft boundary (Sekimoto & Leibler Reference Sekimoto and Leibler1993; Beaucourt, Biben & Misbah Reference Beaucourt, Biben and Misbah2004; Skotheim & Mahadevan Reference Skotheim and Mahadevan2004, Reference Skotheim and Mahadevan2005; Urzay, Llewellyn Smith & Glover Reference Urzay, Llewellyn Smith and Glover2007; Urzay Reference Urzay2010). It might have important implications for advected biological entities, such as red blood cells (Grandchamp et al. Reference Grandchamp, Coupier, Srivastav, Minetti and Podgorski2013) and vesicles (Abkarian, Lartigue & Viallat Reference Abkarian, Lartigue and Viallat2002). Only recently, the associated dynamical repulsion from an immersed soft interface has been studied experimentally. A preliminary qualitative observation was reported in the context of smart lubricants and elastic polyelectrolytes (Bouchet et al. Reference Bouchet, Cazeneuve, Baghdadli, Luengo and Drummond2015). Then a study involving the sliding of an immersed macroscopic cylinder along an inclined plane pre-coated with a thin layer of gel, showed quantitatively an effective reduction of friction induced by the EHD lift force (Saintyves et al. Reference Saintyves, Jules, Salez and Mahadevan2016). Subsequently, the same effect was observed in the trajectories of micrometric spherical beads within a microfluidic channel coated with a biomimetic polymer layer (Davies et al. Reference Davies, Débarre, El Amri, Verdier, Richter and Bureau2018), and through the sedimentation of a macroscopic sphere near a pre-tensed suspended elastic membrane (Rallabandi et al. Reference Rallabandi, Oppenheimer, Zion and Stone2018). Finally, direct measurements of the EHD lift force for two types of elastic materials have been performed at small scales, using surface force apparatus and atomic force microscopy, respectively (Vialar et al. Reference Vialar, Merzeau, Giasson and Drummond2019; Zhang et al. Reference Zhang, Bertin, Arshad, Raphael, Salez and Maali2020).
Despite the increasing number of EHD studies involving spherical probes, the soft-lubrication interactions of a free spherical object immersed in a viscous fluid and moving near an elastic substrate still have to be calculated. In the present article, we aim at filling this gap by deriving a soft-lubrication perturbation theory, in order to compute all the forces and torque for this problem, at first order in dimensionless compliance.
The article is organized as follows. First, we introduce the soft-lubrication framework for a sphere translating near a soft planar surface, in both normal and tangential directions. The substrate deformation is assumed to follow the constitutive response of a linear elastic semi-infinite material. Then we follow a perturbative approach, assuming the substrate deformation to be small with respect to the fluid-gap thickness, which allows us to find the normal and tangential forces as well as the torque experienced by the sphere, at first order in dimensionless compliance. Finally, we discuss the rotation of the sphere, before providing concluding remarks. Besides, in Appendices A–D, the EHD forces are computed analytically using the Lorentz reciprocal theorem, while the procedure introduced in the main text is reproduced for the compressible and incompressible responses of a thin material.
2. Model
The system is depicted in figure 1. We consider a sphere of radius 
$a$, immersed in a Newtonian fluid of dynamic shear viscosity 
$\eta$ and density 
$\rho$. The sphere is moving with a tangential velocity 
$\boldsymbol {u}(t) = u(t)\,\boldsymbol {e}_x$ directed along the 
$x$-axis (by definition of the latter axis), where 
$\boldsymbol {e}_j$ denotes the unit vector along 
$j$. In this first part, we assume that the sphere does not rotate, i.e. the angular velocity reads 
$\boldsymbol {\varOmega }=\boldsymbol {0}$. The sphere is placed at a time-dependent distance 
$d(t)$ (thus a 
$\dot {d} \boldsymbol {e}_z$ normal velocity of the sphere) of an isotropic and homogeneous linear elastic substrate of Lamé coefficients 
$\lambda$ and 
$\mu$, with a reference undeformed flat surface in the 
$xy$-plane at 
$z=0$. We suppose that the sphere–wall distance is small with respect to the sphere radius, such that the lubrication approximation is valid. The fluid inertia is neglected here. Specifically, we assume 
$Re (d/a) \ll 1$, with the Reynolds number 
$Re = \rho u a /\eta$. Furthermore, we suppose that the typical time scale of variation of the sphere velocity is much larger than the diffusion time scale of vorticity that scales as 
$d^2/(\eta /\rho )$, such that the flow is described by the steady Stokes equations. This amounts to assuming that 
$\lvert \dot {u}/u \rvert \ll \eta /(\rho d^2)$ and 
$\lvert \ddot {d}/\dot {d} \rvert \ll \eta /(\rho d^2)$. We stress that slippage at solid surfaces modifies the lubrication pressure as well as the EHD interaction (Vinogradova & Feuillebois Reference Vinogradova and Feuillebois2000), which is of importance for flows at the nanoscale (Bocquet & Charlaix Reference Bocquet and Charlaix2010). Here, we ignore this effect and no-slip boundary conditions are assumed at both the sphere and wall surfaces. Finally, the system is equivalent to a sphere at rest near a wall translating with a 
$-(d(t)\,\boldsymbol {e}_z+\boldsymbol {u}(t))$ velocity. In such a framework, the fluid velocity field can be written as
\begin{equation} \boldsymbol{v}(\boldsymbol{r},z,t) = \frac{\boldsymbol{\nabla} p(\boldsymbol{r},t)}{2\eta} (z - h_0(r,t)) (z-\delta(\boldsymbol{r},t)) - \boldsymbol{u}(t)\,\frac{h_0(r,t)-z}{h_0(r,t)-\delta(\boldsymbol{r},t)}, \end{equation}
where 
$\boldsymbol {r}=(r,\theta )$ is the position in the tangential plane 
$xy$, 
$\boldsymbol {\nabla }$ is the two-dimensional gradient operator on 
$xy$, 
$\delta (\boldsymbol {r},t)$ is the substrate deformation, and 
$z = h_0(r,t)$ is the sphere surface. Near contact, the latter can be approximated by its parabolic expansion 
$h_0(r,t) \simeq d(t) + r^2/(2a)$. Volume conservation further leads to the Reynolds equation
\begin{equation} \partial_t h(\boldsymbol{r}, t) = \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\frac{h^3(\boldsymbol{r}, t)}{12\eta}\,\boldsymbol{\nabla} p(\boldsymbol{r}, t) + \frac{h(\boldsymbol{r}, t)}{2}\,\boldsymbol{u}(t) \right), \end{equation}
where 
$h(\boldsymbol {r}, t) = h_0(r,t) - \delta (\boldsymbol {r}, t)$ is the fluid-gap thickness. In this first part, we assume that the constitutive elastic response is linear and instantaneous, and that the substrate is a semi-infinite medium, such that the deformation reads (Davis et al. Reference Davis, Serayssol and Hinch1986)
\begin{equation} \delta(\boldsymbol{r}, t) ={-}\frac{(\lambda+2\mu)}{4{\rm \pi} \mu(\lambda+\mu)} \int_{\mathbb{R}^2}\textrm{d}^2\boldsymbol{x}\,\frac{p(\boldsymbol{x},t)}{\lvert\boldsymbol{r} - \boldsymbol{x} \rvert}. \end{equation}We non-dimensionalize the problem through
$$\begin{gather} h(\boldsymbol{r}, t) = d^* H(\boldsymbol{R}, T), \quad \boldsymbol{r} = \sqrt{2ad^*} \boldsymbol{R}, \quad d(t) = d^* D(T), \quad \delta(\boldsymbol{r}, t) = d^* {\rm \Delta}(\boldsymbol{R}, T), \end{gather}$$
$$\begin{gather}p(\boldsymbol{r}, t) = \frac{\eta u^* \sqrt{2ad^*}}{d^{*2}}\, P(\boldsymbol{R}, T), \quad \boldsymbol{u}(t) = u^* U(T)\,\boldsymbol{e}_x, \quad \boldsymbol{v} = u^* \boldsymbol{V},\quad t = \frac{\sqrt{2ad^*}}{u^*} T, \end{gather}$$
where 
$d^*$ and 
$u^*$ are characteristic fluid-gap distance and tangential velocity, respectively. The governing equations are then
$$\begin{gather} 12\partial_T H(\boldsymbol{R}, T) = \boldsymbol{\nabla} \boldsymbol{\cdot} \left(H^3(\boldsymbol{R}, T)\,\boldsymbol{\nabla} P(\boldsymbol{R}, T) + 6 H(\boldsymbol{R}, T)\,\boldsymbol{U}(T) \right), \end{gather}$$
$$\begin{gather}H(\boldsymbol{R}, T) = D(T) + R^2 - {\rm \Delta}(\boldsymbol{R}, T) \end{gather}$$and
\begin{equation} {\rm \Delta}(\boldsymbol{R}, T) ={-}\kappa \int_{\mathbb{R}^2} \textrm{d}^2\boldsymbol{X}\, \frac{P(\boldsymbol{X},T)}{4{\rm \pi} \lvert\boldsymbol{R} - \boldsymbol{X} \rvert}, \end{equation}where we introduced the dimensionless compliance
\begin{equation} \kappa=\frac{2\eta u^* a(\lambda+2\mu)}{d^{*2}\mu(\lambda+\mu)}. \end{equation}
The latter is the only dimensionless parameter in the problem. When 
$\kappa$ is small with respect to unity, it corresponds to the ratio between two length scales: the typical substrate deformation 
$\delta \sim {2\eta u^* a(\lambda +2\mu )}/{d^*\mu (\lambda +\mu )}$ induced by a tangential velocity 
$u^*$, and the typical fluid-gap thickness 
$d^*$. Throughout the article, we focus on the small-deformation regime of soft-lubrication where 
$\kappa \ll 1$ (Essink et al. Reference Essink, Pandey, Karpitschka, Venner and Snoeijer2021).
Figure 1. Schematic of the system. A rigid sphere of surface 
$\mathcal {S}_0$ is freely moving in a viscous fluid, near a soft wall of surface 
$\mathcal {S}_{{w}}$ in the flat undeformed state. The lubrication pressure field deforms the latter, which induces an EHD coupling, with forces and torque exerted on the sphere. Note that the surface deformation is magnified for clarity, but that we restrict the analysis to the 
$\delta \ll d$ case.
3. Perturbation theory
We perform a perturbation analysis at small 
$\kappa$ (Sekimoto & Leibler Reference Sekimoto and Leibler1993; Beaucourt et al. Reference Beaucourt, Biben and Misbah2004; Skotheim & Mahadevan Reference Skotheim and Mahadevan2004, Reference Skotheim and Mahadevan2005; Urzay et al. Reference Urzay, Llewellyn Smith and Glover2007; Urzay Reference Urzay2010; Salez & Mahadevan Reference Salez and Mahadevan2015; Pandey et al. Reference Pandey, Karpitschka, Venner and Snoeijer2016; Rallabandi et al. Reference Rallabandi, Saintyves, Jules, Salez, Schönecker, Mahadevan and Stone2017; Saintyves et al. Reference Saintyves, Rallabandi, Jules, Ault, Salez, Schönecker, Stone and Mahadevan2020; Zhang et al. Reference Zhang, Bertin, Arshad, Raphael, Salez and Maali2020), as follows:
$$\begin{gather} H(\boldsymbol{R}, T) = H_0(\boldsymbol{R},T) + \kappa\, H_1(\boldsymbol{R}, T) + O(\kappa^2), \end{gather}$$
$$\begin{gather}P(\boldsymbol{R}, T) = P_0(\boldsymbol{R},T) + \kappa\,P_1(\boldsymbol{R}, T) + O(\kappa^2), \end{gather}$$
where the subscript 
$0$ corresponds to the solution for a rigid wall, with 
$H_0(\boldsymbol {R},T) = D(T) + R^2$.
3.1. Zeroth-order solution: rigid wall
At zeroth order 
$O(\kappa ^0)$, (2.6) reads
\begin{equation} 12 \dot{D} =\boldsymbol{\nabla} \boldsymbol{\cdot} \left(H_0^3\,\boldsymbol{\nabla} P_0 + 6 H_0\boldsymbol{U} \right). \end{equation}In polar coordinates, (3.3) can be rewritten as
\begin{align} \mathcal{L}\cdot P_0 = R^2 \partial_R^2 P_0 + \left(R + \frac{6R^3}{D+R^2} \right) \partial_R P_0 + \partial_\theta^2 P_0 = \frac{R^2}{(D + R^2)^3} \left(12 \dot{D} - 12 R \cos\theta U\right), \end{align}
where 
$\mathcal {L}$ is a linear operator. We solve this equation using an angular-mode decomposition:
\begin{equation} P_0(\boldsymbol{R},T) = P_0^{(0)}(R,T) + P_0^{(1)}(R,T) \cos\theta, \end{equation}where the two coefficients are solutions of the ordinary differential equations
\begin{equation} R^2\,\frac{\textrm{d}^2 P_0^{(0)}}{\textrm{d}R^2} + \left( R + \frac{6R^3}{D+R^2} \right) \frac{\textrm{d} P_0^{(0)}}{\textrm{d}R} = 12\, \frac{R^2 \dot{D}}{(D + R^2)^3}, \end{equation}
\begin{equation} R^2\,\frac{\textrm{d}^2 P_0^{(1)}}{\textrm{d}R^2} + \left( R + \frac{6R^3}{D+R^2} \right) \frac{\textrm{d} P_0^{(1)}}{\textrm{d}R} - P_0^{(1)} ={-}12\,\frac{R^3 U}{(D + R^2)^3}. \end{equation}
In accordance with the boundary conditions 
$P(R\rightarrow \infty ) = 0$ and 
$P(R = 0) < \infty$, the solution is thus
\begin{equation} P_0(\boldsymbol{R},T) ={-}\frac{3 \dot{D}}{2 (D+R^2)^2} + \frac{6 R U \cos\theta}{5(D+ R^2)^2}. \end{equation}
The first-order substrate deformation 
$H_1$ can then be computed from (2.8) at order 
$O(\kappa )$:
\begin{equation} H_1(\boldsymbol{R}, T) = \int_{\mathbb{R}^2} \textrm{d}^2\boldsymbol{X}\, \frac{P_0(\boldsymbol{X},T)}{4{\rm \pi} \lvert\boldsymbol{R} - \boldsymbol{X} \rvert}. \end{equation}
Using e.g. the spatial Fourier transform 
$\tilde {H}_1(\boldsymbol {K}) = \int _{\mathbb {R}^2} H_1(\boldsymbol {R})\exp ({-\textrm {i}\boldsymbol {R}\boldsymbol {\cdot }\boldsymbol {K}})\,\textrm {d}^2\boldsymbol {R}$, we find
\begin{align} H_1(\boldsymbol{R}, T) ={-}\frac{3\dot{D}}{8\sqrt{D}}\, \frac{\mathcal{E}({-}R^2/D)}{D+ R^2} + \frac{3U}{10R\sqrt{D}} \left(- \frac{D\,\mathcal{E}({-}R^2/D)}{D+R^2} + \mathcal{K}({-}R^2/D) \right) \cos \theta, \end{align}
where 
$\mathcal {K}$ and 
$\mathcal {E}$ are the complete elliptic integrals of the first and second kinds (Abramowitz & Stegun Reference Abramowitz and Stegun1964). The dimensionless substrate deformations are plotted in figure 2. In figure 2(a), the sphere is moving tangentially to the substrate with a unit velocity 
$U = 1$. The deformation exhibits a dipolar symmetry, with a negative sign (i.e. the substrate is compressed) at the front. Besides, the isotropic term generated by a sphere moving normally to the substrate is shown in figure 2(b). In particular, for a sphere approaching the substrate, the latter is compressed.
Figure 2. Dimensionless deformation fields at the free surface of the soft substrate, for a sphere placed at a unit distance 
$D = 1$ and for two modes of motion: (a) the sphere velocity is directed tangentially to the substrate, along the 
$x$-axis (see black arrow), and is fixed to a unit value 
$U = 1$; (b) the sphere is approaching the substrate normally (see black cross) with a unit velocity 
$\dot {D} = -1$.
3.2. First-order solution
We can now compute the first-order pressure field 
$P_1$ from (2.6) at order 
$O(\kappa )$:
\begin{equation} 12\partial_T H_1 = \boldsymbol{\nabla} \boldsymbol{\cdot}\left(H_0^3\,\boldsymbol{\nabla} P_1 + 3H_0^2 H_1\, \boldsymbol{\nabla} P_0 + 6 H_1\boldsymbol{U} \right). \end{equation}
Invoking the same linear operator 
$\mathcal {L}$ as in (3.4), we can rewrite (3.11) as
\begin{equation} \mathcal{L} \cdot P_1 = \frac{R^2}{H_0^3} \left(12 \partial_T H_1 - \boldsymbol{\nabla} \boldsymbol{\cdot} \left[3H_0^2 H_1 \boldsymbol{\nabla}\,P_0 + 6 H_1\boldsymbol{U}\right]\right). \end{equation}We then expand all the terms in the right-hand side of (3.12), and we perform once again the angular-mode decomposition:
\begin{equation} \mathcal{L} \cdot P_1 = F_0(R,T) + F_1(R,T) \cos\theta + F_2(R,T) \cos2\theta, \end{equation}where we have introduced the auxiliary functions
\begin{align} F_0(R,T) ={}& \frac{18R^2 U^2}{25 D^{1/2}(D+R^2)^6} \left[({-}10D^2 +2DR^2)\,\mathcal{E}\left(-\frac{R^2}{D}\right) \right. \nonumber\\ &\left.{}+ (8D^2+7DR^2-R^4)\,\mathcal{K}\left(-\frac{R^2}{D}\right)\right] \nonumber\\ &{}+ \frac{9R^2 \dot{D}^2}{4 D^{3/2}(D+R^2)^6} \left[(13D^2 + 3R^2D + 2R^4)\,\mathcal{E}\left(-\frac{R^2}{D}\right) \right. \nonumber\\ &\left.{}+ ({-}4D^2 - 5R^2D-R^4)\,\mathcal{K}\left(-\frac{R^2}{D}\right)\right] \nonumber\\ &{}-\frac{9R^2 \ddot{D}\,\mathcal{E}\left(-\dfrac{R^2}{D}\right)}{2 D^{1/2}(D+R^2)^4} \end{align}and
\begin{align} F_1(R,T) ={}& -\frac{27 R U\dot{D} }{5 D^{1/2}(D+R^2)^6} \left[ ({-}2D^2 +7DR^2+R^4)\,\mathcal{E}\left(-\frac{R^2}{D}\right) \right. \nonumber\\ &\left.{}+ 2(D+R^2)(D-R^2)\,\mathcal{K}\left(-\frac{R^2}{D}\right)\right] \nonumber\\ &{}-\frac{18R \dot{U}}{5 D^{1/2}(D+R^2)^4} \left[- D\,\mathcal{E}\left(-\frac{R^2}{D}\right) + (D + R^2)\, \mathcal{K}\left(-\frac{R^2}{D}\right)\right]. \end{align}
Note that we have not provided 
$F_2$ as it does not contribute to the forces and torque. Also note that by setting 
$D(T) = 1$ in the latter expressions, we self-consistently recover the expression of Zhang et al. (Reference Zhang, Bertin, Arshad, Raphael, Salez and Maali2020). Invoking the angular-mode decomposition 
$P_1(\boldsymbol {R},T) = P_1^{(0)}(R,T) + P_1^{(1)}(R,T)\cos \theta + P_1^{(2)}(R,T) \cos 2\theta$, we get, in particular,
$$\begin{gather} R^2\,\frac{\textrm{d}^2 P_1^{(0)}}{\textrm{d}R^2} + \left( R + \frac{6R^3}{D+R^2} \right) \frac{\textrm{d} P_1^{(0)}}{\textrm{d}R} = F_0(R,T), \end{gather}$$
$$\begin{gather}R^2\,\frac{\textrm{d}^2 P_1^{(1)}}{\textrm{d}R^2} + \left( R + \frac{6R^3}{D+R^2} \right) \frac{\textrm{d} P_1^{(1)}}{\textrm{d}R} - P_1^{(1)} = F_1(R,T). \end{gather}$$
Using scaling arguments, we can write the two relevant first-order pressure components 
$P_1^{(i)}$ as
\begin{equation} P_1^{(0)} = \frac{U^2}{D^{7/2}}\,\phi_{U^2}(R/\sqrt{D}) + \frac{\dot{D}^2}{D^{9/2}}\,\phi_{\dot{D}^2}(R/\sqrt{D}) + \frac{\ddot{D}}{D^{7/2}}\,\phi_{\ddot{D}}(R/\sqrt{D}) \end{equation}and
\begin{equation} P_1^{(1)} = \frac{U\dot{D} }{D^4}\,\phi_{U \dot{D}}(R/\sqrt{D}) + \frac{\dot{U}}{D^3}\,\phi_{\dot{U}}(R/\sqrt{D}), \end{equation}
where the 
$\phi _i$ are five dimensionless scaling functions that depend on the self-similar variable 
$R/\sqrt {D}$ only. Equations (3.18) and (3.19) can be solved numerically with a Runge–Kutta algorithm and a shooting parameter in order to ensure the boundary condition 
$P_1(R \rightarrow \infty, \theta, T) = 0$. All the scaling functions are plotted in figures 3 and 4.
As a remark, we recall that the substrate deformation is induced by the flow, and that at first order it is linear in the velocity field (see (3.10)). Moreover, the volume-conservation equation involves the time derivative of the fluid-layer thickness, and thus in particular the time derivative of the substrate deformation. As a consequence, when calculating the first-order EHD pressure field, we find terms (and thus forces and torques) that are proportional to the acceleration components 
$\ddot {D}$ and 
$\dot {U}$ of the sphere. At first sight, these original inertial-like features might seem inconsistent with steady Stokes flows, but are in fact independent of the fluid density and solely induced by the intimate EHD coupling.
3.3. Forces and torque
The force 
$\boldsymbol {F}$ exerted by the fluid on the sphere is given by
\begin{equation} \boldsymbol{F}= \int_{\mathcal{S}_0} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma} \, \mathrm{d}s, \end{equation}
where 
$\boldsymbol {\sigma } = - p \boldsymbol {I} + \eta (\boldsymbol {\nabla } \boldsymbol {v} + \boldsymbol {\nabla } \boldsymbol {v}^\text {T})$ is the fluid stress tensor, 
$\boldsymbol {n}$ is the unit vector normal to the sphere surface and pointing towards the fluid, and 
$\boldsymbol {I}$ is the identity tensor. Within the lubrication approximation, the fluid stress tensor reads 
$\boldsymbol {\sigma } \simeq - p \boldsymbol {I} + \eta \boldsymbol {e}_z \partial _z \boldsymbol {v}$. One can then evaluate the normal force as
\begin{align} F_z = \int_{\mathbb{R}^2} p(\boldsymbol{r})\,\textrm{d}^2 \boldsymbol{r} ={}& -\frac{6{\rm \pi} \eta a^2 \dot{d}}{d} + 0.41623\,\frac{\eta^2 u^2(\lambda+2\mu)}{\mu(\lambda+\mu)} \left(\frac{a}{d}\right)^{5/2} \nonumber\\ &{}-41.912\,\frac{\eta^2 \dot{d}^2(\lambda+2\mu)}{\mu(\lambda+\mu)} \left(\frac{a}{d}\right)^{7/2} + 18.499\,\frac{\eta^2 \ddot{d} a(\lambda+2\mu)}{\mu(\lambda+\mu)} \left(\frac{a}{d}\right)^{5/2}, \end{align}
where the prefactors have been estimated numerically using (3.18). We recover in particular the classical Reynolds force 
$-6{\rm \pi} \eta a^2 \dot {d}/d$ at zeroth order, i.e. near a rigid wall. We stress that tangential motions do not induce any normal force at zeroth order in 
$\kappa$, as the corresponding pressure field is antisymmetric in 
$x$ (see (3.4)). In contrast, such motions do induce a lift force at first order in 
$\kappa$, due to the symmetry breaking of the contact geometry associated with the elastic deformation. Interestingly as well, normal motions generate a viscous adhesive force at first order in compliance (Wang, Feng & Frechette Reference Wang, Feng and Frechette2020). Besides, an original EHD force proportional to the sphere's normal acceleration is also found, as discussed previously. Finally, in Appendix A, and following previous works (Rallabandi et al. Reference Rallabandi, Saintyves, Jules, Salez, Schönecker, Mahadevan and Stone2017; Daddi-Moussa-Ider et al. Reference Daddi-Moussa-Ider, Rallabandi, Gekle and Stone2018; Rallabandi et al. Reference Rallabandi, Oppenheimer, Zion and Stone2018; Masoud & Stone Reference Masoud and Stone2019), we use the Lorentz reciprocal theorem in order to recover the prefactors of (3.21) analytically, which gives, respectively, 
${243 {\rm \pi}^3}/{12800 \sqrt {2}} \approx 0.41623$, 
${3915 {\rm \pi}^3}/{2048 \sqrt {2}} \approx 41.912$ and 
${27{\rm \pi} ^3}/{32\sqrt {2}} \approx 18.499$. We note that the latter is in agreement with the result of the linear-response theory derived in Leroy & Charlaix (Reference Leroy and Charlaix2011). Furthermore, we recover the lift prefactor (0.416) obtained previously numerically (Zhang et al. Reference Zhang, Bertin, Arshad, Raphael, Salez and Maali2020), as well as analytically in a recently published work (Kargar-Estahbanati & Rallabandi Reference Kargar-Estahbanati and Rallabandi2021).
Similarly, the tangential force reads
\begin{equation} \boldsymbol{F}_\parallel{=} \int_{\mathbb{R}^2} \left({-}p(\boldsymbol{r},t)\,\frac{\boldsymbol{r}}{a} - \eta \partial_z \boldsymbol{v} \right)_{z = h_0(r,t)}\textrm{d}^2\boldsymbol{r}. \end{equation}
Using symmetry arguments, we can show that the tangential force is directed along 
$x$, i.e. 
$\boldsymbol {F}_\parallel =F_x \boldsymbol {e}_x$. At small 
$\kappa$, we further expand it as 
$F_x \simeq F_{x,0} + \kappa F_{x,1}$, where 
$F_{x,0}$ is the viscous drag force applied on a sphere near a rigid plane wall, and 
$\kappa F_{x,1}$ is the first-order EHD correction. The zeroth-order term cannot be evaluated using the lubrication model introduced in the previous section, because the integral in (3.22) diverges, as the shear term 
$\eta \partial _z \boldsymbol {v}$ scales as 
$\sim r^{-2}$ at large 
$r$. An exact calculation has been performed using bispherical coordinates and provides a solution in the form of a series expansion (O'Neill Reference O'Neill1964). Asymptotic-matching methods have also been employed in order to get the asymptotic behaviour at small 
$d/a$ (Goldman, Cox & Brenner Reference Goldman, Cox and Brenner1967; O'Neill & Stewartson Reference O'Neill and Stewartson1967), which reads 
$F_{x,0} \approx 6 {\rm \pi}\eta a u (\frac {8}{15}\log ({d}/{a}) - 0.95429 )$ (see Chaoui & Feuillebois (Reference Chaoui and Feuillebois2003) for a high-precision expansion). We note that the sphere's normal velocity does not contribute to the zeroth-order tangential force, as expected by symmetry.
The first-order EHD correction can be computed with the present model, as the correction pressure field and shear stress scale as 
${\sim }r^{-5}$, at large 
$r$. It reads
\begin{align} F_{x,1} ={}& 2{\rm \pi} \eta u^* a \int_0^\infty \left[{-}2 R P_1^{(1)} - \frac{H_0}{2} \left(\partial_R P_1^{(1)} + \frac{P_1^{(1)}}{R} \right) \right.\nonumber\\ &\left.{}-\frac{H_1^{(1)}}{2}\,\partial_R P_0^{(0)} - \frac{H_1^{(0)}}{2} \left(\partial_R P_0^{(1)} + \frac{P_0^{(1)}}{R} \right) +2 \frac{U H_1^{(0)}}{H_0^2}\right] R\,\textrm{d}R, \end{align}
where 
$H_1^{(i)}$ is the amplitude of the 
$i$th mode in the angular-mode decomposition of 
$H_1$. Evaluating the latter integral numerically, we find
\begin{align} F_x \approx{}& 6 {\rm \pi}\eta a u \left(\frac{8}{15}\log\left(\frac{d}{a}\right) - 0.95429 \right) -10.884\, \frac{\eta^2 u \dot{d}(\lambda+2\mu)}{\mu(\lambda+\mu)} \left( \frac{a}{d}\right)^{5/2} \nonumber\\ &{}+ 0.98661\,\frac{\eta^2 \dot{u}a(\lambda+2\mu)}{\mu(\lambda+\mu)}\left( \frac{a}{d}\right)^{3/2}. \end{align}
We stress that the latter equation is not an exact truncated expansion. In Appendix B, we use again the Lorentz reciprocal theorem in order to compute the first-order EHD force, and we obtain the following analytical expressions for the coefficients of (3.24): 
$-({3177 {\rm \pi}^3}/{6400 \sqrt {2}}) \simeq -10.884$ and 
${9 {\rm \pi}^3}/{200 \sqrt {2}} \simeq 0.98661$, respectively.
The torque exerted by the fluid on the sphere, with respect to its centre of mass, is given by
\begin{equation} \boldsymbol{T}= \int_{\mathcal{S}_0} a\boldsymbol{n} \times (\boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma}) \, \mathrm{d}s. \end{equation}
The latter is directed along the 
$y$-direction for symmetry reasons, i.e. 
$\boldsymbol {T} =T_y \boldsymbol {e}_y$. At small 
$\kappa$, we further expand it as 
$T_{y} \simeq T_{y,0} + \kappa T_{y,1}$. For the same reason as with the the viscous drag force near a rigid wall, the viscous torque near a rigid wall cannot be computed within the lubrication model. Using asymptotic-matching methods (O'Neill & Stewartson Reference O'Neill and Stewartson1967; Chaoui & Feuillebois Reference Chaoui and Feuillebois2003), it is found to be 
$T_{y,0} \approx 8{\rm \pi} \eta u a^2 (-\frac {1}{10}\log (\frac {d}/{a}) - 0.19296 )$. In contrast, the first-order EHD correction can be computed with the present model, and reads
\begin{align} T_{y,1} ={}& -2\eta u^* a^2 {\rm \pi}\int_0^\infty \left[\frac{H_0}{2} \left(\partial_R P_1^{(1)} + \frac{P_1^{(1)}}{R} \right) + \frac{H_1^{(1)}}{2}\,\partial_R P_0^{(0)} \right.\nonumber\\ &\left.{}+ \frac{H_1^{(0)}}{2}\left(\partial_R P_0^{(1)} + \frac{P_0^{(1)}}{R} \right) +2 \frac{U H_1^{(0)}}{H_0^2}\right]R \, \textrm{d}R. \end{align}Evaluating the latter integral numerically, we find
\begin{align} T_y \approx{}& 8{\rm \pi} \eta u a^2 \left(-\frac{1}{10}\log\left(\frac{d}{a}\right) - 0.19296 \right) + 10.884\, \frac{\eta^2 u a\dot{d}(\lambda+2\mu)}{\mu(\lambda+\mu)} \left( \frac{a}{d}\right)^{5/2} \nonumber\\ &{}-0.98661\,\frac{\eta^2 \dot{u}a^2(\lambda+2\mu)}{\mu(\lambda+\mu)}\left( \frac{a}{d}\right)^{3/2}. \end{align}
The EHD torque in (3.27) is thus the same as the EHD tangential force in (3.24), up to a dimensional prefactor 
$-a$, as already observed for the EHD interactions of a rigid cylinder near a soft surface (Salez & Mahadevan Reference Salez and Mahadevan2015).
So far, we have focused on the particular case of a semi-infinite elastic material. In Appendices C and D, we apply the same soft-lubrication approach to other elastic models describing thin substrates, which are also widespread in practice. We find similar expressions, but with different numerical prefactors and scalings with the sphere-wall distance.
4. Rotation
We now add the rotation of the sphere, with angular velocity 
$\boldsymbol {\varOmega }(t)$ in the 
$xy$-plane (see figure 1), to the previous translational motion. We define 
$\beta$ as the angle between 
$\boldsymbol {\varOmega }$ and the 
$x$-axis. We stress that 
$\boldsymbol {\varOmega }$ is not necessarily orthogonal (i.e. 
$\beta = {\rm \pi}/2$) to the translation velocity. We discard the rotation along the 
$z$-axis (e.g. for a spinner), because it does not induce any soft-lubrication coupling. Finally, the system is equivalent to a purely rotating sphere with angular velocity 
$\boldsymbol {\varOmega }(t)$, near a wall translating with a 
$-\boldsymbol {u}(t)$ velocity. In such a framework, the fluid velocity field at the sphere surface is 
$\boldsymbol {v} = -\boldsymbol {\varOmega } \times a\boldsymbol {n}$, and thus 
$\boldsymbol {v} \simeq -\boldsymbol {\varOmega } \times a\boldsymbol {e}_z$. All together, the fluid velocity field is modified as
\begin{align} \boldsymbol{v}(\boldsymbol{r},z,t) ={}& \frac{\boldsymbol{\nabla} p(\boldsymbol{r},t)}{2\eta} (z - h_0(r,t))(z-\delta(\boldsymbol{r},t)) - \boldsymbol{u}(t)\,\frac{h_0(r,t)-z}{h_0(r,t)- \delta(\boldsymbol{r},t)} \nonumber\\ &{}+ a\,\boldsymbol{\varOmega}(t)\times \boldsymbol{e}_z\, \frac{z-\delta(\boldsymbol{r},t)}{h_0(r,t)-\delta(\boldsymbol{r},t)}, \end{align}and the Reynolds equation becomes
\begin{equation} \partial_t h(\boldsymbol{r}, t) = \boldsymbol{\nabla} \boldsymbol{\cdot} \left(\frac{h^3(\boldsymbol{r}, t)}{12\eta}\,\boldsymbol{\nabla} p(\boldsymbol{r}, t) + \frac{h(\boldsymbol{r}, t)}{2}\left[\underbrace{\boldsymbol{u}(t) -a\,\boldsymbol{\varOmega}(t)\times \boldsymbol{e}_z}_{\tilde{\boldsymbol{u}}}\right]\right). \end{equation}
The problem is thus formally equivalent to the one of a sphere that is purely translating with effective velocity 
$\boldsymbol {\tilde {u}} (t)= \boldsymbol {u}(t) -a\,\boldsymbol {\varOmega }(t)\times \boldsymbol {e}_z$. Therefore, we can directly apply the results from the previous sections, and write all the forces and torque exerted on the sphere, as
\begin{align} F_z ={}& -\frac{6{\rm \pi}\eta a^2\dot{d}}{d} + \frac{243 {\rm \pi}^3}{12800 \sqrt{2}}\, \frac{\eta^2 \lvert \boldsymbol{u} -a\boldsymbol{\varOmega}\times \boldsymbol{e}_z \rvert^2}{\mu} \left( \frac{a}{d}\right)^{5/2} - \frac{3915 {\rm \pi}^3}{2048 \sqrt{2}}\,\frac{\eta^2 \dot{d}^2}{\mu} \left(\frac{a}{d}\right)^{7/2} \nonumber\\ &{}+ \frac{27{\rm \pi}^3}{32\sqrt{2}}\,\frac{\eta^2 \ddot{d} a}{\mu} \left( \frac{a}{d}\right)^{5/2}, \end{align}
\begin{align} \boldsymbol{F}_\parallel{=}{}& 6 {\rm \pi}\eta a \boldsymbol{u} \left[\frac{8}{15}\log \left(\frac{d}{a} \right) - 0.95429 \right] +6{\rm \pi}\eta a^2 \boldsymbol{e}_z \times \boldsymbol{\varOmega} \left[\frac{2}{15} \log \left( \frac{d}{a} \right) + 0.25725 \right] \nonumber\\ &{}-\frac{3177 {\rm \pi}^3}{6400 \sqrt{2}}\,\frac{\eta^2 (\boldsymbol{u} -a\boldsymbol{\varOmega}\times \boldsymbol{e}_z) \dot{d}}{\mu} \left( \frac{a}{d}\right)^{5/2} +\frac{9 {\rm \pi}^3}{200 \sqrt{2}}\, \frac{\eta^2 (\boldsymbol{\dot{u}} -a\boldsymbol{\dot{\varOmega}} \times \boldsymbol{e}_z) a}{\mu} \left( \frac{a}{d}\right)^{3/2} \end{align}and
\begin{align} \boldsymbol{T}_\parallel{=}{}& 8 {\rm \pi}\eta a^2 \boldsymbol{e}_z\times \boldsymbol{u} \left[-\frac{1}{10}\log \left(\frac{d}{a} \right) - 0.19296 \right] +8{\rm \pi}\eta a^3 \boldsymbol{\varOmega} \left[\frac{2}{5} \log \left( \frac{d}{a} \right) - 0.37085\right] \nonumber\\ &{}+\frac{3177 {\rm \pi}^3}{6400 \sqrt{2}}\,\frac{\eta^2 (\boldsymbol{u} -a\boldsymbol{\varOmega} \times \boldsymbol{e}_z) a\dot{d}}{\mu} \left( \frac{a}{d}\right)^{5/2} -\frac{9 {\rm \pi}^3}{200 \sqrt{2}}\, \frac{\eta^2(\boldsymbol{\dot{u}} -a\boldsymbol{\dot{\varOmega}}\times \boldsymbol{e}_z)a^2}{\mu} \left( \frac{a}{d}\right)^{3/2}, \end{align}
where we have invoked the force and torque induced by the rotation of a sphere near a rigid wall (Goldman et al. Reference Goldman, Cox and Brenner1967; Urzay Reference Urzay2010) and where the analytical prefactors are computed in Appendices A and B. We stress that the expressions of the EHD forces and torque for a sphere purely translating near thin elastic substrates, as derived in Appendices C and D, can be generalized to further include the sphere's rotation by similarly following the transformation 
$\boldsymbol {u}(t) \rightarrow \boldsymbol {u}(t) -a\,\boldsymbol {\varOmega }(t)\times \boldsymbol {e}_z$.
5. Conclusion
We developed a soft-lubrication model in order to compute the EHD interactions exerted on an immersed sphere undergoing both translational and rotational motions near various types of elastic walls. The deformation of the surface was assumed to be small, which allowed us to employ a perturbation analysis in order to obtain the leading-order EHD forces and torque. The obtained interaction matrix exhibits a form that is qualitatively similar to the one found for a two-dimensional cylinder moving near a thin compressible substrate (Salez & Mahadevan Reference Salez and Mahadevan2015). In both cases, the EHD coupling is nonlinear and generates quadratic terms in the sphere velocity, thus breaking the time-reversal symmetry of the Stokes equations. In addition, original inertial-like terms proportional to the acceleration of the sphere are found – despite the assumption of steady flows. Therefore, while the quantitative details such as numerical prefactors and exponents differ in three dimensions and when using more realistic constitutive elastic responses, we expect that the typical zoology of trajectories identified previously (Salez & Mahadevan Reference Salez and Mahadevan2015) will also hold for spherical objects – and will even be extended with the added degree of freedom. As such, the asymptotic predictions obtained here may open new perspectives in colloidal science and biophysics, through the understanding and control of the emerging interactions within soft confinement or assemblies.
Acknowledgements
We thank A. Maali, Z. Zhang, H. Stone and B. Rallabandi for stimulating discussions.
Funding
The authors acknowledge funding from the Agence Nationale de la Recherche (ANR-21-ERCC-0010-01 EMetBrown) and from the Jean Langlois foundation.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Lorentz reciprocal theorem: normal force
In this appendix, we compute the first-order normal EHD force using the Lorentz reciprocal theorem for Stokes flows (Rallabandi et al. Reference Rallabandi, Saintyves, Jules, Salez, Schönecker, Mahadevan and Stone2017; Daddi-Moussa-Ider et al. Reference Daddi-Moussa-Ider, Rallabandi, Gekle and Stone2018; Rallabandi et al. Reference Rallabandi, Oppenheimer, Zion and Stone2018; Masoud & Stone Reference Masoud and Stone2019), in order to recover analytically the numerical prefactors obtained in the main text. To do so, we introduce the model problem of a sphere moving in a viscous fluid and towards an immobile, rigid, planar surface. We note 
$\boldsymbol {\hat {V}}_\perp = -\hat {V}_\perp \boldsymbol {e}_z$, the velocity at the particle surface 
$\mathcal {S}_0$, and we assume a no-slip boundary condition at the undeformed wall surface 
$\mathcal {S}_{w}$ located at 
$z = 0$ (see figure 1). The viscous stress and velocity fields of the model problem follow the steady, incompressible Stokes equations 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\hat {\sigma }}_\perp = \boldsymbol {0}$ and 
$\boldsymbol {\nabla } \boldsymbol {\cdot } \boldsymbol {\hat {v}}_\perp =0$, and we use the lubrication approximation. In this framework, the stress tensor is 
$\boldsymbol {\hat {\sigma }}_\perp \simeq -\hat {p}_\perp \boldsymbol {I} + \eta \boldsymbol {e}_z \partial _z \boldsymbol {\hat {v}}_\perp$, with
\begin{equation} \hat{p}_\perp(\boldsymbol{r}) = \frac{3\eta \hat{V}_\perp a}{\hat{h}^2(\boldsymbol{r})},\quad \boldsymbol{\hat{v}}_\perp(\boldsymbol{r},z) = \frac{\boldsymbol{\nabla}\hat{p}_\perp (\boldsymbol{r})}{2\eta}\,z(z-\hat{h}(\boldsymbol{r})),\quad \hat{h}(\boldsymbol{r}) = d + \frac{\boldsymbol{r}^2}{2a}. \end{equation}The Lorentz reciprocal theorem states that
\begin{equation} \int_\mathcal{S} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\sigma} \boldsymbol{\cdot} \boldsymbol{\hat{v}}_\perp \, \mathrm{d}s = \int_\mathcal{S} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\hat{\sigma}}_\perp \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}s, \end{equation}
where 
$\mathcal {S} = \mathcal {S}_0 + \mathcal {S}_\text {w} + \mathcal {S}_\infty$ is the total surface bounding the flow, and 
$\mathcal {S}_\infty$ is the surface located at 
$\boldsymbol {r} \rightarrow \infty$. The latter does not contribute here. Using the boundary conditions for the model problem, we get
\begin{equation} \boldsymbol{\hat{V}}_\perp\boldsymbol{\cdot} \boldsymbol{F} ={-}\hat{V}_\perp F_z = \int_{\mathcal{S}} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\hat{\sigma}}_\perp \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}s. \end{equation}
To find the force exerted on the sphere in the real problem, one needs to specify the boundary conditions for the real velocity field. Here, we assume that the sphere does not rotate, and we describe the flow in the translating reference frame of the particle. The no-slip boundary condition thus reads 
$\boldsymbol {v} = \boldsymbol {0}$ on 
$\mathcal {S}_0$. We further assume a small deformation of the wall, so that the velocity field at the undeformed wall surface can be obtained using the Taylor expansion:
\begin{align} \boldsymbol{v}\vert_{z = 0} &= \boldsymbol{v}\vert_{z = \delta} - \delta \partial_z \boldsymbol{v}_0\vert_{z = 0} \nonumber\\ & ={-}u\boldsymbol{e}_x - \dot{d}\boldsymbol{e}_z + (\partial_t - u\partial_x ) \delta \boldsymbol{e}_z- \delta \partial_z \boldsymbol{v}_0\vert_{z = 0}, \end{align}
where 
$\boldsymbol {v}_0$ is the zeroth-order velocity field near a rigid surface. Using results from the main text, we find
\begin{equation} \partial_z \boldsymbol{v}_0\vert_{z = 0} ={-}\frac{3\dot{d}r}{\left(d+\dfrac{r^2}{2a}\right)^2}\, \boldsymbol{e}_r + \frac{2u}{5\left(d+\dfrac{r^2}{2a}\right)} \left(\left(7-\frac{6d}{d+\dfrac{r^2}{2a}}\right) \cos\theta \boldsymbol{e}_r - \sin\theta\boldsymbol{e}_\theta \right). \end{equation}Combining (A1a–c) and (A4), we find the normal force:
\begin{equation} F_z = \frac{1}{\hat{V}_\perp} \int_{\mathbb{R}^2} \left(\hat{p}_\perp (-\dot{d} + \partial_t \delta - u\partial_x \delta)+\eta \partial_z \boldsymbol{\hat{v}}_\perp\vert_{z = 0}\cdot \partial_z \boldsymbol{v}_0 \vert_{z = 0} \delta \right)\textrm{d}\boldsymbol{r}. \end{equation}After some algebra, and computing the integral in Fourier space, we recover the same expression as in (3.21), which reads
\begin{align} F_z ={}& -\frac{6{\rm \pi} \eta a^2 \dot{d}}{d} + A\,\frac{\eta^2 u^2(\lambda+2\mu)}{\mu(\lambda+\mu)} \left( \frac{a}{d}\right)^{5/2} - B\,\frac{\eta^2 \dot{d}^2(\lambda+2\mu)}{\mu(\lambda+\mu)} \left( \frac{a}{d}\right)^{7/2} \nonumber\\ &{}+ C\,\frac{\eta^2 \ddot{d} a(\lambda+2\mu)}{\mu(\lambda+\mu)}\left( \frac{a}{d}\right)^{5/2}, \end{align}
where the numerical coefficients 
$A, B, C$ can be found analytically as
$$\begin{gather} A = \frac{9{\rm \pi}}{25\sqrt{2}} \int_0^\infty k^2\,K_0(k)\left({-}2K_1(k) + k\,K_2(k) \right) k \,\textrm{d}k = \frac{243 {\rm \pi}^3}{12\,800 \sqrt{2}}, \end{gather}$$
$$\begin{gather}B =9{\rm \pi}\sqrt{2} \int_0^\infty k^2\, K_1(k)\left(K_2(k)-\frac{k\,K_3(k)}{8}\right) k \,\textrm{d}k = \frac{3915 {\rm \pi}^3}{2048 \sqrt{2}}, \end{gather}$$
$$\begin{gather}C = \frac{9\sqrt{2}{\rm \pi}}{2} \int_0^\infty k^2\,K_1^2(k) \, \textrm{d}k = \frac{27{\rm \pi}^3}{32\sqrt{2}}, \end{gather}$$
and where 
$K_i$ is the modified Bessel function of the second kind of order 
$i$ (Abramowitz & Stegun Reference Abramowitz and Stegun1964).
Appendix B. Lorentz reciprocal theorem: tangential force
In order to compute the tangential force acting on the particle, we apply the Lorentz reciprocal theorem, but we introduce a different model problem with respect to the previous section. We consider a sphere translating parallel to a rigid immobile substrate, with a velocity 
$\hat {V}_\parallel$ along the 
$x$-axis, and no-slip boundary conditions at both the sphere and substrate surfaces. The velocity and stress fields are denoted 
$\boldsymbol {\hat {\sigma }}_\parallel$ and 
$\boldsymbol {\hat {v}}_\parallel$, respectively, and are solutions of the Stokes equations. The lubrication approximation is used here. The solution reads
\begin{equation} \hat{p}_\parallel(\boldsymbol{r}) = \frac{6\eta \hat{V}_\parallel r\cos\theta}{5\hat{h}^2(\boldsymbol{r})},\quad \boldsymbol{\hat{v}}_\parallel(\boldsymbol{r},z) = \frac{\boldsymbol{\nabla}\hat{p}_\parallel(\boldsymbol{r})}{2\eta}\,z(z-\hat{h}(\boldsymbol{r})) + \boldsymbol{\hat{V}}_\parallel\,\frac{z}{\hat{h}(\boldsymbol{r})}, \end{equation}as shown in the main text. The Lorentz reciprocal theorem leads to
\begin{equation} \boldsymbol{\hat{V}}_\parallel \boldsymbol{\cdot} \boldsymbol{F} = \hat{V}_\parallel F_x = \int_{\mathcal{S}} \boldsymbol{n} \boldsymbol{\cdot} \boldsymbol{\hat{\sigma}}_\parallel \boldsymbol{\cdot} \boldsymbol{v} \, \mathrm{d}s. \end{equation}
Using the lubrication expression of the stress tensor of the model problem, 
$\boldsymbol {\hat {\sigma }}_\parallel \simeq -\hat {p}_\parallel \boldsymbol {I} + \eta \boldsymbol {e}_z \partial _z \boldsymbol {\hat {v}}_\parallel$, we get an expression for the tangential force as
\begin{equation} F_x = \frac{1}{\hat{V}_\parallel} \int_{\mathbb{R}^2} \left[ -\eta \partial_z \boldsymbol{\hat{v}}_\parallel{\cdot} u(t)\, \boldsymbol{e}_x -\hat{p}_\parallel(\partial_t - u(t)\,\partial_x)\delta - \eta (\partial_z \boldsymbol{\hat{v}}_\parallel{\cdot} \partial_z\boldsymbol{v}_0) \delta\right] \textrm{d}\boldsymbol{r}. \end{equation}
For the same reason as the one invoked in the main text, the zeroth-order tangential drag force (i.e. the integral of 
$-\eta \partial _z \boldsymbol {\hat {v}}_\parallel \cdot u(t) \boldsymbol {e}_x$ in (B4)) cannot be computed here as the integral diverges within the lubrication approximation. In contrast, the first-order EHD force is well defined in the lubrication framework and can be computed in Fourier space using Parseval's theorem, leading to
\begin{equation} \kappa F_{x,1} ={-} \frac{3177 {\rm \pi}^3}{6400 \sqrt{2}}\, \frac{\eta^2 u \dot{d}(\lambda+2\mu)}{\mu(\lambda+\mu)} \left(\frac{a}{d}\right)^{5/2} + \frac{9 {\rm \pi}^3}{200 \sqrt{2}}\, \frac{\eta^2 \dot{u}a(\lambda+2\mu)}{\mu(\lambda+\mu)}\left( \frac{a}{d}\right)^{3/2}. \end{equation}Appendix C. Thin compressible substrate
In this appendix, we derive the EHD interactions exerted on a sphere immersed in a viscous fluid and near a thin compressible substrate of thickness 
$h_{sub}$. The deformation field follows the Winkler foundation
\begin{equation} \delta(\boldsymbol{r}, t) ={-}\frac{h_{sub}}{(2\mu + \lambda)}\,p(\boldsymbol{r},t), \end{equation}
which is valid for substrates of thickness smaller than the typical extent of the pressure field, namely the hydrodynamic radius 
$\sqrt {2ad}$ (Leroy & Charlaix Reference Leroy and Charlaix2011; Chandler & Vella Reference Chandler and Vella2020; Kargar-Estahbanati & Rallabandi Reference Kargar-Estahbanati and Rallabandi2021). We perform the same asymptotic expansion as the one in the main text, defining the Winkler dimensionless compliance as (Salez & Mahadevan Reference Salez and Mahadevan2015)
\begin{equation} \kappa^{W} = \frac{\sqrt{2}h_{sub}\eta u^* a^{1/2}}{d^{*5/2} (2\mu + \lambda)}. \end{equation}The first-order substrate deformation, or equivalently here the zeroth-order pressure, reads
\begin{equation} H_1^{W}(\boldsymbol{R}, T) = P_0(\boldsymbol{R}, T) = \frac{3 \dot{D}}{2 (D+R^2)^2} + \frac{6 R U \cos\theta}{5(D+ R^2)^2}. \end{equation}The first-order deformation fields are plotted in figure 5(a,b) for tangential and normal motions of the sphere, respectively. The deformation exhibits the same structure as the one in figure 2 for semi-infinite substrates, but the lateral extent of the deformation is narrower. This is expected as the deformation response induced by a given applied pressure is local for a thin compressible layer (see (C1)), while semi-infinite substrates display a non-local response due to the convolution of the pressure with their Green's function (see (2.3)). The first-order pressure correction follows the same type of equation as in the main text:
\begin{equation} \mathcal{L}\cdot P_1^{W} = F_0^{W}(R,T) + F_1^{W}(R,T) \cos\theta + F_2^{W}(R,T) \cos2\theta, \end{equation}with
\begin{align} F_0^{W}(R,T) ={}& -\frac{144 R^2 U^2}{25 (D+R^2)^7} [ D^2 - 6DR^2 +R^4 ] \nonumber\\ &{}+ \frac{18R^2 \dot{D}^2}{ (D+R^2)^7} [ 5D-4R^2 ] -\frac{18R^2 \ddot{D}}{ (D+R^2)^5} \end{align}and
\begin{equation} F_1^{W}(R,T) = \frac{216 R^3 U\dot{D} }{5 (D+R^2)^7} [{-}5D + R^2 ] + \frac{72R \dot{U}}{5 (D+R^2)^5}. \end{equation}
We note that 
$F_2^{W}$ does not contribute for the forces and torque. The isotropic component of the pressure can be found analytically, using polynomial fractions, as
\begin{align} P_1^{{W},(0)}(R,T) =\frac{9}{125}\,\frac{7-5Y^2}{(1+Y^2)^5}\,\frac{U^2}{D^4} - \frac{3}{40}\,\frac{71+55Y^2+30Y^4}{(1+Y^2)^5}\,\frac{\dot{D}^2}{D^5} +\frac{3}{2}\,\frac{1}{(1+Y^2)^3}\, \frac{\ddot{D}}{D^4}, \end{align}
where 
$Y = R/D^{1/2}$ is the self-similar variable. However, the first angular component of the pressure does not exhibit such an analytical solution, and is thus found by numerical integration of two scaling functions. Its general expression reads
\begin{equation} P_1^{{W},(1)}(R,T) =\frac{U\dot{D}}{D^{9/2}}\,\phi_{U\dot{D}}^{W} \left(\frac{R}{D^{1/2}}\right)+\frac{\dot{U}}{D^{7/2}}\,\phi_{\dot{U}}^{W}\left(\frac{R}{D^{1/2}}\right). \end{equation}Following the same calculation as in the main text, we find the normal force as
\begin{align} F_z^{W} ={}& -\frac{6{\rm \pi} \eta a^2 \dot{d}}{d} + \frac{48{\rm \pi}}{125}\, \frac{\eta^2 u^2 h_{sub}}{a(2\mu +\lambda)} \left(\frac{a}{d}\right)^3 - \frac{72{\rm \pi}}{5}\,\frac{\eta^2 \dot{d}^2 h_{sub}}{a(2\mu +\lambda)} \left(\frac{a}{d}\right)^4 \nonumber\\ &{}+ \frac{6 {\rm \pi}\eta^2 \ddot{d} h_{sub}}{(2\mu +\lambda)} \left(\frac{a}{d}\right)^3. \end{align}
We stress that the prefactors 
$48{\rm \pi} /125$ and 
$6{\rm \pi}$ are in agreement with the results in Urzay et al. (Reference Urzay, Llewellyn Smith and Glover2007) and Leroy & Charlaix (Reference Leroy and Charlaix2011), respectively. Similarly, the force along 
$x$ reads
\begin{align} F_{x}^{W} ={}& 6 {\rm \pi}\eta a u \left(\frac{8}{15}\log\left(\frac{d}{a}\right) - 0.95429 \right) - \frac{484 {\rm \pi}}{125}\,\frac{\eta^2u\dot{d}h_{sub}}{a(2\mu+\lambda)} \left(\frac{a}{d}\right)^3 \nonumber\\ &{}+\frac{12 {\rm \pi}}{25}\,\frac{\eta^2\dot{u}h_{sub}}{(2\mu+\lambda)} \left(\frac{a}{d}\right)^2, \end{align}The torque can be evaluated as well, and reads
\begin{align} T_{y}^{W} ={}& 8{\rm \pi} \eta u a^2 \left(-\frac{1}{10}\log\left(\frac{d}{a}\right) - 0.19296 \right) + \frac{484 {\rm \pi}}{125}\,\frac{\eta^2u\dot{d}h_{sub}}{(2\mu+\lambda)} \left(\frac{a}{d}\right)^3 \nonumber\\ &{}-\frac{12 {\rm \pi}}{25}\,\frac{\eta^2\dot{u}ah_{sub}}{(2\mu+\lambda)} \left(\frac{a}{d}\right)^2. \end{align}
Figure 5. Dimensionless deformation fields at the free surface of two soft substrates, for a sphere placed at a unit distance 
$D = 1$ and for two modes of motion. In (a,b), the substrate's mechanical response follows the Winkler foundation (see (C1)). In (c,d), the substrate's mechanical response is the one of a thin incompressible layer (see (D1)).
All the prefactors for the EHD corrections of the tangential force and torque have been found analytically with the Lorentz reciprocal theorem (see Appendix B). Finally, following the approach in the main text, it is straightforward to generalize (C9), (C10) and (C11) in order to incorporate rotation.
Appendix D. Thin incompressible substrate
In this appendix, we suppose that the substrate of thickness 
$h_{sub}$ is incompressible, i.e. of Poisson ratio 
$\nu = 1/2$, which means that the first Lamé coefficient 
$\lambda$ is infinite. In this situation, the Winkler foundation is not valid. Instead, the mechanical response of a thin substrate follows the relation (Leroy & Charlaix Reference Leroy and Charlaix2011; Chandler & Vella Reference Chandler and Vella2020)
\begin{equation} \delta(\boldsymbol{r}, t) = \frac{h^3_{sub}}{3\mu}\,\nabla^2 p(\boldsymbol{r},t), \end{equation}
where 
$\nabla ^2$ denotes the two-dimensional Laplacian operator in the 
$(x,y)$-plane. We perform the same asymptotic expansion as in the main text, defining the thin-incompressible dimensionless compliance as
\begin{equation} \kappa^{t-i} = \frac{\eta u^* h_{sub}^3 }{ 3\sqrt{2}\mu d^{*7/2} a^{1/2} }. \end{equation}The first-order substrate deformation reads
\begin{equation} H_1^{t-i}(\boldsymbol{R}, T) ={-}\nabla^2 P_0(\boldsymbol{R}, T) ={-}\frac{12 \dot{D}(D-2R^2)}{ (D+R^2)^4} + \frac{48 R U (2D-R^2)\cos\theta}{5(D+ R^2)^4}. \end{equation}The deformation fields are plotted in figure 5(c,d), for tangential and normal motions of the sphere, respectively. The first-order pressure correction follows the same type of equation as in the main text:
\begin{equation} \mathcal{L}\cdot P_1^{t-i} = F_0^{t-i}(R,T) + F_1^{t-i}(R,T) \cos\theta + F_2^{t-i}(R,T) \cos2\theta, \end{equation}with
\begin{align} F_0^{t-i}(R,T) ={}& \frac{1152 R^2 U^2 \left(R^2-2 D\right) \left({+}2 D^2-11 R^2 D+2 R^4\right)}{25 \left(D+R^2\right)^9} \nonumber\\ &{}+\frac{432 R^2 \left(2 D \left(D-5 R^2\right)+3 R^4\right) \dot{D}^2}{\left(D+R^2\right)^9} +\frac{144 R^2 \left(2 R^2-D\right) \ddot{D}}{\left(D+R^2\right)^7} \end{align}and
\begin{equation} F_1^{t-i}(R,T) ={-}\frac{2592 R^3 \dot{D}U \left(7 D^2-12 R^2 D+R^4\right) }{5 \left(D+R^2\right)^9} -\frac{576 R^3\dot{U} \left({-}2 D+R^2\right) }{5 \left(D+R^2\right)^7}. \end{equation}
We note that 
$F_2^{t-i}$ does not contribute to the forces and torque. The isotropic component of the pressure can be found analytically, using polynomial fractions, as
\begin{align} P_1^{{t-i},(0)}(R,T) =\frac{288 \left(7 Y^4-21 Y^2+17\right)}{875 \left(1+Y^2\right)^7}\, \frac{U^2}{D^5} + \frac{126 Y^2-198}{7 \left(1+Y^2\right)^7}\,\frac{\dot{D}^2}{D^6} + \frac{36}{5 \left(1+Y^2\right)^5}\,\frac{\ddot{D}}{D^5}, \end{align}
where 
$Y = R/D^{1/2}$ is the self-similar variable. However, the first angular component of the pressure does not exhibit such an analytical solution, and is thus found by numerical integration of two scaling functions. Its general expression reads
\begin{equation} P_1^{{t-i},(1)}(R,T) =\frac{U\dot{D}}{D^{11/2}}\,\phi_{U\dot{D}}^{t-i} \left(\frac{R}{D^{1/2}}\right)+\frac{\dot{U}}{D^{9/2}}\,\phi_{\dot{U}}^{t-i} \left(\frac{R}{D^{1/2}}\right). \end{equation}Following the same calculation as in the main text, we find the normal force as
\begin{align} F_z^{t-i} ={}& -\frac{6{\rm \pi} \eta a^2 \dot{d}}{d} + \frac{432 {\rm \pi}}{875}\, \frac{\eta^2 u^2 h_{sub}^3}{a^3\mu} \left(\frac{a}{d}\right)^4 - \frac{192 {\rm \pi}}{35}\,\frac{\eta^2 \dot{d}^2 h_{sub}^3}{a^3\mu} \left(\frac{a}{d}\right)^5 \nonumber\\ &{}+ \frac{12 {\rm \pi}}{5}\,\frac{ \eta^2 \ddot{d} h_{sub}^3}{a^2\mu} \left(\frac{a}{d}\right)^4. \end{align}
We stress that the prefactor 
$12{\rm \pi} /5$ is consistent with the linear-response theory in Leroy & Charlaix (Reference Leroy and Charlaix2011). Similarly, the force along 
$x$ reads
\begin{align} F_{x}^{t-i} ={}& 6 {\rm \pi}\eta a u \left(\frac{8}{15}\log\left(\frac{d}{a}\right) - 0.95429 \right) - \frac{64 {\rm \pi}}{35}\,\frac{\eta^2u\dot{d}h_{sub}^3}{a^3\mu} \left(\frac{a}{d}\right)^4 \nonumber\\ &{}+\frac{32 {\rm \pi}}{125}\,\frac{\eta^2\dot{u}h_{sub}^3}{a^2\mu} \left(\frac{a}{d}\right)^3. \end{align}The torque can be evaluated as well, and reads
\begin{align} T_{y}^{t-i} ={}& 8{\rm \pi} \eta u a^2 \left(-\frac{1}{10}\log\left(\frac{d}{a}\right) - 0.19296 \right) + \frac{64 {\rm \pi}}{35}\,\frac{\eta^2u\dot{d}h_{sub}^3}{a^2\mu} \left(\frac{a}{d}\right)^4 \nonumber\\ &{}-\frac{32 {\rm \pi}}{125}\,\frac{\eta^2\dot{u}h_{sub}^3}{a \mu} \left(\frac{a}{d}\right)^3. \end{align}
Here again, the prefactors of the transverse force and torque are found using the Lorentz reciprocal theorem, as discussed above in Appendices A, B and C. We stress that the thin-incompressible limit is mathematically valid for strictly incompressible substrates, but that its range of application is limited in practice. Indeed, usual elastomers and gels, which are considered as almost incompressible, have a Poisson ratio close to 
$\nu \simeq 0.49$, and thus a tiny but finite compressibility. In recent studies (Saintyves et al. Reference Saintyves, Jules, Salez and Mahadevan2016; Rallabandi et al. Reference Rallabandi, Saintyves, Jules, Salez, Schönecker, Mahadevan and Stone2017; Saintyves et al. Reference Saintyves, Rallabandi, Jules, Ault, Salez, Schönecker, Stone and Mahadevan2020), it has been observed that the mechanical response of very thin incompressible elastic substrates is better described by the Winkler foundation (i.e. thin and compressible) than the thin-incompressible limit discussed here. This observation has then been established on solid theoretical grounds for the EHD lift (Chandler & Vella Reference Chandler and Vella2020), and is intimately rooted in the structure of the elastic Green's function. An empirical scaling, based on the numerical calculation of the EHD lift coefficient, has been derived subsequently (Kargar-Estahbanati & Rallabandi Reference Kargar-Estahbanati and Rallabandi2021) and suggests that the thin-incompressible model is valid for thicknesses comprised in the range 
$\frac {\sqrt {7}}{3}(1/2-\nu )^{1/2} \ll h_{sub}/\sqrt {2ad} \leq 0.12$. For 
$\nu = 0.49$, the lower bound of the latter range is 
$0.088$, which confirms that the validity window of the thin-incompressible model is limited.
