Regimes of Soft Lubrication

Elastohydrodynamic lubrication, or simply soft lubrication, refers to the motion of deformable objects near a boundary lubricated by a fluid, and is one of the key physical mechanisms to minimise friction and wear in natural and engineered systems. Hence it is of particular interest to relate the thickness of the lubricant layer to the entrainment (sliding/rolling) velocity, the mechanical loading exerted onto the contacting elements, and properties of the elastic boundary. In this work we provide an overview of the various regimes of soft lubrication for two-dimensional cylinders in lubricated contact with compliant walls. We discuss the limits of small and large entrainment velocity, which is equivalent to large and small elastic deformations, as the cylinder moves near thick or thin elastic layers. The analysis focusses on thin elastic coatings, both compressible and incompressible, for which analytical scaling laws are not yet available in the regime of large deformations. By analysing the elastohydrodynamic boundary layers that appear at the edge of the contact, we establish the missing scaling laws - including prefactors. As such, we offer a rather complete overview of physically relevant limits of soft lubrication.


Introduction
The introduction of a viscous liquid within the narrow gap between two moving bodies prevents solid-solid contact, reduces friction, minimizes wear, and allows precise motion control; a process known as lubrication. Lubrication has played a pivotal role in major engineering milestones, from revolutionising transportation of heavy objects on wooden sledges in ancient Egypt (Dowson 1998) to the efficient and reliable functioning of bearings in applications on many scales ranging up to the huge rolling bearings in todays wind turbines. The underlying hydrodynamic framework of fluid flow in thin gaps is in fact relevant to a range of phenomena across different length scales: motion of drops/bubbles in microfluidic devices (Bretherton 1961;Byun 2013), moving blood cells in capillaries (Fitz-Gerald 1969;Secomb et al. 1986), and biomechanics of synovial joints (Hou et al. 1992) to name a few. However, in most of the engineering applications as well as natural systems, one or both of the moving boundaries are deformable. As a result, lubrication flow is coupled to elasticity through the non-local soft lubrication equations.
Since the realization of the essential coupling between deformation and flow for film formation, the study of elasto-hydrodynamic lubrication in engineering (Tribology) has established as a separate field of research. Its importance has only increased with the current need of reduced material and lubricant in view of sustainability demands. More complex time-varying and extreme operating conditions, such as higher loads and higher temperatures, lead to increased deformation and thinner films which calls for better understanding and improved prediction capability for engineering and design. Only in very specific asymptotic cases analytic solutions could be derived, and numerical solution is often needed. Pioneering work was done by Higginson 1959, 1966) and Hamrock & Dowson (1977), where the complexity of the interaction between deformation and flow led to numerical stability problems at large deformations. Faster computers and improved numerical algorithms (Hamrock et al. 2004;Venner & Lubrecht 2000) allowed the study of realistic problems related to roller bearings, gears, cam-followers, and seals. Many studies aimed at the derivation of empirical formulas for film thickness prediction under steady conditions from the numerical solutions. Experimental validation was obtained from the further development of optical interferometry based techniques (Gohar & Cameron 1963;Cann et al. 1996). At present, complex cases of contact geometries with time varying loading conditions, roughness moving through the contact, contacts with very limited lubrication supply, grease lubricated contacts, and coated surfaces are considered, see Wang (2019). In spite of the plethora of problems studied and results presented the fundamental understanding of the physical phenomena and the appropriate scaling laws developed only relatively slowly. A recent overview is given by (Greenwood 2020).
A key feature of soft lubrication is the emergence of a non-inertial lift force, which facilitates lubrication by maintaining the gap between the moving surfaces (Sekimoto & Leibler 1993;Skotheim & Mahadevan 2004). For micron sized particles flowing in a channel with soft walls, this force causes radial migration of particles (Davies et al. 2018) and promotes particle sorting in microfluidic devices (Geislinger & Franke 2014). A precise measurement of this elastohydrodynamic force can function as a contact-free tool to probe the local rheology of soft materials (Leroy et al. 2012;Wang et al. 2015). For sufficiently soft elastomers, this force is strong enough to sustain the weight of millimeter sized metal cylinders which slide along a deformable wall (Saintyves et al. 2016), and may result in intricate trajectories (Salez & Mahadevan 2015;Rallabandi et al. 2017).
Naturally, it is of fundamental interest to explain how the gap thickness evolves from the coupling between fluid pressure and elastic deformation. The answer depends on details of the physical limit under consideration: Since deformation of a soft layer crucially depends on the geometry, and boundary conditions at the free surface, scaling laws relating the lift force to the lubrication layer are very different in different physical limits. Indeed, recent microfluidic experiments involving polymer brushes, cross linked elastomers, and ionic microgels as soft boundaries have probed elastohydrodynamic interactions in a few different limits, and experimentally found very distinct scaling laws (Davies et al. 2018;Zhang et al. 2020;Vialar et al. 2019).
Specifically, the elastic layer can either be thick or thin, the material can be compressible or incompressible, while the elastic deformation can be large or small compared to the lubricant thickness. A well-studied limit is when the soft layer is thin and compressible, so that deformation is proportional to the local pressure, while assuming deformations to be small (Skotheim & Mahadevan 2005;Urzay et al. 2007;Salez & Mahadevan 2015). This Figure 1. Schematic representation of a lubricated contact with a thin, compressible elastic layer, showing the relevant length scales of the problem: elastic layer thickness d, contact size a, the fluid film thickness h(x) and the elastic deformation u(x) defined in downward direction. One of the edges has been magnified to show the boundary layer width , over which lubrication forces smooth out the elastic deformation.
allows for a perturbative approach to solve the coupled elastohydrodynamic equations. However, most of the soft materials like elastomers and hydrogels are incompressible, and their response is described by non-local integral equations. Furthermore, if particles are squeezed past a soft layer, a Hertz-like contact forms, where deformation is large as compared to the gap thickness. In this case, the effect of lubrication is primarily dominant within a narrow boundary layer at inlet/outlet regions near the contact edges (Bissett 1989;Snoeijer et al. 2013). This lubricated Hertzian limit is not only relevant for classical tribological applications, but also plays a role in soft matter systems such as the slippage of soft microgel pastes (Meeker et al. 2004). Naturally, different scaling laws become relevant in these various limits, but an overview of how these different cases emergeand when these are applicable -is still lacking.
Here we present a unified overview of different regimes of soft lubrication. Our asymptotic analysis will focus on several important cases that were not considered previously, most notably the lubricated Hertz-like limit for cylinders on thin elastic layers on a rigid substrate. We also explore the transition to the small deformation limit, using numerical simulations. The paper is organised as follows. In section 2, we develop the theoretical and numerical framework and discuss the numerous length scales of the problem. Here we emphasise how the relative magnitudes of these length scales determine the different asymptotic regimes of soft lubrication. In section 3 we resolve in detail the lubrication on thin compressible layers in the limit of large deformation, and establish the scaling laws through similarity solutions. In the final section 4 this is complemented by the results for thin incompressible layers. Hence, we provide a complete overview of all regimes, summarised in a table in the concluding section, and discuss experimental consequences.

