2.1. Review: inviscid Burgers equation

In this subsection, we review the inviscid Burgers equation and the corresponding similarity solution associated with shock formation. The approach of Eggers & Fontelos (Reference Eggers and Fontelos2009) is followed closely. The inviscid Burgers equation for velocity $u=u(x,t)$ is given by

(2.2) \begin{equation} \frac {\partial u}{\partial t} + u\frac {\partial u}{\partial x}=0, \end{equation}

where $x$ denotes the spatial coordinate and $t$ denotes time. Consider the formation of a shock at time $t_*$ , occurring at $x=x_*$ with $u_*=u(x_*,t_*)$ . The mechanism behind the shock formation can be described as follows. Close enough to the singularity in space and time, (2.2) is to leading order just a constant speed advection equation $(\partial u/\partial t)+u_*(\partial u/\partial x)=0$ , which on its own would not lead to shock formation. Instead, the next-order nonlinearity $(u-u_*)(\partial u/\partial x)$ gives rise to the shock formation.

It is convenient when there is formation of a front to perform analysis in the frame of the translating front. Then, in order to investigate the behaviour near the shock in space and time, the variables $u':=u-u_*$ , $\overline {x}:=x-x_*+u_*\tau$ and $\tau :=t_*-t$ are considered. The definition of $\tau$ is chosen such that as $t\rightarrow t_*^-$ , then $ \,\tau \rightarrow 0^+$ . Substitution of the definitions into (2.2) gives

(2.3) \begin{equation} \frac {\partial u'}{\partial \tau }-u'\frac {\partial u'}{\partial \overline {x}}=0. \end{equation}

A local similarity solution of the shock as $\tau \rightarrow 0^+$ can be found by considering the self-similarity ansatz (Eggers & Fontelos Reference Eggers and Fontelos2009)

(2.4) \begin{equation} u'=\tau ^{\alpha -1}F(\xi ), \end{equation}

where $\xi :=\overline {x}\tau ^{-\alpha }$ and the scaling comes from the balance between $\partial u'/\partial \tau$ and $u'(\partial u'/\partial \overline {x})$ . Note that the constant $\alpha \geqslant 1$ (since $u'$ does not blow up) cannot be deduced from scaling analysis alone and is deduced from stability considerations. In other words, the self-similarity is of the second kind. Then, (2.3) and (2.4) give

(2.5) \begin{equation} (\alpha -1)F-\alpha \xi \frac {\textrm{d}F}{\textrm{d}\xi }-F\frac {\textrm{d}F}{\textrm{d}\xi }=0, \end{equation}

which can be integrated when $\alpha -1 \neq 0$ to give

(2.6) \begin{equation} \xi = -F - K F^{\frac {\alpha }{\alpha -1}} ,\end{equation}

for some constant $K\gt 0$ ; the case $\alpha = 1$ gives $F=-\xi$ , which cannot be matched to the solution away from the singularity (Eggers & Fontelos Reference Eggers and Fontelos2009). There is then a discrete number of possible $\alpha$ given by the constraint of $F$ being smooth everywhere ( $\alpha /(\alpha -1)$ needs to be an odd integer). Out of the possible $\alpha$ , it can be shown that only $\alpha = {3}/{2}$ gives a stable similarity solution (Eggers & Fontelos Reference Eggers and Fontelos2009). In summary, unravelling definitions, the similarity solution for the Burgers equation, (2.2), is given by

(2.7) \begin{equation} u=u_*+(t_*-t)^{\frac {1}{2}}F\left (\frac {x-x_*+u_*(t_*-t)}{(t_*-t)^{\frac {3}{2}}}\right )\!, \end{equation}

where $\xi = -F - K F^3$ for some constant $K\gt 0$ .

2.2. Derivation of the similarity solution in the IM regime

In this subsection, we give the derivation of the similarity solution associated with shock formation in the IM regime. Throughout the derivation, analogies to the Burgers equation are given (see § 2.1). Consider the formation of a shock at time $t_*$ , occurring at $r = r_*$ with

(2.8) \begin{equation} (u_*,h_*,\varGamma _*):=(u(r_*,t_*),h(r_*,t_*), \varGamma (r_*,t_*)). \end{equation}

The singularity values $(u_*,h_*,\varGamma _*,r_*,t_*)$ are obtained numerically (see supplementary material) and are taken to be known. In order to investigate the behaviour near the shock, let

(2.9) \begin{equation} (u',h',\varGamma ',r', \tau ):=(u-u_*,h-h_*,\varGamma -\varGamma _*, r-r_*, t_*-t). \end{equation}

Then, (2.1), (1.3b ) and (1.3c ) give

(2.10a) \begin{align} \frac {\partial u'}{\partial \tau }- (u_*+u')\frac {\partial u'}{\partial r'} +\frac {2 \frac {\textrm{d} \sigma }{\textrm{d} \varGamma }(\varGamma _*+\varGamma ')}{h_*+h'}\frac {\partial \varGamma '}{\partial r'}&= 0 , \end{align}

(2.10b) \begin{align} \frac {\partial h'}{\partial \tau }-(u_*+u') \frac {\partial h'}{\partial r'}-\frac {(h_*+h')(u_*+u')}{r_*+r'}-(h_*+h')\frac {\partial u'}{\partial r'}&= 0, \end{align}

(2.10c) \begin{align} \frac {\partial \varGamma '}{\partial \tau }-(u_*+u') \frac {\partial \varGamma '}{\partial r'}-\frac {(\varGamma _*+\varGamma ')(u_*+u')}{r_*+r'}-(\varGamma _*+\varGamma ')\frac {\partial u'}{\partial r'}&= 0. \end{align}

As $\tau \rightarrow 0^+$ , the shock forms and hence the derivative terms become singular (e.g. $|\partial u'/\partial r'|\rightarrow \infty$ ). Then, finite terms such as $h_*u_*/r_*$ are not leading-order terms. To leading order, as $\tau \rightarrow 0^+$ , (2.10) simplifies to

(2.11) \begin{equation} \frac {\partial }{\partial \tau }\left [ \begin{array}{c} u' \\ h'\\ \varGamma '\\ \end{array} \right ]= \left [ \begin{array}{ccc} u_* & 0 & \frac {V^2}{\varGamma _*} \\[4pt] h_* & u_* & 0 \\ \varGamma _* & 0 & u_*\\ \end{array} \right ] \frac {\partial }{\partial r'}\left [ \begin{array}{c} u' \\ h'\\ \varGamma '\\ \end{array} \right ]\!, \end{equation}

where a velocity $V$ is defined

(2.12) \begin{equation} V:=\left (-\frac {2\frac {\textrm{d} \sigma }{\textrm{d}\varGamma }(\varGamma _*)\varGamma _*}{h_*}\right )^{\frac {1}{2}}. \end{equation}

In this work, we only consider $(\textrm{d}\sigma /\textrm{d} \varGamma )|_{\varGamma = \varGamma _*}\lt 0$ , as is physically relevant for surfactants, and hence $V\gt 0$ . The leading-order Burgers equation (2.2) is an advection equation and consequently a change of coordinates is applied. Similarly, the leading-order matrix system (2.11) is simply a linear advection equation, where the matrix has eigenvalues $u_*,u_*+V,u_*-V$ with corresponding eigenvectors $(0,1,0), (V,h_*,\varGamma _*),(V,-h_*,-\varGamma _*)$ . Thus, in order to consider local shock formation, we change coordinates to be along a particular characteristic. The rest of this section considers shock formation along the $u_*+V$ characteristic, with space–time coordinates ( $\overline {r}, \tau$ ), where $\overline {r}:=r'+(u_*+V)\tau$ . The $u_*$ characteristic has an eigenvector that has a zero $\varGamma$ component and is not considered. The similarity solution for the shock formation along the $u_*-V$ characteristic can be derived with an analogous method.

Now, we keep account of the order of magnitude of the terms appearing in (2.10) with respect to $\tau$ as $\tau \rightarrow 0^+$ . As mentioned, to leading order, (2.10) is a linear advection equation and hence $u',h',\varGamma ' =\textit {O}(\tau ^{\beta })$ for some constant $\beta$ . Consequently, like the Burgers equation, the nonlinearity gives rise to the singularity and hence there is a balance $(\partial u'/\partial \tau )\sim u'(\partial u'/\partial \overline {r})$ . Then, letting $\overline {r} =\textit {O}(\tau ^{\alpha })$ for the characteristic width of the shock region, $\beta = \alpha - 1$ . Thus, from the leading-order expansion (2.11), since the $u_*+V$ characteristic has eigenvector $(V,h_*,\varGamma _*)$ , we may write

(2.13a) \begin{align} u'&=V(A(\overline {r},\tau )+q_1 \tau ^{\alpha -1}+f_1(\overline {r},\tau )) ,\end{align}

(2.13b) \begin{align} h'&=h_*(A(\overline {r},\tau )+q_2 \tau ^{\alpha -1}+f_2(\overline {r},\tau )) ,\end{align}

(2.13c) \begin{align} \varGamma '&=\varGamma _*(A(\overline {r},\tau )+q_3 \tau ^{\alpha -1}+f_3(\overline {r},\tau )), \end{align}

where $A=\textit {O}(\tau ^{\alpha -1})$ and $f_1,f_2,f_3$ are correction terms much smaller than $\textit {O}(\tau ^{\alpha -1})$ and $q_1,q_2,q_3$ are constants. Since $A$ is a function to be found, without loss of generality, $q_3 = -q_1$ (the explicit change of variables would be $(\tilde {A}, \tilde {q}_i) = (A+((q_1+q_3)/2)\tau^{\alpha-1}, q_i - (q_1+q_3)/2)$ for $i = 1,2,3$ ). Then, substituting (2.13), correct to $\textit {O}(\tau ^{\alpha -2})$ , into (2.10) gives (see Appendix A)

(2.14a) \begin{align} V\frac {\partial (A+q_1\tau ^{\alpha -1})}{\partial \tau }+V^2 \left (q_2-q_1-\frac {2\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}q_1\right )\tau ^{\alpha -1}\frac {\partial A}{\partial \overline {r}}&& \nonumber \\ +\frac {2\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*}A\frac {\partial A}{\partial \overline {r}} + V^2 \left (\frac {\partial f_1}{\partial \overline {r}}- \frac {\partial f_3}{\partial \overline {r}}\right )&= 0 , \end{align}

(2.14b) \begin{align} h_*\frac {\partial (A+q_2 \tau ^{\alpha -1})}{\partial \tau }-Vh_*(q_1+q_2)\tau ^{\alpha -1}\frac {\partial A}{\partial \overline {r}}-2h_*VA\frac {\partial A}{\partial \overline {r}}&&\nonumber \\ +Vh_*\left (\frac {\partial f_2}{\partial \overline {r}} - \frac {\partial f_1}{\partial \overline {r}}\right )&= 0, \end{align}

(2.14c) \begin{align} \varGamma _*\frac {\partial (A-q_1 \tau ^{\alpha -1})}{\partial \tau }-2\varGamma _*VA\frac {\partial A}{\partial \overline {r}}+V\varGamma _*\left (\frac {\partial f_3}{\partial \overline {r}} - \frac {\partial f_1}{\partial \overline {r}}\right )&= 0. \\[9pt] \nonumber \end{align}

Rearranging and combining (2.14a ) and (2.14c ) gives

(2.15) \begin{equation} \frac {\partial A}{\partial \tau }+\frac {V}{2}\left (q_2-q_1-\frac {2\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}q_1\right )\tau ^{\alpha -1}\frac {\partial A}{\partial \overline {r}}-V\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_* V^2}\right )A\frac {\partial A}{\partial \overline {r}}=0, \end{equation}

