2.1 Formulation of system for similarity solution

While the system can be simulated as above for any $\epsilon $ , a chief advantage of considering the simplified model (1.5) with initial condition (1.6) is that it permits a similarity solution. This similarity solution may be calculated numerically by solving a single boundary value problem with much greater precision than the original PDE, which will be particularly valuable when looking at small values of the height $\epsilon $ in front of the initial discontinuity. We also note that, while exact for the simplified model (1.5), this similarity solution also represents the early-time behaviour of more complex cross-diffusion models of the form (1.3) evolving from an initial condition (1.6), and so captures essential behaviour of this class of models as well (see Section 4 for further discussion).

To construct similarity solutions, we first define the similarity variable $\eta $ , and similarity profiles f and g by

(2.1) $$ \begin{align} \eta = \frac{x-x^*_0}{t}, \quad u(x,t) = tf(\eta), \quad v(x,t) = g(\eta). \end{align} $$

The factors t in (2.1) result from requiring the terms in (1.5) to have the same dependence on t, while the scale of v cannot change as it must match to the fixed values behind and ahead of the initial discontinuity in the large $\eta $ limit. From (1.5), the profiles f and g then satisfy the ordinary differential equation (ODE) system

(2.2) $$ \begin{align} f - \eta \frac{{d} f}{{d} \eta} = g, \quad -\eta \frac{{d} g}{{d} \eta} = \frac{{d} }{{d} \eta}\bigg(f\frac{{d} g}{{d} \eta}\bigg), \quad -\infty < \eta < \infty \end{align} $$

with far-field conditions consistent with the initial condition (1.6):

(2.3) $$ \begin{align} f, g \to \begin{cases} 1, & \eta \to -\infty , \\ \epsilon, & \eta \to \infty. \end{cases} \end{align} $$

As $t\to 0^+$ , the profiles contract around the point $x=x^*_0$ ; thus, a similarity solution satisfying these conditions will approach the correct initial condition (1.6).

Looking closer at the far-field conditions (2.3), if we expand f and g close to a constant in each far-field $|\eta |\to \infty $ , we find that (2.2) supports the behaviour

(2.4) $$ \begin{align} f \sim C_0 + C_1\eta + C_2(\eta) \exp\bigg(-\frac{\eta^2}{2C_0}\bigg), \quad |\eta| \to\infty , \end{align} $$

for $C_0$ , $C_1$ constants, and the Wentzel–Kramers–Brillouin (WKB) prefactor $C_2$ is a function of $\eta $ (whose exact form is unimportant). To satisfy the far-field conditions, we require $C_0 = 1$ for ${\eta \to -\infty }$ and $C_0 = \epsilon $ for $\eta \to \infty $ , while the algebraically growing mode $C_1$ must be zero in both cases. The remaining WKB mode decays exponentially in each case. Thus, each far-field condition effectively applies two constraints on the system, resulting in four boundary conditions for the third-order system (2.2).

To find solutions then, an additional degree of freedom is required. This is provided by the fact that the system (2.2) is singular at the point at $\eta = 0$ . Expanding f and g near $\eta =0$ , we find that f may be nondifferentiable; in particular,

(2.5) $$ \begin{align} f = a_0 + a_1\eta\log|\eta| + a_2\eta + \cdots, \quad g = a_0 - a_1\eta + \cdots, \end{align} $$

where $a_0$ , $a_1$ , $a_2$ are arbitrary, and all higher coefficients may be determined in terms of these three constants. For g to be differentiable across the singular point, $a_0$ and $a_1$ must take the same value on either side of $\eta =0$ . However, $a_2$ may take different values on either side, say

$$ \begin{align*} a_2 = \begin{cases} a_2^-, & \eta < 0, \\ a_2^+, & \eta> 0. \end{cases} \end{align*} $$

This extra degree of freedom allows us to satisfy both far-field conditions. A consequence of the singular point is that even for time greater than zero, u will be nondifferentiable at the point $x_0$ at which the initial discontinuity existed. The PDE system (1.5) permits this because there is no diffusion of u.

2.2 Numerical calculation of similarity solutions