Problem formulation
We consider the two-dimensional problem of a rigid cylinder moving at constant velocity parallel to a deformable substrate, as sketched in fig. 1. Here we consider the translation of the cylinder, but we remark that the analysis is identical for the case that the cylinder is rotating (Pandey et al. 2016).
Importantly, the analysis is organised in terms of the length scales that appear in the problem ( fig. 1). The hierarchy of the length scales is a crucial aspect of the lubrication problem, as the relative sizes of these lengths determine the asymptotic limit. Hence, the organisation of lengths scale will determine the scaling laws that relate the load and velocity to the thickness of the lubrication layer.
A first aspect of the length scales is that we assume the cylinder radius, R, to be much larger than the width of the contact zone, a, so that the shape can locally be approximated as a parabola. The lubricant layer thickness, h(x), and the elastic deformation, u(x), are then related by ( fig. 1) Here c ≡ u(0) − h(0) is the indentation of the bottom of the cylinder, measured with respect to the undeformed elastic substrate. Another obvious length scale is the thickness d of the elastic layer. Finally, we anticipate the emergence of a dynamic length scale , sketched in the inset of fig. 1, which acts as a narrow boundary layer at the edges of the contact (Snoeijer et al. 2013). We will find that the relative magnitudes of a, d, , u, and the compressibility of the elastic layer (quantified by the Poisson ratio ν), will determine how the lubricant layer h scales with the sliding velocity. We remark that, apart from d, none of these lengths are known a priori, but follow from a consistent analysis of the problem. Now that the geometry is in place, we can turn to the mechanical equations. We start with the flow of liquid within the narrow gap, which is described by the Stokes equation, ∇p = η∇ 2 v, where p is pressure, η is viscosity and v is the velocity field. We consider a reference system where the cylinder is stationary, and the substrate moves with a velocity V to the right. Since both the film thickness and its gradient are small compared to the width of the contact, using no slip boundary conditions the velocity profile within the gap can be written in the lubrication limit, Incompressibility of the liquid requires that, at steady-state, the flux Q = u u−h v x dz remains constant. As a result, integration of (2.2) gives the Reynolds equation (Reynolds 1886), where h * = 2Q/V . Next, we couple the pressure p(x) to the elastic deformation u(x) of the solid. This is most conveniently done by transforming the problem to the Fourier domain (f ). Namely, the elastic deformation of the substrate in Fourier domain is given bỹ where the Green's function K(q) reads (Hannah 1951;Wang 2019), (2.5) In this expression ν is the Poisson ratio, G is the shear modulus, while d is the layer thickness. For a given pressure, the backward transform to u(x) is to be performed numerically. However, the above kernel can be simplified for various limiting cases, leading to analytical expressions of u(x) in real-space. We distinguish between two limits defined by the quantity qd, a ratio of the substrate thickness to the typical wavelength q of the deformation. In the short wavelength limit (qd 1) the substrate's response locally reduces to that of an elastic half space K(q) ∼ |q| −1 . In the opposite limit of large