which upon a change of variables

(2.16) \begin{equation} \overline {r}_{s}:=\frac {\overline {r}-\frac {V}{2}\left (q_2-q_1-\frac {2\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}q_1\right )\alpha ^{-1}\tau ^{\alpha }}{V\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_* V^2}\right )} ,\end{equation}

gives

(2.17) \begin{equation} \frac {\partial A}{\partial \tau }-A\frac {\partial A}{\partial \overline {r}_{s}}=0. \end{equation}

Therefore, we have arrived at precisely the inviscid Burgers (2.3). Physically, $\overline {r}_{s}$ is the spatial coordinate in the frame of the moving shock. As per § 2.1, $\alpha = {3}/{2}$ and $A$ has a similarity solution given by $A = \tau ^{(1/2)}F(\eta )$ , where $\eta :=\overline {r}_s\tau ^{-(3/2)}$ and

(2.18) \begin{equation} \eta = -F-KF^3 ,\end{equation}

for some constant $K\gt 0$ . Furthermore, by ensuring that the similarity solution may be matched to the outer region away from the local singularity region, it is possible to deduce that $q_1 = q_2 = 0$ (see Appendix B).

Then, unravelling all the definitions, the similarity solutions for the finite-time singularity for shock formation along the $u_*+V$ characteristic (see Appendix C for $u_*-V$ case) are given by

(2.19a) \begin{align} u &= u_* + V (t_*-t)^{\frac {1}{2}}F\left (\frac {r-r_*+(u_*+V)\left (t_*-t\right ) }{\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}\right )V(t_*-t)^{\frac {3}{2}}}\right )\!, \end{align}

(2.19b) \begin{align} h &= h_*+h_*(t_*-t)^{\frac {1}{2}}F\left (\frac {r-r_*+(u_*+V)\left (t_*-t\right )}{\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}\right )V(t_*-t)^{\frac {3}{2}}}\right )\!, \end{align}

(2.19c) \begin{align} \varGamma &= \varGamma _* + \varGamma _*(t_*-t)^{\frac {1}{2}}F\left (\frac {r-r_*+(u_*+V)\left (t_*-t\right )}{\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}\right )V(t_*-t)^{\frac {3}{2}}}\right )\!, \end{align}

where $F$ is given by (2.18). Corrections are $\textit {O}(t_*-t)$ as the shock forms with $t \rightarrow t_*^{-}$ . In summary, like the Burgers equation, the shock region has width $\textit {O}((t_*-t)^{{3}/{2}})$ and $u-u_*,h-h_*,\varGamma -\varGamma _* = \textit {O}((t_*-t)^{{1}/{2}})$ in the shock as $t \rightarrow t_*^{-}$ .

The method of translating coordinates along eigenvectors and reducing to Burgers equation, as shown above, is generalisable to a wider class of hyperbolic coupled partial differential equations and will be the subject of a future manuscript.

The details of shock formation may now be discussed. From (2.18), we have that $\textrm{d}F/\textrm{d}\eta =-(1+3KF^2)^{-1}$ and hence