We now describe the method of numerically calculating solutions of (2.2). To achieve high accuracy and observe the change in profiles as $\epsilon $ is varied, we use the numerical continuation package AUTO-07p [Reference Doedel, Champneys, Dercole, Fairgrieve, Kuznetsov, Oldeman, Paffenroth, Sandstede, Wang and Zhang13]. To obtain a solution with the appropriate boundary conditions (2.4), we will apply a homotopy continuation method, in which an initial trial solution is computed with explicit error terms in the boundary conditions, which are then each reduced to zero in turn via numerical continuation. This process has been applied to the computation of similarity solutions previously with much success [Reference Dallaston, Tseluiko, Zheng, Fontelos and Kalliadasis11, Reference Tseluiko, Baxter and Thiele27].

To correctly treat the singular point at $\eta =0$ , we formulate the system as two sets of ODEs, one starting at $\eta = -\Delta \eta $ , just to the left of the singular point and extending to $\eta = -L$ in the negative direction, and one starting from $\Delta \eta $ and extending to ${\eta = R}$ in the positive direction. Here, $\Delta \eta \ll 1$ and $R,L \gg 1$ ; in our calculations, we take ${\Delta \eta = 10^{-4}}$ , $R=8$ , and $L=5$ ; the exponentially decaying term in (2.5) means that the far-field condition is approached rapidly in each direction, so L and R do not have to be very large. Defining the scaled flux variable $h(\eta ) = u(\eta ) v'(\eta )$ , we then have the system

$$ \begin{align*} \frac{{d} f^\pm}{{d} \eta} &= \frac{f^\pm - g^\pm}{\eta}, \\ \frac{{d} g^\pm}{{d} \eta} &= \frac{h^\pm}{f^\pm}, \\ \frac{{d} h^\pm}{{d} \eta} &= -\eta \frac{h^\pm}{f^\pm}, \end{align*} $$

where $\pm $ superscripts refer to the positive and negative $\eta $ -directions, respectively.

Close to the singular point, we use the expansion (2.5) as boundary conditions:

$$ \begin{align*} \left. \begin{array}{ll} f^\pm = a_0 \pm a_1 \Delta\eta \log(\Delta\eta) \pm a_2^\pm \Delta\eta , \\ g^\pm = a_0 \pm a_1\Delta\eta , \\ h^\pm = a_0a_1, \end{array}\!\! \right\} \quad \eta = \pm\Delta\eta. \end{align*} $$

Here, $a_0$ , $a_1$ and $a_2^\pm $ are parameters that will be varied to satisfy the far-field conditions. In the far field, we apply

$$ \begin{align*} f^+ - g^+ &= \mathcal E^+, \quad f^+ = C_0^+, \quad \eta = R, \\ f^- - g^- &= \mathcal E^-, \quad f^- = C_0^-, \quad \eta = -L. \end{align*} $$

Here, $\mathcal E^{\pm }$ represents error terms from the desired boundary conditions (2.3) (in which f and g are equal), while $C_0^\pm $ represents each far-field value. The aim of the homotopy continuation procedure will be to reduce each of $\mathcal E^\pm $ to zero, set $C_0^- = 1$ and $C_0^+ = \epsilon $ for any given $\epsilon $ .

We proceed as follows. We start with the exact linear solution to (2.2), $f = \eta /2$ , ${g = 1}$ , with domain size parameters $R = L = 1$ (the parameter $\Delta \eta = 10^{-4}$ is held constant throughout the continuation procedure). This initial solution then satisfies the imposed boundary conditions, given the parameter values:

$$ \begin{align*} a_0 = 1, \quad a_1 = 0, \quad a_2^\pm = \mp \frac{1}{2}, \quad C_0^- = \frac{3}{2}, \quad C_0^+ = \frac{1}{2}, \quad \mathcal E^\pm = \mp \frac{1}{2}. \end{align*} $$

We then apply numerical continuation, first reducing the left-hand error $\mathcal E^-$ to zero and the left-hand value $C_0^-$ to one, while letting $a_0, a_2^-$ and $\mathcal E^+$ vary. We then reduce the right-hand error $\mathcal E^+$ to zero and $C_0^+$ to $\epsilon =0.1$ , letting each of $a_0, a_1, a_2^\pm $ vary. Finally, we increase the domain size L, R until the domain size has negligible effect on the solution. At the end of this process, we have a similarity solution for the given $\epsilon $ . We may then numerically continue in $\epsilon $ to find the solution for any such $\epsilon $ .