Kernel expansion Deformation relation
Half space (qd 1) Thin layer (qd 1), Thin layer (qd 1), Table 1. Overview of asymptotes for the Fourier space Green's function K(q), given in (2.5). The two columns give the expansions of the Green's function and the corresponding pressure-deformation relations in real space.
The analysis can be closed by either imposing a vertical load exerted on the cylinder, or equivalently by fixing its vertical position. In our calculations we will impose the total vertical force L, expressed in terms of the fluid pressure, (2.6) The problem is defined by the system of three equations, (2.1, 2.3, 2.4), which describe the three unknown fields h(x), p(x) and u(x). Furthermore the constraint (2.6) enables us to find the vertical position c of the cylinder. Combined with the boundary conditions p(±∞) = 0 and u(±∞) = 0, the problem is fully defined.

Regimes
Since our goal is to identify and describe the various limits of the elastohydrodynamic problem, the relevant length scales need to be carefully compared. These length scales as defined in figure 1 are given by the contact size a, the elastic layer thickness d, the boundary layer width , the typical surface displacement u 0 , and the sought-after lubricant thickness h 0 . The reader can already consult table 2 for a complete overview of the various regimes.
Small versus large deformations. Two key limits of the problem come from the comparison between the deformation u 0 and the fluid film thickness h 0 . When the sliding velocity of the cylinder is high, the cylinder experiences a strong upward force due to viscous stresses that separates it from the substrate. This gives the small deformation limit, where h 0 u 0 (Skotheim & Mahadevan 2005;Pandey et al. 2016). By contrast, for relatively low sliding velocity the cylinder indents the substrate with only a thin layer of fluid in between the two, so that h 0 u 0 . This is what will be referred to as the large deformation limit (Bissett 1989;Snoeijer et al. 2013); remark that we still work in the context of linear elasticity, so that typical strains remain small.
Short versus long wavelengths. In the large deformation limit, a narrow lubrication boundary layer develops at the edge of the contact which regularizes the jump in pressure gradient caused by elasticity. Comparing the width of this boundary layer to the substrate thickness d, we can further distinguish two limits in the regime of large deformations. When d, the elastic problem is in the long wavelength limit, in which the thin-layer approximation (qd 1) is valid over the full width of the contact. In the short wavelength limit ( d), the thin elastic layer approximation does not hold at the edge of the contact, and substrate locally responds like a half space (qd 1) within the boundary layer.
Compressible versus incompressible elastic layers. For each of the three above mentioned limits, one can choose a suitable response function given in table 1. For a comparatively thin layer, d , one needs to make an additional distinction between a compressible and an incompressible material. This distinction is made since the leading order term of the Taylor expansion of the kernel scales with (1 − 2ν), which vanishes in the incompressible case (ν = 1 2 ) making the next order term dominant. In the rest of the paper, we develop solutions of the elastohydrodynamic equations and in table 2 we summarise these solutions through the scaling laws for the different regimes discussed above.