(2.20) \begin{equation} \max \left |\frac {\textrm{d}F}{\textrm{d}\eta }\right |=1, \end{equation}

which in turn gives that

(2.21) \begin{equation} \max \left |\frac {\partial u}{\partial r}\right |=\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}\right )^{-1}(t_*-t)^{-1} \text{ as } t \rightarrow t_*^{-}, \end{equation}

where such an expression is useful for verification (see § 2.3) since (2.21) does not contain any undetermined coefficients (recall that the singularity values such as $\varGamma _*$ are known numerically). Additionally, (2.21) shows that the nonlinearity of the surface tension isotherm has an effect on the shock formation. More precisely, if the surface tension isotherm is locally concave $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)|_{\varGamma = \varGamma _*}\lt 0$ , then the shock is weakened, and if the surface tension isotherm is locally convex $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)|_{\varGamma = \varGamma _*}\gt 0$ , the shock is strengthened.

Another useful expression, to be used in § 2.3 below for the purpose of additional numerical verification, is that

(2.22) \begin{equation} \max \left |\frac {\partial ^2 u}{\partial r^2}\right |=\frac {25 \sqrt {15}}{108}K^{\frac {1}{2}}V^{-1}\left (1-\frac {\frac {\textrm{d}^2 \sigma }{\textrm{d} \varGamma ^2}(\varGamma _*)\varGamma _*^2}{h_*V^2}\right )^{-2}(t_*-t)^{-\frac {5}{2}}\text{ as } t \rightarrow t_*^{-}, \end{equation}

since $\max |\textrm{d}^2F/\textrm{d}\eta ^2|=(25 \sqrt {15}/108)K^{{1}/{2}}$ from (2.18). Equation (2.22) is also used to obtain an estimate for the constant $K$ .

2.3. Verification of the similarity solution

In this subsection, we verify the similarity solution given by (2.19). In particular, a linear isotherm $\sigma (\varGamma ) = - \varGamma$ with initial surfactant distribution $\varGamma _i(r) = e^{-r^2}$ is considered.

Since it is the derivatives $(\partial u/\partial r,\partial h/\partial r,\partial \varGamma /\partial r)$ that become singular, rather than the variables $(u,h,\varGamma )$ themselves, the collapse to the similarity solution is harder to deduce systematically from seeing the profiles of $u,h$ and $\varGamma$ alone. Instead, we may verify the similarity solution using log–log plots to find the exponents. Since we wish to show two power-law behaviours, $u-u_* \sim (t_*-t)^{{1}/{2}}$ and $\overline {r}=r-r_*+(u_*+V)(t_*-t)\sim (t_*-t)^{{3}/{2}}$ , it is sufficient to check the two predictions (2.21) and (2.22). In figure 4(a) we plot $\max |\partial u/\partial r|$ versus $t_*-t$ and in figure 4(b) $\max |\partial ^2 u/\partial r^2|$ versus $t_*-t$ is reported. Numerical solutions of the thin-film equations are shown by solid curves and analytical similarity predictions (2.21) and (2.22) are shown by dashed lines. There are no fitting parameters in figure 4(a). In figure 4(b), $K \approx 0.083$ is chosen as the best fit, which gives an estimate for $K$ for the example shown. Figure 4 shows that, indeed, $(\partial u/\partial r)\sim (t_*-t)^{-1}$ and $(\partial ^2 u/\partial r^2)\sim (t_*-t)^{- (5/2)}$ . The agreement breaks down at around $t_*-t\approx 5 \times 10^{-4}$ due to numerical inaccuracies, arising from issues such as not knowing $t_*$ analytically (see supplementary material). This is a standard feature of log–log plots verifying numerical similarity solutions (Eggers & Fontelos Reference Eggers and Fontelos2015).

Figure 4. Systematic verification of the similarity solution for an example with a linear isotherm $\sigma (\varGamma )=-\varGamma$ and initial surfactant distribution $\varGamma _i(r) = e^{-r^2}$ . The solid curves are the numerical solutions of the IM thin-film equations and the dashed lines are the similarity solution predictions (2.21, 2.22). (a) Log–log plot of $\max |\partial u/\partial r|$ (no fitting parameters). (b) Log–log plot of $\max |\partial ^2 u/\partial r^2|$ , where $K \approx 0.083$ is chosen as the best fit for the similarity solution prediction (2.22).

The collapse of the profiles for $u, h$ and $\varGamma$ is shown in figure 5, where the solid curves plotted have colours corresponding to $-\log (t_*-t)$ . Figure 5(a) shows the horizontal velocity $u$ evolving over time, where the spatial coordinate is given by $\overline {r}=r-r_*+(u_*+V)(t_*-t)$ . Figure 5(b) shows the rescaled horizonal velocity $(u-u_*)V^{-1}(t_*-t)^{-(1/2)}$ versus the rescaled spatial coordinate $\overline {r}V^{-1}(t_*-t)^{- ({3}/{2})}$ , which is the expected similarity coordinate as predicted by (2.19a ) . Similarly, figure 5(c,e) shows the variables $h$ and $\varGamma$ and figure 5(d, f) displays the rescaled thickness and surfactant concentration versus the rescaled spatial coordinate. The dashed curve in figure 5(b, d, f) is $x = -y - Ky^3$ (where the horizontal axis is $x$ and vertical axis is $y$ ). The coefficient $K\approx 0.083$ is that estimated from $\max |\partial ^2 u/\partial r^2|$ in figure 4(b), as above. We see in figure 5(b,d, f), respectively, that the dependent variables $u, h$ and $\varGamma$ collapse as $t \rightarrow t_*^-$ , in good agreement with (2.19). With a more accurate numerical procedure, one could take smaller $t_*-t$ and also know $u_*, h_*, \varGamma _*,r_*, t_*$ more accurately (see supplementary material).

Figure 5. Shock formation with isotherm $\sigma (\varGamma )=-\varGamma$ and initial surfactant distribution $\varGamma _i(r) = e^{-r^2}$ . A finite-time shock forms at $t=t_*$ ( $\approx 1.136$ ). Colour bar: the solid curves have colour according to $-\log (t_*-t)$ . (a) The horizontal velocity $u$ with the spatial coordinate given by $\overline {r} = r-r_*+(u_*+V)(t_*-t)$ , which is the coordinate system moving with the inflection point. (b) The appropriately rescaled horizontal velocity versus the appropriately rescaled spatial coordinate (i.e. the similarity solution). (c,d,e, f) Analogous to (a–b) but for thickness $h$ and surfactant concentration $\varGamma$ . The dashed curves in (b), (d) and ( f) are the curve $x = -y-Ky^3$ where the horizontal axis is $x$ , the vertical axis is $y$ and $K \approx 0.083$ is a constant numerically estimated using relation (2.22) and figure 4(b).

2.4. Nonlinear effects of the surface tension isotherm