In Figure 2(c,d), we plot the resulting similarity solution for $\epsilon = 0.1$ . Comparing this similarity solution to the PDE solution profiles collapsed onto the self-similar coordinates, we see the profiles agree exactly; thus, confirming the similarity solution is an exact solution of the original PDE for the given initial condition. For the remainder of the study, we thus focus purely on the calculation of the similarity solution itself.

2.3 Numerical results for small $\epsilon $

Now, having a method of calculating similarity solutions, we can readily examine how solutions depend on the height $\epsilon $ . In Figure 3, we show how the constants $a_0$ and $a_1$ depend on $\epsilon $ over the range $\epsilon \in (0,1)$ . These constants may be considered characteristic parameters describing a solution to the system (1.5), given

$$ \begin{align*} u(x^*_0,t) = v(x^*_0,t) = a_0, \quad \frac{\partial v}{\partial x}(x^*_0,t) = -\frac{a_1}{t}, \end{align*} $$

where $x^*_0$ is the location of the initial discontinuity.

Figure 3

Coefficients from the expansion (2.5) of similarity solutions (2.2) over front height $\epsilon $ : (a) the value $a_0$ ; (b) $a_1$ and (c) $a_2^\pm $ . The parameters $a_0$ and $a_1$ show a turning point, below which $a_0$ and $a_1$ appear to slowly return to $1$ and $0$ , respectively. The value of $a_2^+$ increases in magnitude (slowly) as $\epsilon \to 0$ . The crosses marked on the curves in panels (a,b) correspond to the values $\epsilon = 10^{-1}, 10^{-2}, 10^{-3}$ shown in Figure 4(a,b), while the asterisks correspond to the values $\epsilon = 0.268, 0.0078$ shown in Figure 4(c,d).

As $\epsilon $ decreases, at first $a_0$ decreases and $a_1$ increases. However, there is a turning point at approximately $\epsilon = 0.1$ below which $a_0$ begins to increase and $a_1$ to decrease. Focusing on small $\epsilon $ , in Figure 4(a,b), we plot solution profiles for $\epsilon = 10^{-1}, 10^{-2}, 10^{-3}$ . For small $\epsilon $ , the similarity solutions are qualitatively different on either side of the singular point $\eta =0$ . In particular, the solution profiles are quite flat to the left of the singular point, but to the right, a steep front forms. Close to the bottom of the front, the solution reaches a point where the slopes of both f and g rapidly flatten as f and g approach $\epsilon $ . However, the values of $a_2^\pm $ differ (see Figure 3(c)), in particular with $a_2^+ < 0$ increasing in magnitude as $\epsilon $ becomes small. This contrast indicates the magnitude of derivatives of a solution to the left and the right of the singular point become markedly different.

Figure 4

Similarity profiles (a) f and (b) g as $\epsilon $ becomes small. An apparent interface at $\eta =\eta ^*$ is seen, where the magnitude of f becomes of order $\epsilon $ and rapidly flattens out. In panels (c,d), two solutions for different values may have the same value of $a_0$ and, thus, the same profile for $\eta < 0$ .

The most notable feature for solutions when $\epsilon $ is small is the formation of what could be considered an apparent interface at a location $\eta ^*> 0$ . This is not a true interface, as for $\epsilon $ positive, $v\to \epsilon $ never vanishes, but the slope of g appears almost vertical near this point, and the change in slope of f and g occurs very rapidly. In a solution to the PDE (1.5), this apparent interface will occur at $x^*(t) = x_0^* + \eta ^*t$ ; thus, the value of $\eta ^*$ in the similarity solution represents the speed of propagation of such a point.

To characterize the location of $\eta ^*$ over $\epsilon $ , we define it as the point at which the slope of g becomes smallest (that is, most negative). In Figure 5, we plot this quantity, along with the difference between $a_0$ and its asymptotic value $1$ , and the slope $a_1$ . As each of these quantities is seen to tend to zero as $\epsilon \to 0$ , but slowly, we plot these against the logarithmic variable $[\log (1/\epsilon )]^{-1/2}$ . Here, we see initial signs that these variables are of this logarithmic order as $\epsilon \to 0$ , which will motivate the asymptotic analysis in the next section (for comparison, we also plot the results of the asymptotic analysis that we will undertake). Note that in this plot, we take $\epsilon $ down to $10^{-6}$ , which is at approximately the limit of the numerical solution procedure. For this value, $[\log (1/\epsilon )]^{-1/2} \approx 0.27$ is still relatively large, which increases the difficulty of deducing the behaviour in this limit.

Figure 5