Non-dimensionalization
While various of the cases reported in table 2 are already available in the literature, we are not aware of any lubrication analysis of large deformations on thin elastic layers. Therefore, we will from now on consider lubricated contacts with thin elastic coatings, and we primarily focus on the case of compressible layers. We introduce dimensionless variables that are suitable for this specific case (in the long-wave limit of the elastic deformation), and identify the single dimensionless number that governs the solution.
We consider the geometric relation (2.1) along with the deformation on thin compressible layers as given in table 1. Combining these two equations we get, (2.7) The small deformation limit has been solved before in Skotheim & Mahadevan (2005), but here we consider the opposite limit. To non-dimensionalize the above equation for this limit we introduce the following scales, where variables with an overbar are dimensionless. From this, one infers a natural characteristic horizontal scale, the contact size a, as (2.10) Here we further introducedc = c/u 0 , based on the natural vertical scale (2.11) Using the relation between pressure and deformation (2.10), we find that Reynolds equation (2.3) can be written as which contains the single dimensionless parameter of the problem, This parameter can be interpreted as the dimensionless velocity. The lubrication problem for a compressible elastic layer in the long-wave limit, is thus described by a first order ODE for the flowing thicknessh(x), given by (2.12). Note, however, that there are two boundary conditions to be imposed, namely the pressure must vanish atx = ±∞, which can be expressed in terms ofh using (2.10). The two boundary conditions for the first order ODE (2.12) can indeed be satisfied by tuning the value ofh * . In addition, the coefficientc must be tuned such that (2.14) in order to impose the desired value of the load.

Numerical methods
To confirm the analytical results and to assess the crossover between the asymptotic regimes, we employ two numerical methods. The first method is numerically solving (2.12) in Wolfram Mathematica, along with the corresponding boundary conditions and constraints. We use a shooting method for a range of values forc, finding the pair ofh * and λ that allows to satisfy the boundary conditions,p(x → ±∞) = 0 and ∞ −∞p (x)dx = 4 3 . Typical examples are given in fig. 3. These accurately resolve the problem as long as the solution is self-consistent with the long-wave expansion.
However, we will find that in the limit of very small λ, the obtained solution is no longer consistent with the long-wave expansion (of the elastic solid) that underlies (2.12), we remark that the long-wave approximation for the liquid is satisfied for all cases considered. In this case, we turn to finite element simulations, where the two-dimensional plane strain problem is implemented using the library oomph-lib (Heil & Hazel 2006). A typical example of the numerical simulation is shown in figure 2. The compressible elastic layer is simulated by a non-linear solid model using the Neo-Hookean constitutive law with a Poisson ratio of ν = 0.3. Since the strains in this problem always remain small, the simulations reduce to linear elasticity and the choice of the non-linear rheology does not influence the result. The domain has a width W = 10 5 u 0 (|x| W/2) and height d = 30 u 0 , and consists of refineable quad elements. The cylinder has a radius of R = 10 6 u 0 , making for a contact size of a 1.4 × 10 3 u 0 . The fluid pressure is imposed on the top boundary of the elastic layer using surface elements that simulate the Reynolds equation (2.3). The pressure in the fluid is constrained to zero at both sides of the domain and the nodes at the bottom of the elastic layer are fixed in both the vertical and horizontal direction. The nodes in both the left and right side of the layer are fixed only in the horizontal direction and are subject to a no-shear condition in the vertical direction.

Thin compressible layer
In this section we present results for the lubricated contact with a thin, compressible, elastic substrate. We first illustrate the emergence of three distinct asymptotic regimes by presenting numerical results, both from the long-wave expansion of (2.12) and based on the finite element simulations. Then, we turn to a detailed asymptotic analysis of the two regimes that involve large deformations, respectively corresponding to long and short wavelengths. Figure 3 reports a number of film thickness profilesh = h/u 0 for lubricated contacts, for various values of the dimensionless sliding velocity λ as defined by (2.13). The results presented in the figure were obtained from numerical integration of (2.12), and also includes the "dry" solution, corresponding to the case λ = 0 where there is no motion. The profiles in fig. 3 clearly reveal a qualitative change in the shape of the contact, upon increasing the velocity. At small λ, the profiles exhibit a large elastic deformation where the cylinder pushes deep into the soft coating. The liquid profile closely resembles that of the "dry" solution, except near the inlet (x = −1) and the outlet (x = 1) regions where the effect of lubrication is prominent. At these two edges of the contact region, we encounter a narrow boundary layer of width , as sketched in the inset of fig. 1. These boundary layers will be analysed in detail below, as they are essential for computing the entrained liquid thickness. At large λ, the liquid thickness becomes very large and the profile is smooth over the entire horizontal range. Gradually the gap approaches the shape of the rigid cylinder, since the elastic deformation of the coating becomes small in comparison to the lubricant thickness.