Here, some details of two examples with nonlinear effects of the surface tension isotherm are discussed. The first example is the Langmuir isotherm, $\sigma _{\textit{conca}v\textit{e}}(\varGamma )=\varGamma _{\infty }\log (1-(\varGamma /\varGamma _{\infty }))$ with $\varGamma _{\infty }=1.1$ , corresponding to a concave isotherm. The prefactor of $\sigma _{\textit{conca}v\textit{e}}(\varGamma )$ is chosen so that $(\textrm{d} \sigma _{\textit{conca}v\textit{e}}/\textrm{d} \varGamma )|_{\varGamma = 0}=-1$ , which is how the problem was non-dimensionalised (see 1.2) and the choice of $\varGamma _{\infty }=1.1$ is so that the effects of concavity can be seen clearly. The second example is convex, $\sigma _{\textit{con}v\textit{ex}}(\varGamma )=-\varGamma _0 \tanh (\varGamma /\varGamma _0)$ with $\varGamma _0 = 0.1$ . Again, the prefactor is chosen so that $(\textrm{d} \sigma _{\textit{con}v\textit{ex}}/\textrm{d} \varGamma )|_{\varGamma = 0}=-1$ and the choice $\varGamma _0 = 0.1$ is so that the effects of convexity can be seen clearly.

A log–log plot of $\max |\partial u/\partial r|$ as $t \rightarrow t_*^-$ for the two isotherms $\sigma _{\textit{conca}v\textit{e}}$ (blue) and $\sigma _{\textit{con}v\textit{ex}}$ (pink) is presented in figure 6. Recall that the thin-film equations (solid curves) solved are given by (2.1), (1.3b ) and (1.3c ) and the similarity solution prediction (dashed lines) for $\max |\partial u/\partial r|$ is given by (2.21). The results shown in figure 6 demonstrate that, indeed, the similarity solution (2.21) accurately captures the nonlinear effects of the surface tension isotherm. The prediction $\max |\partial u/\partial r|=(t_*-t)^{-1}$ (black dashed) is the prediction for a locally linear isotherm ( $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)_{\varGamma = \varGamma _*}=0$ ) and so we can see from figure 6 that, as expected, local concavity $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)_{\varGamma = \varGamma _*}\lt 0$ weakens the shock and local convexity $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)_{\varGamma = \varGamma _*}\gt 0$ strengthens the shock, by changing the prefactor of the blow up of $\max |\partial u/\partial r|$ as $t \rightarrow t_*^-$ .

Figure 6. Effects of the nonlinearity of surface tension isotherms. Log–log plot of $\max |\partial u/\partial r|$ as predicted by solution of the thin-film (T-F) equations (2.1), (1.3b ), (1.3c ) for an example concave isotherm as described in the text (blue solid). Analogous to an example convex isotherm as described in the text (pink solid). Initial surfactant distribution is $\varGamma _i = e^{-r^2}$ . Corresponding similarity solution (S) predictions (2.21) are shown in dashed lines for the concave (blue dashed) and convex (pink dashed) cases. For comparison, the line $\max |\partial u/\partial r |=(t_*-t)^{-1}$ is shown also (black dashed). Inset: magnified view of the same log–log plots.

For the initial surfactant deposition $\varGamma _i = e^{-r^2}$ considered in this section, the effects of the change of $\max |\partial u/\partial r|$ are moderate. For the concave example discussed, $\max |\partial u/\partial r|\approx 0.85(t_*-t)^{-1}$ as $t \rightarrow t_*$ , i.e. only a 15 % decrease of the shock strength, which is due to the term $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)|_{\varGamma = \varGamma _*}\varGamma _*^2/(h_*V^2)$ being small in magnitude in (2.21). For the convex example discussed, $\max |\partial u/\partial r|\approx 1.03(t_*-t)^{-1}$ as $t \rightarrow t_*$ . Different initial surfactant deposition profiles $\varGamma _i$ may give stronger effects. In particular, it might be mathematically interesting to consider what would happen if $(\textrm{d}^2 \sigma /\textrm{d} \varGamma ^2)|_{\varGamma = \varGamma _*}\varGamma _*^2/(h_*V^2) = 1$ , where the derivation given in § 2.2 would have to be amended to account for higher-order terms, since the nonlinear term in (2.15) would vanish.

3.1. Scalings

The scaling for the outer regions is the same for the IMCE, IMC and IME regimes and may be derived as follows. Subscripts are used to denote the region under discussion. Let $\Delta r$ denote the characteristic width of a given region. In region III, balancing inertia with the capillary stress gradient and/or viscous extensional stress, $u_{\text{III}} t^{-1} \sim h_{\text{III}}(\Delta r_{\text{III}})^{-3}$ and/or $u_{\text{III}}(\Delta r_{\text{III}})^{-2}$ . Additionally, considering conservation of mass, and far-field matching, $h_{\text{III}} \rightarrow 1$ , gives that $u_{\text{III}}\sim t^{-(1/2)}$ , $h_{\text{III}}\sim 1$ , $\Delta r_{\text{III}} \sim t^{{1}/{2}}$ . Integrating across region II, where there is a quasi-static balance between Marangoni, capillary and/or extensional stress (analogous to Rankine–Hugoniot jump conditions), one deduces that $\varGamma _I \sim h_{\text{III}}^2(\Delta r_{\text{III}})^{-2}$ and/or $h_{\text{III}} u_{\text{III}}(\Delta r_{\text{III}})^{-1}$ , which implies that $\varGamma _I \sim t^{-1}$ . Since $\int _0^{\infty }\varGamma r \textrm{d}r$ is constant, $\Delta r_{\text{I}} \sim t^{{1}/{2}}$ . Conservation of mass then gives $u_{\text{I}}\sim \Delta r_{\text{I}} t^{-1}$ and hence $u_{\text{I}}\sim t^{-(1/2)}$ . Since $\varGamma h^{-1}$ is conserved in Lagrangian coordinates (D3), then $h_{\text{I}}\sim t^{-1}$ . In summary,

(3.4) \begin{equation} (u_{\text{I}},h_{\text{I}},\varGamma _{\text{I}},\Delta r_{\text{I}})\sim (t^{-\frac{1}{2}},t^{-1},t^{-1},t^{\frac {1}{2}}) ,\end{equation}

and

(3.5) \begin{equation} (u_{\text{III}},h_{\text{III}},\Delta r_{\text{III}})\sim (t^{-\frac{1}{2}},1,t^{\frac {1}{2}}) ,\end{equation}

with $\varGamma _{\text{III}}=0$ .

The summary of the dominant balance in each of the regions is as follows. The balance will be more precisely discussed in the derivation. Region I is where the surfactant concentration is uniform (to leading order), as there are no terms which would balance the Marangoni stress arising from a non-uniform surfactant concentration. In the IMCE and IMC regimes, the Marangoni stress balances the capillary stress in region II. In the IME regime, the Marangoni stress balances the viscous extensional stress in region II. Inertia does not appear in region II due primarily to the presence of sharp gradients (e.g. the capillary term in (3.1) has third-order spatial derivatives). Region III is where inertia balances capillary stress and/or viscous extensional stress, as there are no surfactants.

3.2. Region I

Expecting a similarity solution, region I is given by $0 \leqslant r t^{-(1/2)} \lt \eta _{\!f}$ for some $\eta _{\!f} = \eta _{\!f}(\mathcal{M}, {\textit{Re}})$ . The similarity coordinate is then given by $\eta :=\eta _{\!f}^{-1}rt^{-(1/2)} \in [0,1)$ . With the scalings of region I, conservation of momentum (any one of (3.1–3.3)) to leading order gives

(3.6) \begin{equation} \frac {\partial \varGamma _{\text{I}}}{\partial r}=0. \end{equation}