Apparent interface location $\eta ^*$ , as well as $1-a_0$ (where $a_0 = f(0) = g(0)$ ) and slope ${a_1 = -g'(0)}$ , in the limit where $\epsilon $ becomes small. Numerical results indicate that these solution parameters are all going to zero, although only logarithmically fast, namely, of order $[\log (1/\epsilon )]^{-1/2}$ as $\epsilon \to 0$ . Also plotted are the two-term asymptotic approximations (3.12) for each parameter (note these approximations are valid in the limit $\epsilon \to 0$ ).

Finally, we note some interesting implications of the nonmonotonic nature of the constants $a_0$ and $a_1$ seen in Figure 3. One implication, which can be seen in the small $\epsilon $ solution in Figure 4(a,b), is that a smaller value of $\epsilon $ to the right of the singular point may result in a larger value of g (and thus v) to the left of the singular point. This means that given initial conditions for which the diffusivity u vanishes, the problem (1.5) does not necessarily satisfy a comparison principle; one may have two solutions $v_1$ and $v_2$ , where $v_1 \leq v_2$ at $t=0$ , but $v_1$ does not remain bounded by $v_2$ for $t> 0$ .

The profile to the left of $\eta =0$ is uniquely characterized by $a_0$ , as for a given $a_0$ , both $a_1$ and $a_2^-$ are determined by requiring the $\eta \to -\infty $ condition to be met. A second implication of nonmonotonic $a_0$ is thus that two different values of $\epsilon $ can result in two similarity solutions that are identical for $\eta <0$ , while being distinct for $\eta> 0$ . In Figure 4(c, d), we plot the similarity solutions for $\epsilon \approx 0.268$ and $\epsilon \approx 0.0078$ , both corresponding to a value of $a_0 = 0.7$ . Here, we see the exact agreement for $\eta < 0$ . Given the similarity solutions represent exact solutions of the PDE, this implies that two initial conditions to (1.5) which are the same on one side of the singular point $x = x^*_0$ could remain the same for $t>0$ , even if this is not the case on the other side of the singular point.

3.1 Negative $\eta $

When $\eta < 0$ , the profiles f and g are both constrained in magnitude by the far-field condition $f, g \to 1$ as $\eta \to -\infty $ . Given the definition (3.1), we assume the asymptotic expansion

$$ \begin{align*} f(\eta) = 1 + \delta f_1(\eta) + \delta^2 f_2(\eta) + \cdots, \quad g(\eta) = 1 + \delta g_1(\eta) + \delta^2 g_2(\eta) + \cdots. \end{align*} $$

Substituting into the ODE system (2.2) and taking successive terms in $\delta $ , we have

$$ \begin{align*} \frac{{d} ^2g_1}{{d} \eta^2} + \eta \frac{{d} g_1}{{d} \eta} = 0, \quad \frac{{d} ^2g_2}{{d} \eta^2} + \eta \frac{{d} g_2}{{d} \eta} = -\frac{{d} }{{d} \eta}\bigg(f_1\frac{{d} g_1}{{d} \eta}\bigg), \end{align*} $$

and so on, with each $f_j$ found from $g_j$ by integrating

$$ \begin{align*} f_j = -\eta \int \frac{g_j}{\eta^2} \, d\eta, \quad j = 1, 2, \ldots. \end{align*} $$

To satisfy the far-field condition, we require $f_j, g_j\to 0$ as $\eta \to -\infty $ , while from the definition of $\delta $ ,

$$ \begin{align*} f_j(0), g_j(0) = \begin{cases} -1, & j = 1, \\ 0, & j = 2,3,\ldots .\end{cases} \end{align*} $$

We thus have

(3.2a) $$ \begin{align} g_1 = -\bigg(1 + \operatorname{\mathrm{erf}}\bigg(\frac{\eta}{\sqrt 2}\bigg)\bigg), \quad f_1 = \frac{\eta}{\sqrt{2\pi}}\operatorname{\mathrm{Ei}}\bigg(-\frac{\eta^2}2\bigg) - 1 - \operatorname{\mathrm{erf}}\bigg(\frac{\eta}{\sqrt{2}}\bigg), \end{align} $$

where $\operatorname {\mathrm {erf}}$ is the error function and $\operatorname {\mathrm {Ei}}$ is the exponential integral, defined for negative arguments by

$$ \begin{align*} \operatorname{\mathrm{Ei}}(s) = \int_{-\infty}^s \frac{{e}^{s'}}{s'} \, ds'. \end{align*} $$

The next term $g_2$ is complicated but is possible to find explicitly by integrating

$$ \begin{align*} g_2 &= \bigg(1+\operatorname{\mathrm{erf}}\bigg(\frac{\eta}{\sqrt{2}}\bigg)\bigg) \bigg(-\frac{e^{-{\eta^2}/{2}} \eta}{\sqrt{2 \pi }}-\frac{1}{3}+\frac{1}{3\pi }-\frac{\log (2)}{3 \pi }\bigg) -\frac{1}{3}\operatorname{\mathrm{erf}}\bigg(\frac{\eta}{\sqrt{2}}\bigg)^2 \notag \\ &\qquad+\frac{e^{-{\eta^2}/{2}} (\eta^2-1)}{3 \pi }\operatorname{\mathrm{Ei}}\bigg(-\frac{\eta^2}{2}\bigg) +\frac{\operatorname{\mathrm{Ei}}(-\eta^2)}{3 \pi }-\frac{e^{-\eta^2}}{3 \pi }+\frac{1}{3}. \end{align*} $$

The derivative of $g_2$ at $\eta =0$ is needed to give the correction to the slope $a_1$ that matches with that for $\eta>0$ . Differentiating the above, we find that

$$ \begin{align*} a_1 = a_{11}\delta + a_{12}\delta^2 + \cdots, \end{align*} $$

where

(3.3) $$ \begin{align} a_{11} = -\frac{{d} g_1}{{d} \eta}\bigg|_{\eta\to0^-} = \sqrt{\frac{2}{\pi}}, \quad a_{12} = -\frac{{d} g_2}{{d} \eta}\bigg|_{\eta\to0^-} = \frac{5\pi - 2 - 2\log(2)}{3\sqrt{2}\pi^{3/2}}. \end{align} $$

3.2 Positive $\eta $ away from apparent interface

We now consider the interval ${0 < \eta < \eta ^*}$ , where $\eta ^*$ is the apparent interface location. Since $\eta ^*$ is seen to decrease at the same logarithmic rate as $\delta $ , we introduce the scaled independent variable $\xi $ by $\eta = \delta \xi $ , so that (2.2) becomes

(3.4) $$ \begin{align} f - \xi\frac{{d} f}{{d} \xi} = g, \quad \frac{{d} }{{d} \xi}\bigg(f\frac{{d} g}{{d} \xi}\bigg) = -\delta^2\xi\frac{{d} g}{{d} \xi}. \end{align} $$

We take the apparent interface location $\eta ^* = \delta \xi ^*$ , where $\xi ^*$ has a series expansion in $\delta $ :

$$ \begin{align*} \xi^* = \xi_0^* + \delta \xi_1^* + \cdots. \end{align*} $$

Expanding the similarity profiles f and g in this region as

$$ \begin{align*} f = f_0(\xi) + \delta f_1(\xi) + \delta^2 f_2(\xi) + \cdots, \quad g = g_0(\xi) + \delta g_1(\xi) + \delta^2 g_2(\xi) + \delta^3 g_3(\xi) + \cdots, \end{align*} $$

the first terms must satisfy the differential equations

$$ \begin{align*} &f_0 - \xi \frac{{d} f_0}{{d} \xi} = g_0, \quad \frac{{d} }{{d} \xi}\bigg(f_0\frac{{d} g_0}{{d} \xi}\bigg) = 0, \\ &f_1 - \xi \frac{{d} f_1}{{d} \xi} = g_1, \quad \frac{{d} }{{d} \xi}\bigg(f_0\frac{{d} g_1}{{d} \xi} + f_1\frac{{d} g_0}{{d} \xi}\bigg) = 0. \end{align*} $$

At $\xi =0$ , the value and derivative of g must match with those for $\eta \to 0^-$ , giving the boundary conditions

$$ \begin{align*} g_0 = 1, \quad g_1= -1, \quad \frac{{d} g_0}{{d} \xi}= 0, \quad \frac{{d} g_1}{{d} \xi} = 0, \quad \xi = 0. \end{align*} $$

Another boundary condition comes from the apparent interface at $\xi ^*$ . We expect f to be exponentially small in $\delta $ close to this point, so expanding $\xi ^*$ and f, each term in $\delta $ must vanish:

$$ \begin{align*} f_0(\xi_0^*) =0, \quad f_1(\xi_0^*) + \frac{{d} f_0}{{d} \xi}\bigg|_{\xi=\xi^*}\xi_1^* = 0. \end{align*} $$

This set of equations and boundary conditions implies that $f_0$ and $f_1$ are linear, while $g_0$ and $g_1$ are constant:

(3.5) $$ \begin{align} g_0 = 1, \quad g_1 = -1, \quad f_0 = 1 - \frac{\xi}{\xi_0^*}, \quad f_1 = -1 + \frac{\xi_0^* + \xi_1^*}{\xi_0^{*2}}\xi. \end{align} $$

We can continue to higher terms. We will calculate the next corrections $g_2$ and $g_3$ , which satisfy

$$ \begin{align*} \frac{{d} }{{d} \xi}\bigg(f_0\frac{{d} g_2}{{d} \xi}\bigg) = 0, \quad \frac{{d} }{{d} \xi}\bigg(f_0\frac{{d} g_3}{{d} \xi}\bigg) = -\frac{{d} }{{d} \xi}\bigg(f_1\frac{{d} g_2}{{d} \xi}\bigg) \end{align*} $$

with boundary conditions

$$ \begin{align*} g_2 = g_3 = 0, \quad \frac{{d} g_2}{{d} \xi} =9 -a_{11}, \quad \frac{{d} g_3}{{d} \xi} = -a_{12}, \quad \xi=0. \end{align*} $$

Note that the higher corrections to f are not required to calculate g to this order. Integrating and applying the conditions at zero, we find

(3.6a) $$ \begin{align} g_2 &= a_{11} \xi_0^* \log\bigg(1- \frac{\xi}{\xi_0^*}\bigg) , \end{align} $$

(3.6b) $$ \begin{align} g_3 &= a_{11}\frac{\xi_1^*\xi}{\xi_0^*-\xi} + (a_{12}\xi_0^* + a_{11}\xi_1^*)\log\bigg(1-\frac{\xi}{\xi_0^*}\bigg). \end{align} $$

3.3 Inner problem near apparent interface

The expansions (3.5) and (3.6) are singular in the limit $\xi \to \xi ^*$ . There is thus a region close to $\xi =\xi ^*$ where the solution will match onto an inner solution, in which the far-field condition $f \to \epsilon $ is attained. We thus define inner variables z, $\hat f$ and $\hat g$ by

$$ \begin{align*} \xi = \xi^* + \delta^{-2}\epsilon z, \quad f = \epsilon \hat f(z), \quad g = \delta^2 \hat g(z). \end{align*} $$

Here, the scale of f is determined so that $\hat f \to 1$ as $z\to \infty $ , while the scales for z and g are determined by finding the appropriate balance of terms in (3.4). Under this change of variables, the system (3.4) is, to leading order,

(3.7) $$ \begin{align} -\xi^* \frac{{d} \hat f}{{d} z} = \hat g, \quad -\xi^*\frac{{d} \hat g}{{d} z} = \frac{{d} }{{d} z}\bigg(\hat f \frac{{d} \hat g}{{d} z}\bigg). \end{align} $$

Note that a result of solving the system in a small region near $\xi ^*$ is that the system is autonomous. Solutions of (3.7) are effectively travelling wave solutions of (1.5); this inner problem is thus relevant close to an apparent interface not just for the similarity solution, but in more general contexts (see further discussion in Section 4).

To match to the far-field condition, we require $\hat f \to 1$ , $\hat f' \to 0$ as $z\to \infty $ . Integrating (3.7) twice and applying these conditions, we find

(3.8) $$ \begin{align} \frac{{d} \hat f}{{d} z} = -\xi^*\log(\hat f); \end{align} $$

thus, we have an implicit solution as follows:

(3.9) $$ \begin{align} z = z_0 - \frac{1}{\xi^*}\operatorname{\mathrm{li}}(\hat f) = z_0 - \frac{1}{\xi^*}\operatorname{\mathrm{Ei}}(\log(\hat f)), \quad \hat g = \xi^{*2}\log(\hat f), \end{align} $$

where $z_0$ is an arbitrary translational constant and $\operatorname {\mathrm {li}}$ is the logarithmic integral. Note that $\hat g$ vanishes in the limit $z \to \infty $ , which is consistent with the definition of $\hat g$ given $\delta ^2 \gg \epsilon $ ; an additional inner problem would be required to match g onto its far-field condition $g \to \epsilon $ .

What remains is to find the asymptotic relationship between $\delta $ and $\epsilon $ , as well as the asymptotic formula for the apparent interface location $\xi ^*$ . To do so, we need to match the $z\to -\infty $ expansion of the inner problem (3.9) to the $\xi \to \xi ^*$ expansion of the outer problem (3.5), (3.6).

While the asymptotic expansion for the inverse logarithmic integral could be applied directly to (3.9), a simpler way to carry out the matching is to write the relationship between $\hat g$ and $\hat f$ in (3.9) in terms of the outer variables, giving

(3.10) $$ \begin{align} g = \xi^{*2}\delta^2\bigg[\log\bigg(\frac{1}{\epsilon}\bigg) + \log(f)\bigg] + O(\epsilon). \end{align} $$

To have an $O(1)$ term on the right of (3.10) to match that on the left, we require $\log \bigg (\frac {1}{\epsilon }\bigg ) = O(\delta ^{-2})$ . Thus, we let

$$ \begin{align*} \log\bigg(\frac{1}{\epsilon}\bigg) = \ell_0\delta^{-2} + \ell_1\delta^{-1} + \cdots, \quad \xi^* = \xi^*_0 + \delta \xi^*_1 + \cdots. \end{align*} $$

Substituting these expansions into (3.10), along with those for f and g, (3.5), (3.6), we have from each order:

(3.11a) $$ \begin{align} g_0 &= \xi_0^{*2} \ell_0, \end{align} $$

(3.11b) $$ \begin{align} g_1 &= \xi_0^{*2}\ell_1 + 2\xi_0^*\xi_1^* \ell_0, \end{align} $$

(3.11c) $$ \begin{align} g_2 &= \xi_0^{*2}\log(f_0) + \mathrm{constant}, \end{align} $$

(3.11d) $$ \begin{align} g_3 &= \xi_0^{*2}\frac{f_1}{f_0} + 2\xi_0^*\xi_1^* \log(f_0) + \mathrm{constant}. \end{align} $$

Here, the constant terms in the latter two equations involve higher order corrections ( $\xi ^*_2$ , $\ell _2$ and so forth).

The terms in the expansions for $\xi _0$ and $\log ({1}/{\epsilon })$ are now found as follows. The constants $\xi _0^*$ and $\xi _1^*$ are found by comparing the logarithmic terms in (3.11)(c,d) with those in (3.6). The more singular term in $g_3\, \propto f_0^{-1}$ then matches automatically. The constants $\ell _0$ and $\ell _1$ are then found by comparing the constant terms in (3.11)(a,b) to (3.5). We thus find

$$ \begin{align*} \xi_0^* = a_{11}, \quad \xi_1^* = a_{12}, \quad \ell_0 = \frac{1}{a_{11}^2}, \quad \ell_1 = -\bigg(\frac{1}{a_{11}^2} + \frac{2a_{12}}{a_{11}^3} \bigg), \end{align*} $$

where $a_{11}$ and $a_{12}$ are the terms in the expansion of the slope at $\eta =0$ (3.3).

Finally, since $\epsilon $ is the independent parameter of the original problem, we invert the series expansions to find the front at $\eta ^* = \delta \xi ^*$ . Thus,

(3.12a) $$ \begin{align} \delta \sim &\ell_0^{1/2}\log\bigg(\frac{1}{\epsilon}\bigg)^{-1/2} + \frac{\ell_1}{2}\log\bigg(\frac{1}{\epsilon}\bigg)^{-1} \notag \\ &\quad = \sqrt{\frac{\pi}{2}} \log\bigg(\frac{1}{\epsilon}\bigg)^{-1/2} - \frac{4\pi - 1 - \log 2}{6}\log\bigg(\frac{1}{\epsilon}\bigg)^{-1}, \end{align} $$

(3.12b) $$ \begin{align} \eta^* \sim &\xi_0^*\ell_0^{1/2} \log\bigg(\frac{1}{\epsilon}\bigg)^{-1/2} + \bigg(\frac{\xi_0^*\ell_1}{2} + \xi_1^*\ell_0\bigg) \log\bigg(\frac{1}{\epsilon}\bigg)^{-1} \notag\\ &\quad = \log\bigg(\frac{1}{\epsilon}\bigg)^{-1/2} - \frac{\sqrt\pi}{2\sqrt 2} \log\bigg(\frac{1}{\epsilon}\bigg)^{-1}, \end{align} $$

(3.12c) $$ \begin{align} a_1 \sim & a_{11}\ell_0^{1/2} \log\bigg(\frac{1}{\epsilon}\bigg)^{-1/2} + \bigg(\frac{a_{11}\ell_1}{2} + a_{12}\ell_0\bigg) \log\bigg(\frac{1}{\epsilon}\bigg)^{-1} \notag \\ &\quad = \log\bigg(\frac{1}{\epsilon}\bigg)^{-1/2} - \frac{\sqrt\pi}{2\sqrt 2} \log\bigg(\frac{1}{\epsilon}\bigg)^{-1}. \end{align} $$

Note that the expansions for the front location, $\eta ^*$ , and the slope at zero, $a_1$ , are the same to this order.

3.4 Comparison with numerical results

We of course want to compare the above asymptotic result with the numerical results of Section 2. We plot the expansions (3.12) in Figure 5 along with the numerical calculations of the same quantities. The logarithmic dependence of the solution parameters on $\epsilon $ is certainly consistent with their observed numerical behaviour, even if the agreement is not very close at the level that we can achieve numerically. A major challenge in comparing the asymptotic solution directly to numerical results is the slow convergence of the logarithmic terms. For $\epsilon = 10^{-6}$ , which is around the smallest value we can attain with our numerical method, $[\log ({1}/{\epsilon })]^{-1/2} \approx 0.27$ , which is not very small. At this value, the correction terms in the asymptotic formulae (3.12) do improve the leading order result. For $\epsilon = 10^{-6}$ , the leading order asymptotic estimates of $\delta $ in (3.12a) and $\eta ^*$ in (3.12b) are respectively

$$ \begin{align*} \delta \sim \sqrt{\frac{\pi}{2}}\bigg[\log\bigg(\frac{1}{\epsilon}\bigg)\bigg]^{-1/2} \approx 0.34, \quad \eta^* \sim \bigg[\log\bigg(\frac{1}{\epsilon}\bigg)\bigg]^{-1/2} \approx 0.27. \end{align*} $$

Including the $O([\log ({1}/{\epsilon })]^{-1})$ correction terms leads to estimates $\delta \approx 0.206$ and ${\eta ^* \approx 0.224}$ . By comparison, the values estimated from the numerical solution in Section 2.3 are

(3.13) $$ \begin{align} \delta = 1-a_0 \approx 0.227, \quad \eta^* \approx 0.241. \end{align} $$

However, by using the numerical estimates for $\delta $ and $\eta ^*$ in (3.13) instead of asymptotic approximations, we can find encouraging agreement with the solution profiles in each region. In Figure 7(a), we plot the asymptotic solution for $\eta < 0$ and $0 < \eta < \eta _0$ . Here, we use the $O(\delta )$ expansions (3.2a) for $\eta < 0$ , and for $\eta> 0$ , we use the approximations for f and g that satisfy the conditions at $\eta =0$ and include their first nonuniform terms, that is,

(3.14) $$ \begin{align} g \approx 1 - \delta + \delta^2\sqrt{\frac{2}{\pi}}\log\bigg(1-\frac{\eta}{\eta^*}\bigg), \quad f \approx (1-\delta)\bigg(1- \frac{\eta}{\eta^*}\bigg). \end{align} $$

Figure 7

(a) Comparison between the numerically computed similarity solution with $\epsilon =10^{-6}$ , and approximate profiles given by (3.2a) for $\eta < 0$ and (3.14). Due to the very slow convergence of the logarithmic terms as $\epsilon \to 0$ , numerical estimates of $\delta $ and $\eta ^*$ are used in fitting these profiles. (b) Asymptotic and numerical profiles in the inner region near $\eta ^*$ . Since the asymptotic solution is accurate to $O(\epsilon )$ in this region, the agreement is excellent.

In Figure 7(b), we compare the solution close to $\eta ^*$ , scaled into variables z, $\hat f$ and $\hat g$ , with the inner solution (3.9). To do so, we also have to numerically estimate the translational constant $z_0$ (this constant is exponentially small in $\delta $ , so in terms of the asymptotic analysis, cannot not be estimated from matching to any finite number of terms of the outer problem). Otherwise, the agreement is practically perfect, which is to be expected given the error in the inner region is of order $\epsilon $ , rather than $\delta $ .