Numerical results
We now take a closer look at the deformations over the full range of λ for the longwave description and the finite element simulations -both pertaining to thin compressible coatings. In figure 4 we therefore plot the fluid layer thickness at the centerh 0 = h 0 /u 0 , as a function of the dimensionless velocity λ. The layer thickness is furthermore compensated by the scaling λ 1/2 that will be derived analytically below. From the figure we can see that for values of λ 10 −4 the two different numerical methods agree with one another -confirming the validity of the long-wave expansion at moderate and large values of λ. The long-wave data (represented as stars), exhibit two asymptotic regimes: h 0 ∼ λ 1/2 (dashed line) and h 0 ∼ λ 4/7 (dash-dotted line). However, it is also clear that at the smallest λ, the long-wave description breaks down, since the numerical solution of the true elastohydrodynamic problem (represented by triangles) gives yet another scaling law h 0 ∼ λ 3/5 (solid line). We will see below that this limit emerges when the width of the boundary layer becomes comparable to the thickness of the coating d, so that the long-wave description is no longer justified.
In summary, we thus find the emergence of 3 regimes at small, intermediate and large values of λ. The analysis for large λ has already been performed by Skotheim x/a h/u0 λ = 1.8 × 10 1 λ = 6.9 × 10 0 λ = 1.0 × 10 0 λ = 5.9 × 10 −2 Dry solution Figure 3. Numerical film thickness profilesh = h/u0 as a function ofx = x/a, for the case of a compressible thin elastic layer. Lubricated contacts correspond to λ > 0, and are obtained from numerical integration of (2.12). The case λ = 0 is a "dry" contact, given by the analytical expression of (3.1b).  & Mahadevan (2004,2005). Hence, below we derive the intermediate asymptotics (using the long-wave description at small λ), and the ultimate asymptotics λ 1 based on the short-wave description.

Long wavelength elasticity: Dry solution
Before turning to the dynamical case, we first provide the dry (static) solution of the long-wave approximation, given by (2.10). Formally, this corresponds to λ = 0, for which no flow takes place. This dry solution therefore corresponds to the situation where the pressure vanishes outside the contact region,p = 0 for |x| > 1, and the gap thickness vanishes within the contact,h = 0 for |x| < 1. The corresponding dry solution to (2.10) reads,p where Θ is the Heaviside step function. Expanding these near the edges of the contact, located atx = ∓1, we find This linear behaviour of the gap can be observed in figure 3, at the edge of the dry solution. This asymptotics is of particular importance to the dynamical case at small speed, where λ is small but nonzero. In that limit, lubrication effects quickly decay outside the narrow boundary layers near the edge of the contact, so that the static solution is approached. Specifically, the boundary layer solutions must therefore be matched to the linear asymptotics given by (3.2).

Long wavelength elasticity: Lubricated solution
Once the cylinder starts to move, a thin film of liquid is entrained within the contact region and in particular, the left-right symmetry of the deformation is broken (fig. 3). Liquid is squeezed into the contact at the inlet, atx = −1, and is pushed out at the outlet, atx = 1. To analyze the behavior of (2.12) at the two lubricated edges, we look for similarity solutions of the form Here the + sign refers to the inlet and − to the outlet. In terms of the similarity variables, the matching condition (3.2) becomes H(ξ) ∓2ξ for ξ → ∓∞. This λ-independent farfield condition requires that α = β. The Reynolds lubrication equation (2.12) written in terms of the above similarity variables, gives Anticipating that β > 0, i.e. thin boundary layers for λ 1, the term ∼ λ β is subdominant with respect to the constant. For the remaining equation to be independent of λ, we need α − β = 1 − 2α. Combined with the matching condition of α = β, we thus find the exponents The solution to this equation must match to H(ξ → ∓∞) = ∓2ξ, as discussed before. Furthermore, the gap should converge to a film of constant thickness inside the contact, leading to ∂H ∂ξ (ξ → ±∞) = 0.