The solution to (3.6), (1.3b ) and (1.3c ) using global conservation of surfactant $\int _0^{\infty }\varGamma r \textrm{d}r = {1}/{2}$ shows that

(3.7) \begin{equation} u_{\text{I}}=\eta _{\!f} t^{-\frac{1}{2}}\frac {\eta }{2}, \, h_{\text{I}}=\eta _{\!f}^{-2}t^{-1}f(\eta ), \, \varGamma _{\text{I}}=\eta _{\!f}^{-2}t^{-1}, \end{equation}

where $\eta \in [0,1)$ and $f$ is given by (using (D4) and (D5))

(3.8) \begin{equation} 1=\varGamma _i\left (\left (2\int _0^{\eta } \eta 'f(\eta ') \textrm{d} \eta '\right )^{\frac {1}{2}}\right )f(\eta ). \end{equation}

Since $f$ is minimal where $\varGamma _i$ is maximal (3.8), it then follows that $\min f = 1$ and hence

(3.9) \begin{equation} h_{\textit{min}}=\eta _{\!f}^{-2}t^{-1}. \end{equation}

Finally, in order for the thickness profile $h_{\text{I}}$ in region I to match onto region II, $h_{\text{I}}$ is singular as the end of region I is approached, $\eta \rightarrow 1^-$ . By considering the limits $\eta \rightarrow 1^-$ of (3.8), it can be shown that $f \sim (1-\eta )^{-1}$ as $\eta \rightarrow 1^-$ (see Appendix D.1). Then, $h_{\text{I}}\sim t^{-1}(1-\eta )^{-1}= \eta _{\!f} t^{-(1/2)}(\eta _{\!f} t^{{1}/{2}}-r)^{-1}$ as $\eta \rightarrow 1^-$ .

3.3. Region II

Region II is the transition region between regions I and III, where regularisation of the surfactant front occurs. In order to analyse region II, it is useful to consider the reference frame of the shock propagating at spatial location $r_{\!f} =\eta _{\!f} t^{{1}/{2}}$ , with the spatial coordinate given by $\Delta r_{\text{II}}:=r_{\text{II}}-\eta _{\!f} t^{ {1}/{2}}$ and velocity given by $\Delta u_{\text{II}}:=u_{\text{II}}-(1/2)\eta _{\!f} t^{-(1/2)}$ .

The scalings for region II may be derived as follows. By using conservation of mass and surfactant (1.3b , 1.3c ) and matching $h$ (note $h_{\text{I}} \sim \eta _{\!f} t^{-(1/2)}(\eta _{\!f} t^{{1}/{2}}-r)^{-1}$ as $r\eta _{\!f}^{-1}t^{-(1/2)} \rightarrow 1^-$ ) and $\varGamma$ onto region I,

(3.10) \begin{equation} (\Delta u_{\text{II}}, h_{\text{II}},\varGamma _{\text{II}}) \sim (\Delta r_{\text{II}}t^{-1}, (\Delta r_{\text{II}})^{-1}t^{-\frac{1}{2}}, t^{-1}). \end{equation}

Then, we may compare the different terms in the momentum balance (3.1–3.3). First, the inertial term has order $\textit {O}(\Delta u_{\text{II}} t^{-1}) \sim (\Delta r_{\text{II}})t^{-2} \ll t^{-(3/2)}$ (as $\Delta r_{\text{II}}$ is a smaller than $\Delta r_{\text{I}} \sim t^{{1}/{2}}$ ), the Marangoni stress term has order $\textit {O}(\varGamma _{\text{II}} (h_{\text{II}} \Delta r_{\text{II}})^{-1})\sim t^{-(1/2)}$ , the capillary stress term has order $\textit {O}(h_{\text{II}} (\Delta r_{\text{II}})^{-3})\sim (\Delta r_{\text{II}})^{-4}t^{-(1/2)}$ and the viscous extensional stress term has order $\textit {O}(\Delta u_{\text{II}} (\Delta r_{\text{II}})^{-2})\sim (\Delta r_{\text{II}})^{-1}t^{-1}$ . In particular, the inertial term does not appear in region II to leading order.

Then, in the IMC regime, the leading-order momentum balance in region II is between Marangoni and capillary stresses. Hence, the above scalings give $\Delta r_{\text{II}} \sim 1$ . In the IME regime, the leading-order momentum balance in region II must be between Marangoni and extensional stresses, and hence $\Delta r_{\text{II}} \sim t^{-(1/2)}$ . Finally, in the IMCE regime, the leading-order momentum balance in region II is between Marangoni and capillary stresses (since balancing Marangoni and extensional stresses gives $\Delta r_{\text{II}} \sim t^{-(1/2)}$ , which gives the capillary stress term $\sim t^{{3}/{2}} \gg t^{-1}$ ), and hence $\Delta r_{\text{II}} \sim 1$ . In particular, the long-wavelength approximation will fail to hold for the IME regime as $t \rightarrow \infty$ since $\Delta r_{\text{II}}\rightarrow 0$ (see § 5 for a discussion).

3.3.1. First case (IMCE, IMC regimes): Marangoni–capillary stress balance

As above, the scalings discussed for the IMCE and IMC regime are given by

(3.11) \begin{equation} (\Delta u_{\text{II}}, h_{\text{II}}, \varGamma _{\text{II}}, \Delta r_{\text{II}}) \sim (t^{-1},t^{-\frac{1}{2}},t^{-1},1). \end{equation}

Then, the momentum balance (3.1) and (3.2) becomes, at leading order,

(3.12) \begin{equation} \frac {2}{h_{\text{II}}}\frac {\partial \varGamma _{\text{II}}}{\partial \Delta r_{\text{II}}}=\frac {1}{2\mathcal{M}}\frac {\partial ^3 h_{\text{II}}}{\partial (\Delta r_{\text{II}})^3} ,\end{equation}

which integrates with respect to $\Delta r_{\text{II}}$ to give

(3.13) \begin{equation} \varGamma _{\text{II}}-\eta _{\!f}^{-2}t^{-1}=\frac {1}{4\mathcal{M}}\left (h_{\text{II}}\frac {\partial ^2 h_{\text{II}}}{\partial (\Delta r_{\text{II}})^2}-\frac {1}{2}\left (\frac {\partial h_{\text{II}}}{\partial \Delta r_{\text{II}}}\right )^2\right )\!, \end{equation}

where the spatially constant term $\eta _{\!f}^{-2}t^{-1}$ is deduced by matching onto the region I solution (3.7). In particular, since $\varGamma _{\text{III}}=0$ (region III does not contain surfactants),

(3.14) \begin{equation} -\eta _{\!f}^{-2}t^{-1}=\lim _{\text{II}\rightarrow \text{III}}\frac {1}{4\mathcal{M}}\left (h_{\text{II}}\frac {\partial ^2 h_{\text{II}}}{\partial (\Delta r_{\text{II}})^2}-\frac {1}{2}\left (\frac {\partial h_{\text{II}}}{\partial \Delta r_{\text{II}}}\right )^2\right )\!, \end{equation}

where $\lim _{\text{II}\rightarrow \text{III}}$ denotes the limit as region III is approached from region II.

3.3.2. Second case (IME regime): Marangoni–extensional stress balance

As explained, the scalings in the IME regime are given by

(3.15) \begin{equation} (\Delta u_{\text{II}}, h_{\text{II}}, \varGamma _{\text{II}}, \Delta r_{\text{II}}) \sim (t^{-\frac {3}{2}},1,t^{-1},t^{-\frac{1}{2}}). \end{equation}

Then, the momentum balance (3.3) becomes, at leading order,

(3.16) \begin{equation} \frac {2}{h_{\text{II}}}\frac {\partial \varGamma _{\text{II}}}{\partial \Delta r_{\text{II}}}=\frac {4}{{\textit{Re}} } \left (\frac {\partial ^2 \Delta u_{\text{II}}}{\partial \Delta r_{\text{II}}^2}+\frac {1}{h_{\text{II}}}\frac {\partial \Delta u_{\text{II}}}{\partial \Delta r_{\text{II}}}\frac {\partial h_{\text{II}}}{\partial \Delta r_{\text{II}}}+\frac {1}{4h_{\text{II}}t}\frac {\partial h_{\text{II}}}{\partial \Delta r_{\text{II}}}\right )\! , \end{equation}

which integrates with respect to $\Delta r_{\text{II}}$ as

(3.17) \begin{equation} \varGamma _{\text{II}} - \eta _{\!f}^{-2}t^{-1}=\frac {1}{{\textit{Re}}} \left (2h_{\text{II}}\frac {\partial \Delta u_{\text{II}}}{\partial \Delta r_{\text{II}}}+\frac {h_{\text{II}}}{2t}\right )\!. \end{equation}

In particular,

(3.18) \begin{equation} - \eta _{\!f}^{-2}t^{-1}=\lim _{\text{II}\rightarrow \text{III}} \frac {1}{{\textit{Re}}} \left (2h_{\text{II}}\frac {\partial \Delta u_{\text{II}}}{\partial \Delta r_{\text{II}}}+\frac {h_{\text{II}}}{2t}\right )\!, \end{equation}

where $\lim _{\text{II}\rightarrow \text{III}}$ denotes the limit as region III is approached from region II.

3.4. Region III

Region III is the region ahead of the front where there is no surfactants and hence $\varGamma = 0$ . Due to the need to match onto region II, § 3.3, the derivation of the structure of region III also has to consider two different cases. The first case considers the IMCE and IMC regimes, which contain capillary stress terms, where the IMC similarity solutions can directly be seen to be obtainable from the IMCE similarity solutions in the limit ${\textit{Re}} \gg 1$ . The second case considers the IME regime, which does not contain a capillary stress term; as shown in § 3.3 the absence of these higher derivatives changes the details of region II and hence the matching onto region III. Not containing the curvature gradient term $\partial ^3 h/\partial r^3$ also changes the order of the ordinary differential equations (ODEs) of the similarity solutions. Thus, obtaining the IME regime similarity solutions directly from the IMCE similarity solutions is not trivial, but is shown to indeed be the case in § 4.

3.4.1. First case (IMCE, IMC regimes)

The similarity ansatz is given by

(3.19) \begin{equation} u_{\text{III}}=\eta _{\!f}t^{-\frac{1}{2}}U(\eta ), \,h_{\text{III}}=H(\eta ), \end{equation}

where $\eta \in (1,\infty )$ . Substitution of the self-similarity ansatz (3.19) into the IMCE momentum (3.1) gives

(3.20) \begin{align} -\frac {1}{2}U-\frac {1}{2}\eta \frac {\textrm{d}U}{\textrm{d}\eta }+U\frac {\textrm{d}U}{\textrm{d}\eta } &=\frac {1}{2 \eta _{\!f}^4\mathcal{M}}\frac {d }{\textrm{d} \eta }\left (\frac {1}{\eta }\frac {d}{\textrm{d} \eta }\left (\eta \frac {\textrm{d} H}{\textrm{d}\eta }\right )\right )\nonumber \\ &\quad + \frac {4}{{\textit{Re}} \eta _{\!f}^2}\frac {1}{H}\left (\frac {d}{\textrm{d}\eta }\left (\frac {H}{\eta }\frac {d}{\textrm{d} \eta }\left (\eta U\right )\right ) - \frac {1}{2}\frac {U}{\eta }\frac {\textrm{d}H}{\textrm{d}\eta }\right )\!. \end{align}

Then, defining $(J,K):=(\textrm{d}H/\textrm{d}\eta , \textrm{d}^2H/\textrm{d}\eta ^2)$ and using (1.3b ) gives the similarity solutions for the IMCE regime as a system of ODEs for $(U,H,J,K)$

(3.21a) \begin{align} \frac {\textrm{d}U}{\textrm{d}\eta } &=\left (\frac {\eta }{2}-U\right )\frac {J}{H}-\frac {U}{\eta }, \end{align}

(3.21b) \begin{align} \frac {\textrm{d}H}{\textrm{d}\eta }&=J, \end{align}

(3.21c) \begin{align} \frac {\textrm{d}J}{\textrm{d}\eta } &= K, \end{align}

(3.21d) \begin{align} \frac {\textrm{d}K}{\textrm{d}\eta } &= -\frac {K}{\eta }+\frac {J}{\eta ^2}+2\eta _{\!f}^4\mathcal{M}\left (-\frac {U}{2}+\left (U-\frac {\eta }{2}\right )\left (\left (\frac {\eta }{2}-U\right )\frac {J}{H}-\frac {U}{\eta }\right )\right ) \nonumber \\ &\quad -\frac {8\eta _{\!f}^2\mathcal{M}}{{\textit{Re}} }\left (\left (\frac {1}{2}+\frac {U}{2\eta }-\left (\frac {\eta }{2}-U\right )\frac {J}{H}\right )\frac {J}{H} +\left (\frac {\eta }{2}-U\right )\frac {K}{H}\right )\!, \end{align}

where (3.21d ) is obtained by repeatedly using (3.21a ) to eliminate derivatives of $U$ in (3.20). By matching onto region II, the behaviour near the left boundary condition at $\eta =1+\delta$ for $\delta \ll 1$ can be deduced to satisfy (see Appendix E.1)

(3.22a) \begin{align} U(1+\delta ) &= \frac {1}{2} +\dots , \end{align}

(3.22b) \begin{align} H(1+\delta ) &= \sqrt {8\mathcal{M}}\delta - \frac {\eta _{\!f}^2\mathcal{M}}{{\textit{Re}} }\delta ^2 \log \delta + q\delta ^2 +\dots , \end{align}

(3.22c) \begin{align} J(1+\delta ) &= \sqrt {8 \mathcal{M}}-\frac {2\eta _{\!f}^2 \mathcal{M}}{{\textit{Re}} }\delta \log \delta - \frac {\eta _{\!f}^2 \mathcal{M}}{{\textit{Re}} }\delta +2q\delta +\dots , \end{align}

(3.22d) \begin{align} K(1+\delta ) &= -\frac {2\eta _{\!f}^2 \mathcal{M}}{{\textit{Re}} }\log (\delta ) - \frac {3\eta _{\!f}^2 \mathcal{M}}{{\textit{Re}} }+2q + \dots , \end{align}

for some constant $q$ .

In summary, the similarity solution for region III can be obtained numerically via a shooting algorithm as follows. First, guess some values $\eta _{\!f},q$ . Then, consider the left boundary condition (3.22) to be at $\eta =1+\delta$ for some $\delta \ll 1$ (chosen small enough so that the result is independent of $\delta$ ). Then, integrate $\eta \rightarrow \infty$ subject to (3.21). The process is repeated by adjusting the two constants $\eta _{\!f}$ and $q$ towards satisfying the two constraints

(3.23) \begin{equation} H(\infty )= 1 \text{ and } \int _1^{\infty }(H-1)\eta \textrm{d}\eta =\frac {1}{2}, \end{equation}

where the second constraint is from global conservation of mass (see Appendix E.1), which shows that the mass that was located in region I at $t = 0$ must be in region III.

The IMC similarity solutions may be obtained directly from the IMCE similarity solutions in the limit that ${\textit{Re}} \gg 1$ , thereby omitting the ${\textit{Re}}^{-1}$ terms (3.21–3.23); see Appendix F for the exact details.

3.4.2. Second case (IME regime)

The similarity ansatz for the IME regime is still given by (3.19) and when substituted into (3.3) gives

(3.24) \begin{equation} -\frac {1}{2}U-\frac {1}{2}\eta \frac {\textrm{d}U}{\textrm{d}\eta }+U\frac {\textrm{d}U}{\textrm{d}\eta } = \frac {4}{{\textit{Re}} \eta _{\!f}^2}\frac {1}{H}\left (\frac {d}{\textrm{d}\eta }\left (\frac {H}{\eta }\frac {d}{\textrm{d} \eta }\left (\eta U\right )\right ) - \frac {1}{2}\frac {U}{\eta }\frac {\textrm{d}H}{\textrm{d}\eta }\right )\!. \end{equation}

Then, letting $J=\textrm{d}H/\textrm{d}\eta$ and using (1.3b ) gives the similarity solutions for the IME regime as a system of ODEs for $(U,H,J)$

(3.25a) \begin{align} \frac {\textrm{d}U}{\textrm{d}\eta } &= \left (\frac {\eta }{2}-U\right )\frac {J}{H}-\frac {U}{\eta }, \end{align}

(3.25b) \begin{align} \frac {\textrm{d}H}{\textrm{d}\eta }&=J, \end{align}

(3.25c) \begin{align} \frac {\textrm{d}J}{\textrm{d}\eta }&=\left (\frac {1+\frac {U}{\eta }}{2\left (U-\frac {1}{2} \eta \right )}+\frac {J}{H}\right )J \nonumber \\ &\quad + \frac {{\textit{Re}} \eta _{\!f}^2}{4}\left (\frac {HU}{2\left (U-\frac {1}{2} \eta \right )} +\left (U-\frac {\eta }{2}\right )J+\frac {HU}{\eta }\right )\!\!, \end{align}

where (3.25a ) is used repeatedly to eliminate derivatives of $U$ in (3.24). By matching onto region II, the behaviour near the left boundary condition at $\eta =1+\delta$ for $\delta \ll 1$ can be deduced to satisfy (see Appendix E.2)

(3.26a) \begin{align} U(1+\delta ) &= \frac {1}{2}-\frac {1}{2}\delta +\dots , \end{align}

(3.26b) \begin{align} H(1+\delta ) &= 2{\textit{Re}}\eta _{\!f}^{-2}-4\tilde {q}\delta ^{\frac {1}{4}} +\dots , \end{align}

(3.26c) \begin{align} J(1+\delta ) &= -\tilde {q}\delta ^{-\frac {3}{4}}+\dots , \end{align}

for some constant $\tilde {q}$ to be found. For the IME regime, it is difficult to see from (3.25, 3.26) alone that the IME regime may be recovered upon taking the limit $\mathcal{M} \rightarrow \infty$ of the IMCE regime (3.22), which makes sense, given that the regularisation mechanism is fundamentally different in between the IME (§ 3.3.2) and IMCE regimes (§ 3.3.1). However, it is shown in § 4 that the IME regime may in fact be obtained from the IMCE regime in the limit $\mathcal{M} \rightarrow \infty$ .

In summary, the similarity solution for region III can once again be found via the shooting algorithm, for two parameters $\eta _{\!f}, \tilde {q}$ , with the far-field conditions (3.23).

3.5. Initial surfactant concentration distribution with polynomial decay

In this subsection, we discuss the case where the initial surfactant distribution $\Gamma_i(r) \sim r^{-\alpha}$ as $r \rightarrow \infty$ . Since the amount of added surfactant is finite, it follows that $\alpha \gt2$ (consider $\int_0^{\infty}r \Gamma_i(r) dr$ ). The assumption of $\Gamma_i(r)$ decaying faster than polynomial corresponds to $\alpha = \infty$ in the below expressions, which recovers the exact expressions in the rest of the paper. Note that the arguments by Eshima et al. (Reference Eshima, Stone and Deike2025) inadvertently had only dealt with the case $\alpha = \infty$ .

From (3.7), we have that the thickness in region I is given by $h = \eta_f^{-2}t^{-1}f(\eta)$ with $f$ given by (3.8). Then, $f \sim (1-\eta)^{-1-2(\alpha-2)^{-1}}$ as $\eta \rightarrow 1^-$ (analogous method as in Appendix D). Then, the scaling of region II (3.10) becomes

(3.27) \begin{align} (\Delta u_{\text{II}}, h_{\text{II}},\Gamma_{\text{II}}) \sim \big(\Delta r_{\text{II}}t^{-1}, |\Delta r_{\text{II}}|^{-1-2(\alpha-2)^{-1}}t^{-\frac{1}{2}+(\alpha-2)^{-1}}, t^{-1}\big). \end{align}

By considering dominant force balance arguments in region II, as in § 3.3, it can be shown that the dominant balance in region II for the IMC and IMCE regimes is still given by Marangoni and capillary stresses with (3.11) modified to

(3.28) \begin{align} (\Delta u_{\text{II}},h_{\text{II}},\Gamma_{\text{II}},\Delta r_{\text{II}})\sim \big(t^{-1+(2\alpha-2)^{-1}},t^{-\frac{1}{2}+(2\alpha-2)^{-1}},t^{-1},t^{(2\alpha-2)^{-1}}\big). \end{align}

Similarly, the dominant balance in region II for the IME regime is still given by Marangoni and extensional stresses with (3.15) modified to

(3.29) \begin{align} (\Delta u_{\text{II}},h_{\text{II}},\Gamma_{\text{II}},\Delta r_{\text{II}})\sim \big(t^{-\frac{3}{2}+2\alpha^{-1}},1,t^{-1},t^{-\frac{1}{2}+2\alpha^{-1}}\big).\end{align}

Since the dominant force balance in region II does not change, the solutions shown in the text still follows (see the solution column in table 1).

4.1. Comparison with numerical solutions of the thin-film equations

In this section, we verify and discuss the similarity solutions obtained in § 3 using the thin-film equations. The verification is done by checking that the similarity solutions accurately capture the three physical predictions from the thin-film equations in the late-time limit: i) film thinning ( $t^{-1}$ ), ii) front propagation ( $t^{{1}/{2}}$ ) and iii) capillary wave characteristics.

In figure 8 we compare the IMCE thin-film equations (solid curves) with the similarity solutions (dashed curves), where (a,c,e,g) show the time evolution of the minimum thickness $h_{\textit{min}}$ and (b,d, f,h) show the time evolution of the location of the surfactant front $r_{\!f}$ for $(\mathcal{M},{\textit{Re}}) \in \{0.1,1,10\} \times \{1,4,10,40\}$ . As previously derived, the similarity solution predicts $h_{\textit{min}}=\eta _{\!f}^{-2}t^{-1}$ and $r_{\!f} = \eta _{\!f} t^{{1}/{2}}$ . As expected, for late times, $t \gg 1$ , the agreement between the similarity solution and the thin-film equations is good. Similarly, in figure 9 we compare the thickness profile for the thin-film equations at $t = 1000$ (solid curves) with the similarity solutions (dotted curves). The horizontal and vertical axes are given by $rt^{-(1/2)}$ and $h$ . Figure 9(a–l) shows the result of a parameter sweep for $(\mathcal{M},{\textit{Re}}) \in \{0.1,1,10\} \times \{1,4,10,40\}$ . The agreement between the thin-film equations and the similarity solution is once again robust for the whole parameter space. There are no fitting parameters used in figures 8 and 9 and the similarity solutions therefore capture the desired physical results accurately. The same verification between the thin-film equations and the similarity solutions may be done for the IMC and IME equations and are given in Appendix G.