Inlet region
At the inlet (3.5) becomes (3.6) We first wish to find the constant thickness H 0 that is approached when entering the central zone of the contact (corresponding to ξ → ∞). Interestingly, H 0 is not simply equal to H * , but is to be determined from the equation H 0 − H * − 2H 3 0 = 0. For this film thickness to take on a positive real value, it is required that H * 1 9 √ 6. However, under this condition there are still two possible solution branches of H 0 for a given H * .
For the selection of the relevant H 0 , we now investigate whether the solution indeed approaches H 0 for ξ → ∞. We therefore linearise (3.6) using H = H 0 + H 1 (ξ), which gives This shows that the solution exponentially approaches a constant film thickness when H 0 > 1 6 √ 6. This condition selects the branch of acceptable H 0 , which for a given value of H * indeed decays to a flat film below the contact. We need to separately consider the special case where H 0 = 1 6 √ 6 and H * = 1 9 √ 6. In this case the similarity solution can even be found in closed form, This solution exhibits an algebraic decay to the flat film of thickness H 0 as ξ → ∞, and therefore should also be considered. It will turn out that, indeed, the special case of H 0 = 1 6 √ 6 and H * = 1 9 √ 6 provides the relevant solution for the lubrication problem. This will be shown by analysing the "central region", as done below in Sec. 3.3.3. Hence, the correct inlet thickness reads (3.9) while the corresponding value for H * is given by (3.10)

Outlet region
The similarity equation for the outlet differs from (3.6) only in the sign of the last term on the right hand side, (3.11) Anticipating that H * = 1 9 √ 6, we could find a lengthy closed form solution to (3.11) using Mathematica (not given here). Remarkably, the solution approaches a constant thickness H out as ξ → −∞ that turns out to be different from the inlet thickness H in . Indeed, solving ∂H ∂ξ = 0 for H * = 1 9 √ 6, one finds (3.12) which is significantly below H in . This means the liquid film cannot possibly be flat within the central region connecting the inlet and outlet, and we need to treat this region separately.

Central region
As just demonstrated the constant film thickness at the inlet and outlet cannot be equal, these satisfy two different equations for the same value of H * . To find the variation in film thickness across the central zone, another similarity solution has to be introduced, namely,h = λ 1/2 H(x). (3.13) Since this similarity solution does not scale the solution in horizontal direction, this equation is valid when |x| 1 − λ 1/2 , i.e. over the full width of the contact excluding the boundary layers. Substituting this into the lubrication equation (2.12) gives an algebraic equation for the film thickness, (3.14) Since H describes the central zone up to the boundary layers, (3.14) needs to cover the entire domain −1 <x < 1. One verifies that this is only the case when H * 1 9 √ 6. This condition is exactly opposite to the requirement derived from the inlet solution, which only gives physical solutions for H * 1 9 √ 6. Hence, the only possibility for the central zone to connect inlet and outlet is when H * = 1 9 √ 6, as anticipated above. Figure 5 summarizes the key results of this section. The panels (a-c) respectively report the similarity solutions for the narrow boundary layers at (a) the inlet, (c) the outlet and (b) the central region that connects the two. The symbols represent the analytical solutions derived above, and these are plotted along with scaled FEM numerical solutions (solid lines, for various λ). An excellent agreement is observed in all three regions.

Summary and numerical validation
The similarity solution also provides the exact expression for the fluid film thickness h 0 . For this we evaluate (3.14) at x = 0, which gives immediately that the (scaled) film thickness at the center of the contact is equal H * = 1 9 √ 6. In dimensional form it gives This result, including the theoretical prefactor, is used in fig. 4 (dashed line), which shows that the numerical solutions of the long-wavelength approximation indeed converge towards the derived similarity solution. With respect to the full resolution of the elastohydrodynamic problem using finite elements, the solution (3.15) appears as an intermediate (long-wave) asymptotics.
3.4. Short wave elasticity: the ultimate asymptotics for λ 1 The numerical results obtained by finite element simulations, shown in fig. 4, revealed that at very small λ, another limit develops. The reason for this is that the boundary layer width in the long-wave theory was found to scale as ∼ λ 1/2 , so that at sufficiently small λ it becomes comparable to the thickness of the elastic layer -at which point the longwavelength description becomes inconsistent. Hence, the ultimate small λ asymptotics involves the hierarchy of scales d a, and requires the full kernel of (2.5). Below we first consider the corresponding static problem, to identify the appropriate horizontal, vertical and pressure scales. Subsequently, in the next subsection, we use the dry solution to match the boundary layers.