Figure 8. Comparison of the time evolution (a,c,e,g) of the minimum thickness $h_{\textit{min}}$ and (b,d, f,h) the location of the surfactant front $r_{\!f}$ , when comparing the similarity solution (dashed), and the IMCE thin-film equations (solid). For reference, the similarity solution predicts $h_{\textit{min}}=\eta _{\!f}^{-2}t^{-1}$ and $r_{\!f} = \eta _{\!f} t^{{1}/{2}}$ . The comparison is obtained for $(\mathcal{M}, {\textit{Re}}) \in \{0.1,1,10\} \times \{1,4,10,40\}$ . The curves are coloured according to $\mathcal{M} = 0.1,1,10$ (purple, orange and blue–green, respectively) and shaded according to ${\textit{Re}}=1,4,10,40$ (light to dark). As expected, for $t \gg 1$ , the numerical solutions to the thin-film equation approach the similarity solutions.

Figure 9. Comparison of the thickness profiles predicted by the thin-film equations at $t = 1000$ (solid curves) and the similarity solutions (dotted curves) for various $\mathcal{M}, {\textit{Re}}$ . Thin-film equations are solved for $\varGamma _i = e^{-r^2}$ . The horizontal and vertical axes are given by $rt^{-(1/2)}$ and $h$ , where $r$ is the radial coordinate, $t$ is time and $h$ is the thickness. The curves are coloured according to $\mathcal{M} = 0.1,1,10$ (purple, orange, blue–green, respectively) and shaded according to ${\textit{Re}}=1,4,10,40$ (light to dark). (a–l) Parameter sweep for $(\mathcal{M}, {\textit{Re}}) \in \{0.1,1,10\}\times \{1,4,10,40\}$ .