Dry solution
In fact, even the static indentation problem with d a, already involves the full kernel (2.5). This problem goes back to Meijers (1968), who considered the Hertz contact on layers of arbitrary thickness. On the scale of the contact width a, the static solution of the long-wavelength description, (3.1a, 3.1b), is perfectly valid when d a. Near the edges, however, the linear behaviour of the pressure gives way to a square-root dependence, as is illustrated in fig. 6(a). The reason for this is that near the edges, the elasticity is governed by the short-wavelength limit of the the elastic kernel (2.5). Introducing ζ = x ± a, indicating the unscaled distance to the contact edge, this gives the elastic relation which combines (2.1) and the deformation from table 1 in the limit qd 1. It is wellknown that this equation gives (3.17) at the contact edge (Snoeijer et al. 2013), where the scaling properties of A remain to be determined. This behaviour is indeed manifestly different from the linear behaviour obtained from the long-wave expansion (3.2), which in terms of a dimensional scaling law can be written as with a given by (2.9). To estimate the value of the constant A, we now use a matching of the two descriptions, as is depicted schematically in fig. 6(a). This is done by equating (3.17) to (3.18) at a distance |ζ| ∼ d. This matching gives an expression for the parameter A, where A p is a (dimensionless) constant. The value of the constant is obtained from a numerical solution using the finite element method, of the static dry contact problem. The result is given in fig. 6(b), from which we deduce that A p ≈ 1.33. Now that we know the behavior of the static solution near the contact edge, we can identify the proper scalings for p, h and ζ. This is important, since these will subsequently be used to describe the dynamic boundary layers. To establish the scalings, we first take two derivatives of (3.16), which gives where we also performed an integration by parts. Demanding that the three terms in this equation are of the same order, while also respecting the scaling of A in (3.17,3.19), we find the appropriate scalings to be With this, the dimensionless equation describing the elasticity of the layer reads ∂ 2ĥ (3.22) while the asymptotics of the static solution reduces tô p ±2ζ. (3.23)

Lubricated solution
We now turn to the dynamical case, by expressing the Reynolds equation (2.3) in the new variables of (3.21). This gives which is of the same form as before. However, owing to the new scaling, the parameter that gives the dimensionless velocity is different, and reads The analysis of the boundary layer is now given by the system (3.22, 3.24), with matching condition (3.23). This reduced problem is strictly identical to that of the boundary layers on a half-space, as previously described by Snoeijer et al. (2013). Hence the solution is exactly the same, and gives a thicknessĥ ∼ λ 3/5 s , a boundary layer widtĥ ζ ∼ λ 2/5 s , while the pressurep ∼ λ 1/5 s . For detailed derivations we refer to Snoeijer et al. (2013). Note, however, that in dimensional variables the result is different. In particular, we find for the dimensional flowing thickness where the numerical constant 0.4467 was determined in Snoeijer et al. (2013) by solvinĝ h * from the boundary layer problem, while we remind that A p ≈ 1.33. This result is given as the solid line in fig. 4, and indeed accurately describes the ultimate λ 1 asymptotics.
The structure of the boundary layers is further illustrated in fig. 7. There, the similarity solutions from Snoeijer et al. (2013) are reproduced (symbols), and compared to the scaled FEM simulations at very small λ. All curves indeed collapse, confirming the validity of the similarity analysis. Interestingly, for this case the thickness of the lubrication layers leaving the inlet and entering the outlet directly align. The central zone is therefore perfectly flat -this is to be contrasted with the result for the long-wave similarity solutions of fig. 5, for which the inlet and outlet were connected by a nontrivial central region.

Discussion
The aim of this work is to present all scaling laws governing the movement of a cylinder over an elastic substrate in the presence of a fluid. Some of these scaling laws were derived previously in the literature, but focussing on separate regimes -importantly, the limiting cases of thin elastic layers with large deformations had not previously been considered. This important gap is filled in the present paper, thereby offering a complete and comprehensive overview. Table 2 summarises both previous and our newly obtained results. All the different regimes of soft lubrication are organised in terms of the length scales d (thickness of the elastic coating), a (size of the contact), u 0 (scale of elastic deformation) and (boundary layer width near the contact edges). In terms of experimental parameters, the table should be read as follows. For fixed material properties, magnitude of the applied load L determines the typical size of the contact a. Hence, upon specifying the load one selects one of the rows of table 2, comparing the contact size to the elastic layer thickness. Then, at very small velocity V one starts in the left column, where the flowing thickness h 0 u 0 , referring to the case of relatively large elastic deformation. Subsequently, upon increasing the velocity, the boundary layer and the lubrication film h 0 increases, and one progressively moves to the second and third columns, ultimately reaching the small deformation limit (h 0 u 0 ). There are a few important features in this table. At very small velocities (left column), the scaling h 0 ∼ V 3/5 is the same for all cases. In this case, the width of the boundary layer is the smallest scale in the system and thereby these exhibit a universal structure. Differences arise only from matching these boundary layers to different (static) outer solutions. On the contrary, at very large velocity (right column) the scaling laws always involve the combination (ηV ) 2 /L. So, when working at a fixed separation distance instead of at fixed load, one finds the well known scaling law for the lift force L ∼ (ηV ) 2 (Sekimoto & Leibler 1993;Skotheim & Mahadevan 2004;Greenwood 2020). A key result of our analysis is the emergence of an intermediate asymptotics, h 0 ∼ V 1/2 for the thin, compressible layer -closed form analytical solutions for this case are provided.
Not physically possible ( a) 3π 8 no long-wave asymptotics Another new result of our study is the scaling law for a thin incompressible layer subjected to large indentation. Leaving the details in appendix A, here we discuss the key features of this limit. In the limit where the boundary layer represents the smallest scale, d, the only change compared to the compressible case is in the static "dry" solution, leading to different exponents of geometric and material parameters. However, the scaling h 0 ∼ V 3/5 is robust. Next we consider the intermediate range, d a, where thesmall boundary layers are in between the coating thickness and the contact size. Interestingly, this case does not admit a consistent long-wave description, and therefore a proper intermediate asymptotic regime appears to be missing. The resulting longwavelength ODE can be derived -unlike the compressible case, however, it does not admit solutions at small λ that connect the inlet and the outlet region. For completeness, we show in fig. 8 the numerical relation between the central thickness h 0 and the velocity V obtained by finite element simulations. It very nicely confirms the scalings of the bottom row in table 2, corresponding to the small and large speed asymptotics for thin incompressible layers.
In conclusion, our study identifies key geometric parameters in the soft lubrication problem of compressible and incompressible elastic solids. From an experimental perspective, our results point out the importance of: (i) the ratio of elastic coating thickness to contact size, (ii) the Poisson ratio of the layer (compressible versus incompressible), (iii) the ratio of elastic deformation over the lubricant thickness. When interpreting experiments, one must verify the relevant regime, with Table 2 serving as a guide. Our unified framework also connects the two most studied limits of the problem: the lubricated Hertzian contact regime, relevant for many industrial applications, and the large separation (small deformation) regime, important in microfluidic and biological contexts. We anticipate that the results presented here to have important implications beyond the realm of soft lubrication, e.g. in understanding the rheology of suspension flow in channels (Rosti et al. 2019), size dependent migration of particles (Rallabandi et al. Distinct asymptotic limits are observed at low velocity (solid line,h0 ∼ λ 3/5 ) and high velocity (dashed line,h0 ∼ λ 4/9 ). The prefactor of the high velocity asymptote is derived analytically, while for the low velocity asymptote the data was fitted by (A 11), with Ai,p 1.85.

2018)
, and direct measurement of interfacial rheology in soft matter systems (Garcia et al. 2016).
corresponding scales need to be identified from matching to the static solution (now, on thin incompressible layers). In particular, we know the asymptotics to be of the form p A i ±ζ, (A 2) where ζ = x ± χa i is the distance from the edge, and we need to identify A i of the static solution.
For this, we proceed as in Section 3.4. We first consider the static long-wave solution, from (A 1), which gives the relevant horizontal length scale a i , a i = 5 6 √ 3 We introduce the scalings, As before the dry problem can be split into two distinct regions, which gives the conditionsh = 0 where |x| < χ andp = 0 where |x| > χ. Note that the size of the contact is not given by a i , but by χa i . This is because incompressibility causes the contact to be larger than the length over which the cylinder is below the surface of the undeformed layer. Incompressibility dictates that the volume of the layer is conserved,thus ∞ −∞ h − (x 2 − 1) dx = 0, which gives that the size of the contact χ = √ 3. Using this result the dry solution is found to be, h = x 2 − 1 Θ x 2 − χ 2 (A 5a) p = 1 9 9 − 6x 2 +x 4 Θ χ 2 −x 2 . (A 5b) Interestingly, while the pressure is a continuous function, the film thickness is a discontinuous function. Equation (A 1) predicts that the surface of the layer follows the cylinder when |x| < χ, and then suddenly becomes flat. We now need to compare the "true" asymptotics of the statics solution, given by (A 2), to that of the long-wave approximation, which near the edge of the contact (x = ∓χ) is given byp which in dimensional form gives the scaling, Equating this expression to (A 2) at a distance |ζ| ∼ d, we find where A i,p is a numerical constant that remains to be determined. To recover the "universal" boundary layer equations, we scale the problem as,