4.2. Self-consistency of the similarity solutions

We also verify that the IMC and IME regimes are obtained from expected limits of the IMCE regime. Since all three regimes have a front at radial coordinate $r_{\!f}(t)=\eta _{\!f} t^{{1}/{2}}$ as $t \rightarrow \infty$ , we may compare the value of the prefactor $\eta _{\!f}$ . Figure 10(a) shows $\eta _{\!f}(\mathcal{M},{\textit{Re}})$ as obtained from the similarity solution ODEs in the IMCE regime (see § 3.4.1) versus ${\textit{Re}}$ for $\mathcal{M}=0.1,1,10$ (purple, orange, blue–green solid curves) and compares with $\eta _{\!f}(\mathcal{M}=0.1), \eta _{\!f}(\mathcal{M}=1), \eta _{\!f}(\mathcal{M}=10)$ (purple, orange, blue–green dotted lines) as obtained from the similarity solution ODEs in the IMC regime (see Appendix F). Panel (b) shows $\eta _{\!f}(\mathcal{M},{\textit{Re}})$ as obtained from the similarity solution ODEs in the IMCE regime (see § 3.4.1) versus $\mathcal{M}$ for ${\textit{Re}}=1,4,10$ (light grey, grey, dark grey solid curves) and compares with $\eta _{\!f}({\textit{Re}}=1), \eta _{\!f}({\textit{Re}}=4), \eta _{\!f}({\textit{Re}}=10)$ (light grey, grey, dark grey dotted lines) as obtained from the similarity solution ODEs in the IME regime (see § 3.4.2). From figure 10(a), the IMCE regime prediction for $\eta _{\!f}$ agrees well with the IMC regime prediction in the limit ${\textit{Re}} \rightarrow \infty$ ; similarly, panel (b) shows that the IMCE regime prediction for $\eta _{\!f}$ agrees well with the IME regime prediction in the limit $\mathcal{M} \rightarrow \infty$ . In other words, the relation

(4.1) \begin{equation} \lim _{{\textit{Re}}\rightarrow \infty } \eta _{\!f}^{\textit{IMCE}}(\mathcal{M}, {\textit{Re}}) = \eta _{\!f}^{\textit{IMC}}(\mathcal{M}) \text{ and } \lim _{\mathcal{M}\rightarrow \infty } \eta _{\!f}^{\textit{IMCE}}(\mathcal{M},{\textit{Re}}) = \eta _{\!f}^{\textit{IME}}({\textit{Re}}) \end{equation}

has been verified (although no formal proof is given).

Figure 10. Comparison of IMC and IME regimes as limits of the IMCE regime. (a) Solid curves show $\eta _{\!f}(\mathcal{M},{\textit{Re}})$ as obtained from the ODEs found for the similarity solutions in the IMCE regime (see § 3.4.1) against ${\textit{Re}}$ for $\mathcal{M}=0.1,1,10$ (purple, orange, blue–green). Dotted lines show $\eta _{\!f}(\mathcal{M}=0.1), \eta _{\!f}(\mathcal{M}=1), \eta _{\!f}(\mathcal{M}=10)$ as obtained from the similarity solution ODEs in the IMC regime (see Appendix F) . (b) Solid curves show $\eta _{\!f}(\mathcal{M},{\textit{Re}})$ as obtained from the similarity solution ODEs in the IMCE regime against $\mathcal{M}$ for ${\textit{Re}}=1,4,10$ (light grey, grey, dark grey). Dotted lines show $\eta _{\!f}({\textit{Re}}=1), \eta _{\!f}({\textit{Re}}=4), \eta _{\!f}({\textit{Re}}=10)$ (light grey, grey, dark grey) as obtained from the similarity solution ODEs in the IME regime (see § 3.4.2).

The late-time similarity solution has been found for inertial surfactant deposition with or without capillary stress and/or extensional stress in the momentum balance, by considering the three possible dynamical regimes IMCE, IMC and IME.

The scalings of the regions I and III away from the surfactant front are the same for all three regimes. However, the details differ in both regions I and III. For the uniform surfactant region I, the value of the coefficient $\eta _{\!f}$ , which sets the location of the surfactant front ( $r_{\!f} = \eta _{\!f} t^{{1}/{2}}$ ) and the minimum thickness ( $h_{\textit{min}} = \eta _{\!f}^{-2} t^{-1}$ ), is dependent on the precise value of $\mathcal{M}$ and ${\textit{Re}}$ . For the capillary wave region III, as ${\textit{Re}}$ increases, there is less viscous damping of the capillary waves and as $\mathcal{M}$ increases, the thickness profile bump becomes sharper, which makes sense since the Marangoni stress is stronger (see figure 9). The IMC regime and IME regime can be obtained as limits of the IMCE regime (§ 4.1).