Self-similar structure of the entry and exit discontinuities in two-way diffusion equations

Abstract

Two-way diffusion equations arising in kinetic problems relating to electron scattering and in Brownian particle dynamics present singularities absent from conventional diffusion equations. Although calculations by Stein & Bernstein, and Fisch & Kruskal have revealed the formation of entry and exit slope discontinuities at the critical points where the velocity changes sign, the analytical structure of these discontinuities remains unclear. Here we fill this gap via a local similarity variable analysis, illustrated through the two-way diffusion equation $y \partial n/\partial x=\partial ^2 n/\partial y^2$ in $-1 \leq y \leq 1$; $0 \leq x \leq L$, with $n(x,\pm 1)=0$ with various entry conditions $n(0,y)_{y\gt 0}$, and the exit condition $n(L,y)_{y\lt 0}=0$. The similarity variable $\eta =y/x^{1/3}$ permits the analytical characterization of the entry discontinuity, except for constants determined by matching with numerical solutions obtained with two numerical schemes: separation of variables following the construction of Beals, or finite-difference discretization of the transient partial differential equation, which converges in time to a solution almost identical to the separation of variables solution. Although the slope discontinuity depends markedly on the initial condition $n(0,y)_{y\gt 0}$, a simple general similarity solution structure emerges empirically, always involving a spontaneous singular contribution $C |y|^{1/2}$ at $x=0,y\lt 0$. Slow convergence of both numerical solutions near $\{x,y\}=\{0,0\}$ is attributed to the poor eigenfunction representation of the ever-present singular solution component $|y|^{1/2}$. The similarity approach applies equally to other two-way diffusion equations when the coefficient of $\partial n/\partial x$ changes sign linearly with $y$. It can also be extended to situations where this coefficient is discontinuous at the critical points.

1. Introduction

Two-way diffusion equations such as (1.1) arise for instance in the kinetic theory of electron scattering, when the coefficient $u(y)$ multiplying the time-like derivative term $\partial /\partial x$ changes sign at one or more values of the spatial variable $y$ :

(1.1) \begin{equation} u(y) \frac {\partial n}{\partial x} = \frac {\partial }{\partial y} \left [ D(y) \frac {\partial n}{\partial y} \right ]\! . \end{equation}

This sign inversion alters drastically the nature of the solution and calls for special methods of analysis. Indeed, a similar contrast exists between one-way and two-way diffusion equations as with fluid boundary layer problems without and with flow separation, since the separation arises at a point of velocity inversion that complicates drastically the nature of the problem.

Two-way diffusion equations also arise in the transport of ions and charged nanoparticles carried in air by combined electric and flow fields, and in the kinetic theory of Brownian particles. When analysed by separation of variables, two-way diffusion equations have the peculiarity of possessing a degenerate zero eigenvalue. Completeness of the eigenfunction basis then requires the addition of a special solution of the form $x-g(y)$ (Bethe, Rose & Smith 1938; Stein & Bernstein 1976; Fisch & Kruskal 1980). This so-called diffusion solution does not appear in the examples treated in the present study, but it is often of fundamental importance, irrespective of whether the problem is treated by separation of variables or by other means. One generally considers boundary conditions leading to a self-adjoint problem. The domains of positive and of negative $u(y)$ will be referred to, respectively, as the positive or forward and the negative, backward or retrograde regions. Initial conditions are stipulated partially at $x=0$ on the forward region, and partially at the downstream end $x=L$ on the backward region:

(1.2) \begin{align} & n(0,y)=p^+(y) \quad \mbox{where} \quad u(y)\gt 0, \end{align}

(1.3) \begin{align} & n(L,y)=p^-(y) \quad \mbox{where} \quad u(y)\lt 0. \end{align}

The general problem and several concrete examples have been considered in early studies by Bethe et al. (1938), Stein & Bernstein (1976) and Fisch & Kruskal (1980). The latter contributed substantially to the treatment by separation of variables, and also conjectured that the set of eigenfunctions with negative (positive) eigenvalues was complete in the region of positive (negative) $u(y)$ . Beals (1981) proved this conjecture and a number of important additional theorems. He also contributed a method to construct the solution based on the complete set of positive and negative eigenfunctions and proved its existence. A substantial literature cites these developments.

The present investigation focuses on the solution behaviour near the point(s) $y_c$ of velocity inversion, $u(y_c)=0$ , which we shall denote critical point(s). These points are evidently singular enough to change the nature of the differential equation, calling for tools of analysis departing considerably from the case with no such criticality. An additional driver of discontinuous behaviour may be tentatively expected at critical points, as the initial conditions affecting them on their immediate right and left come from distant boundaries. A discontinuity may also be tentatively surmised from Beals’ construction of the solution, where the initial and final conditions, applicable, respectively, in the positive and the negative domain, are unified into a single function $p(y)$ defined over the whole $y$ domain. This function is, from its construction, typically discontinuous at the critical points. Since $p(y)$ is the main driver of the solution to the linear and otherwise homogeneous problem, it is conceivable that a discontinuity in $p(y)$ might be transmitted to the solution $n(x,y)$ .

Although not evident in prior calculations on the electron scattering problem (Bethe et al. 1938), a discontinuity in the $y$ derivative of the ion concentration at the entry and exit sections has been clearly demonstrated for the same problem by Stein & Bernstein (1976), who developed a remarkably accurate finite difference scheme to deal with its specific difficulties. They also challenged their readers by noting that ‘no analytic results have been forthcoming with regard to this singularity’. This theoretical gap still remains true today. Fisch & Kruskal (1980) treated numerically a different problem, selected ad hoc to produce simple eigenfunctions, and also encountered a discontinuity in the slope of the entry and exit profiles. Even though the discontinuous $u(y)$ used by Fisch & Kruskal (1980) differed drastically from the linear variation of $u(y)$ at the critical point in Stein & Bernstein (1976), the two discontinuities found numerically looked rather similar, suggesting a certain universality. In both calculations the slope jumped from zero to a value seen graphically to be large, possibly infinity. Similar slope discontinuities have been subsequently encountered and described by various numerical schemes (Kim & Tranquilli 2008). Other expressions for the coefficient $u(y)$ have been investigated, including singular ones, though no related information on singularity formation has been forthcoming (Fleige 2023). These general considerations and specific precedents on the probable theoretical interest of the solutions to (1.1) near entry and exit critical points provide sufficient stimulus to take the challenge of Stein and Bernstein to determine the analytical structure of these regions.

As additional motivation, we have observed such singularities in our numerical analysis of two-way diffusion problems connected to experimental studies on the motion of a gas carrying ions inside a tube. The gas-driven ion motion is opposed by a uniform axial electric field, such that the combination of both fields creates regions of retrograde ion motion near the tube walls. Other flow fields with velocity inversions supporting similar two-way transport problems are not uncommon. A prototypical example in the absence of electric fields is found at the boundary of the two fluids in the Kelvin Helmholtz instability problem (Batchelor 1967, p. 511). Velocity inversions are readily established in electrokinetic problems in liquids, where one has control of the direction of the electric field relative to a pressure driven flow field. For example, Rubinstein & Zaltzman (2013) have illustrated not only the generation of such flow fields, but also the associated two-way transport problem. These flows in liquids have the potential to create rather rich and complicated structures, as they may generate circulation, vortices and intense mixing, even in porous media with strong viscous dissipation (Mirzadeh et al. 2020). The case considered here of ions moving in a gas requires a wall boundary condition of zero ion concentration, which excludes the diffusion solution. A zero-flux boundary condition may be more appropriate in electrokinetic problems in liquids, in which case the diffusion solution could be active. This does not remove the entry singularities, as illustrated by the work of Stein & Bernstein (1976) and Fisch & Kruskal (1980). As a final point, we note that the whole field of gas dynamics is fundamentally linked to a two-way diffusion problem: the Boltzmann equation, whose degenerate zero eigenvalue is associated with so-called normal or diffusion solutions. These are nothing less than the conservation equations of fluid dynamics and heat and mass transport. Bidirectional diffusion is an inescapable feature of kinetic problems, because the coefficient of the spatial gradient is an independent velocity variable that takes both positive and negative values at every point in space.

We shall later describe two numerical solution schemes to solve (1.1) for the particular case when $u(y)=y$ and $D(y)=1$ , illustrating the imperfect representations given by both approaches near the singular regions. These limitations suggest the need for a local treatment carried out presently via similarity solutions that describe analytically the singular region. Though the present detailed analysis is restricted to cases where $u(y)$ varies linearly with $y$ near the critical point, it may be extended to the problem of Fisch & Kruskal (1980), where $u(y)$ is a step function. In addition to the rich range of singular structures here found, the casuistic encountered in various examples using different initial conditions sheds useful light into the role of the full and partial range solutions to the problem, their connection to the issue of existence, and the reason for the poor convergence of numerical solutions near entry and exit critical points.

The focus of this study is the local structure of the singularity, analysed in § 5. However, the local problem involves unknown constants, whose determination requires the global solution in the vicinity of the singularity. Two numerical methods used extensively here to compute the global field are summarily described in § 2 and § 3 and in detail in Appendix A.

2. Numerical analysis by separation of variables

In this section we follow earlier approaches to solve (1.1) by separation of variables via (2.1), which leads to the eigenvalue problem (2.2):

(2.1) \begin{align} & n = e^{-\lambda x} F(y), \end{align}

(2.2) \begin{align} & u(y) \lambda F(y) + [D(y)F^{\prime}(y)]^{\prime} = 0, \end{align}

with suitable self-adjoint boundary conditions at two fixed $y$ boundaries of the domain.

Although our notation may appear to differ from that of Beals (1981), it merely adapts his projection operators to our numerical approach. We distinguish between the positive eigenvalues $\lambda _j$ and corresponding eigenfunctions $F_{\!j}(y)$ relevant in the forward region, and their negative counterparts, $\mu _{\!j}$ and $G_{\!j}(y)$ relevant in the retrograde region. The ion concentration receives contributions from the negative eigenfunctions at the downstream end, and from the positive eigenfunctions at the upstream end. These penetrate by diffusion from one domain into the other, so that the number concentration at any position is the linear combination of positive and negative eigenfunctions,

(2.3) \begin{equation} n(x,y)= a_j F_{\!j}(y) \exp (-\lambda _j x) + b_{\!j} G_{\!j}(y) \exp [-\mu _{\!j} (x-L)], \end{equation}

where summation over repeated indices is implicit. Convergence of the series in the mean, rather than uniform convergence, is understood throughout the paper.

The partial initial and final conditions (1.2)–(1.3) then imply

(2.4) \begin{align} & p^+(y)= a_j F_{\!j}(y) + b_{\!j} G_{\!j}(y) \exp (\mu _{\!j} L) \quad \mbox{where} \quad u(y)\gt 0,\end{align}

(2.5) \begin{align} & p^-(y)= a_j F_{\!j}(y) \exp (-\lambda _j L) + b_{\!j} G_{\!j}(y) \quad \mbox{where} \quad u(y)\lt 0 . \end{align}

Following Beals (1981), these two statements, each valid in each of the two regions of the $y$ domain, may be combined into a single statement valid over the whole domain by defining a single function $p(y)$ taking the values $p^+(y)$ and $p^-(y)$ in the forward and backward regions, respectively. Equations (2.4) and (2.5) hence give a single entry/exit condition statement similarly as in the single initial condition for a proper Sturm–Liouville problem. The remaining task is just to determine the $a_j$ and $b_{\!j}$ coefficients in (2.3) for given $p(y)$ . Bethe et al. (1938) gave approximate solutions by retaining only the lowest modes. Stein & Bernstein (1976) discarded the eigenfunction series approach for a number of practical reasons such as: (i) the inconvenience of calculating many non-standard eigenfunctions; (ii) the unavailability of orthogonality properties for the half-range problem and the resulting non-diagonal matrices. Fisch & Kruskal (1980) followed the separation of variables method while obviating the first difficulty by selecting a step function for $u(y)$ , resulting in simple eigenfunctions. They determined the expansion coefficients by minimized square errors. Their seemingly arbitrary discontinuous $u(y)$ is in fact physically relevant, as it arises in different two-way diffusion equations modelling counter-current separators (Fitt, Ockendon & Shillor 1985; Hagan & Ockendon 1991).

Today, most of the practical difficulties noted in these pioneering studies have lost much of their sting. For example, the computer program Mathematica calculates fast and accurately a stunning number of known functions. It can also rapidly and precisely calculate and store hundreds of eigenfunctions for problems such as that of Bethe et al., and of Stein and Bernstein, whose eigenfunctions are not standard functions. Furthermore, the determination and inversion of relatively large non-diagonal matrices is no longer particularly demanding with the approach to be described in (A25) and (A26). The lack of orthogonality properties can sometimes be obviated by complex variable methods (Hagan & Losek 1999). In one instance, the velocity distribution $f(u,x)$ of heavy molecules in a light carrier gas undergoing a shock wave was approximately determined while skipping the task of inverting non-diagonal matrices, and apparently without involving a slope discontinuity at the critical point $u=0$ (Fernández-Feria & Fernández de la Mora 1987).

The full details of this calculation method are given in Appendix A.

3. Numerical analysis by time-dependent finite differences

The second numerical treatment of the problem uses a standard finite difference discretization of the time-dependent form of (1.1). A second-order explicit upwind discretization of the $x$ -derivative is combined with an implicit centred discretization of the $y$ -derivatives, and the discretized problem is marched in time, with a time step limited by a Courant–Friedrichs–Lewy condition, until the solution becomes stationary. Time marching for this linear problem amounts to using an iterative scheme to invert the system matrix. Non-uniform $x$ - and $y$ -grids are used to increase the resolution around critical points at the entry and exit sections while keeping moderate grid sizes. A typical grid for $u=y$ , $D=1$ , with $L=0.1$ , has $5000 \times 300$ points and gives minimal $x$ - and $y$ -steps of $1.65 \times 10^{-7}$ and $6.70 \times 10^{-5}$ , respectively, around $x=y=0$ .

4. Numerical results for Couette flow

4.1. The discontinuity

We restrict our analytical and numerical investigations to the following problem:

(4.1) \begin{equation} y \frac {\partial n}{\partial x} = \frac {\partial ^2 n}{\partial y^2} \quad \mbox{with} \quad n(x,\pm 1) = 0, \end{equation}

with certain generalizations sketched occasionally. As a first example we select the following initial and final conditions:

(4.2) \begin{gather} n(0,y)_{y\gt 0}=1, \end{gather}

(4.3) \begin{gather} n(L,y)_{y\lt 0}=0 \quad \mbox{with} \quad L=0.1 . \end{gather}

The same two-way diffusion equation (4.1), with the same boundary condition $n(x,\pm 1)=0$ (Kaper et al. 1984; Keller & Weinberger 1997), or the alternative condition $\partial n/\partial y=0$ at $y=\pm 1$ (Rubinstein & Zaltzman 2013), have been previously considered, though without focus on entry and exit singularities.

Figure 1. ( $a$ – $c$ ) Comparison between the profiles $n(x,y)$ computed by either time evolution of the finite difference equations (dashed red), or the eigenfunction expansion with 112 modes (black) for Couette flow with $L=0.1$ . Here ( $a$ ) $x=\{0, 0.03, 0.06, 0.1\}$ , top to bottom. ( $b$ ) Detail of the entry region with $x = 10^{-5} \{0, 1, 2, 5, 10\}$ . ( $c$ ) Detail of the exit region with $x=\{0.09, 0.095, 0.099, 0.0999, 0.1\}$ . Both solutions are indistinguishable in the figures almost everywhere. However, the finite difference method suggests the existence of a discontinuity in $n(x,y)$ at $y=0$ at the entry and exit points ( $x=0$ and $x=L$ ), not apparent in the eigenfunction series. ( $d$ ) Comparison of the full global solution (A11)–(A12) at $x=0$ (black) for $L=\{0.3, 0.1, 0.05\}$ (top to bottom) with the half-range solution (4.5) (red-dashed) extended to negative $y$ , showing their close coincidence for $L \geq 0.3$ . The half-range solution exhibits the same discontinuous slope as the full-range solution.

Details of the calculation are given in Appendix A. Figure 1 shows numerical results for the finite domain $0 \leq x \leq L$ ; $L=0.1$ . The solutions obtained by finite differences (dashed red) or through an eigenfunction sum with 112 terms (black) are visually indistinguishable almost everywhere. However, the finite difference calculation suggests the existence of a discontinuity at $y=0$ at the entry and exit points ( $x=0$ and $x=L$ ). The truncated eigenfunction expansion matches the other calculation for $x \gt 2 \times 10^{-5}$ , but does not capture the apparent discontinuity. Its inadequacy in the singular region of very small $x$ is confirmed in later calculations increasing the number of modes from 112 to 500.

4.2. The source of the slope discontinuity is not a discontinuous initial condition $p(y)$

Figure 2 shows several other cases computed (with 500 eigenfunctions of each sign) with the same null final condition $n(L,y)_{y\lt 0}=0$ , and the following three alternative initial conditions in $y\gt 0$ and $L=0.1$ :

(4.4) \begin{equation} n(0,y)_{y\gt 0} = [4y(1-y)]^{s+1}, \quad s=\{0, 3, 24\}. \end{equation}

These initial conditions are selected such as to increase steadily the number of $y$ derivatives $\partial ^s n(0,y)/\partial y^s$ that vanish at $y=0^+$ . As can be seen, the entry discontinuity to the left of $y=0$ persists even with $s=24$ . Accordingly, a discontinuity in the solution still arises when the joint initial/final condition function $p(y)$ is highly smooth at $y=0$ . This behaviour will later be confirmed in the case of initial value functions with all derivatives null at $y=0$ , hence with a completely smooth $p(y)$ . The source of the singularity is evidently in the differential equation, rather than in a discontinuous initial/final condition. This conclusion is confirmed in the limit $L \rightarrow \infty$ , as the discontinuity persists when the two boundaries have no influence on each other. The generality of this induction is reinforced by the calculation of Fisch & Kruskal (1980) for a different problem. Their Gaussian initial condition decaying rapidly to zero in the singular point at the entry (hence matching smoothly with the null outlet condition) still gave a prominent entry discontinuity.

Figure 2. Full-range initial ( $x=0$ , black) and final ( $x=L$ , black dashed) profiles for $n(L,y)_{y\lt 0}=0$ ; $n(0,y)_{y\gt 0}=[4y(1-y)]^{s+1}$ (red, dashed), computed with 500 eigenfunctions of either sign for $L=0.1$ . The entry discontinuity in $[\partial n(0,y)/\partial y]_{y=0}$ persists even when the initial and final conditions $p^{\pm }(y)$ merge at $y=0$ with a high level of continuity. $s=\{0, 3, 24\}$ in ( $a$ ), ( $b$ ) and ( $c$ ).

4.3. Relation between the half-range and the full-range solutions

The notions of partial and full range solutions to two-way diffusion equations has played an important role in the proof of existence (Fisch & Kruskal 1980; Beals 1981). Partial completeness implies that any function $p^+(y)$ used as initial condition, $n(0,y)=p^+(y)$ in the forward domain may be represented in terms of the eigenmodes with negative exponents (positive domain). Hence,

(4.5) \begin{equation} p^+(y) = \sum a_j^{\prime} F_{\!j}(y) \quad \mbox{where} \quad u(y)\gt 0 . \end{equation}

Projecting (4.5) on the $F_n$ functions over the positive domain leads to the algebraic system (in matrix form, see Appendix A)

(4.6) \begin{equation} \boldsymbol{P}^{\prime}= \boldsymbol{a}^{\prime} \boldsymbol{\cdot }\boldsymbol{\mathsf{F}} \quad \mbox{or} \quad \boldsymbol{a}^{\prime}=\boldsymbol{P}^{\prime} \boldsymbol{\cdot }\boldsymbol{\mathsf{F}}^{-1}, \end{equation}

with

(4.7) \begin{equation} P_k^{\prime} = \int _+{p^+(y) F_k(y) u(y) \, \textrm {d} y} \quad \mbox{and} \quad F_{\!jk}=\int _0^1 {F_{\!j}(y) F_k(y) u(y) \, \textrm {d} y} . \end{equation}

The half-range solution is accordingly

(4.8) \begin{equation} n^{\prime}(x,y) = \sum a_j^{\prime} F_{\!j}(y) \exp (-\lambda _j x) . \end{equation}

Figure 1(d) includes the 500 eigenmode representation of this half-range solution for the Couette problem (4.1) with $n(0,y)_{y\gt 0}=1$ (dashed red line) determined from (4.6) and extended to the negative domain. The figure shows also the corresponding representation of the full-range problem with $p^-(y)=0$ (black lines). The half-range eigenfunction sum matches well the initial condition $n=1$ at $y\gt 0$ , as expected from the completeness of the positive eigenfunctions for $y\gt 0$ . The figure represents ion concentration profiles at $x=0$ only, and includes full-range calculations for three values of $L$ : 0.3, 0.1, 0.05 (black, top to bottom). One can see that the full-range solution has reached its large $L$ asymptote at $L=0.3$ . Interestingly, this asymptote agrees closely with the half-range solution extended into $y\lt 0$ , which evidently does not depend on $L$ . This appears to confirm Fisch and Kruskal’s (1980) notion of a close relation between the large $L$ asymptote of the full solution and the half-range solution. However, the two are not identical, as the large $L$ limit of the full problem given by (A15)–(A16) in Appendix A differs from (4.6). In fact, the full- and half-range solutions depart substantially from each other for modest values of $L$ . The reasons for the approximate coincidence found at large $L$ between $n$ and $n^{\prime}$ extended into the negative region are given in Appendix A. There we also argue that, in the real situation of large $L$ , the full-range solution may be viewed as a small perturbation of the half-range solution. Therefore, the invertibility of the half-range matrix $\boldsymbol{\mathsf{F}}$ would imply that of the corresponding matrix of the full-range problem. Consequently, the existence of a solution to the half-range problem would imply the existence of a solution to the full-range problem. These numerical calculations thus illustrate one important element of Beals’ general proof of existence (Beals 1981).

Also of considerable interest in figure 1(d) is that the $x \rightarrow 0$ limit of the half-range solution (4.8) extended into the negative domain results in very much the same derivative discontinuity as the full solution. In other words, the matrix $\boldsymbol{\mathsf{F}}$ associated with the positive eigenfunctions embodies the source of the singularity. For all the initial conditions considered in § 5, we have confirmed numerically that the half-range solution exhibits a slope singularity very similar to that in the full solution computed for $L=0.1$ (Appendix C).

It may seem paradoxical that series (4.5) involving the $F_k(y)$ functions that decay exponentially towards negative $y$ would produce the algebraic decay observed at $y\lt 0$ . Nevertheless, the first few $F_k(y)$ modes penetrate sufficiently into the negative domain to enable this behaviour.

We note finally a method introduced by Dita (1986), where the full-range solution may apparently be constructed from the half-range solution. At the end of his article, Dita invokes the Picard–Lindelof theorem allowing the continuation of the half-range solution into the other half-range.

5. Local similarity analysis of the structure of the discontinuity

5.1. Similarity formulation, eigenfunctions and their properties

Neither the discontinuity at $y=0$ , nor its evolution with $x$ near the entry and the exit are easily describable by global numerical analysis of the full domain. Here we use the fact that the governing equation (4.1) admits solutions in terms of the similarity variable

(5.1) \begin{equation} \eta = y/x^{1/3} \end{equation}

to seek solutions in the form $n(x,\eta)$ . The boundary conditions given in terms of $y$ make the problem non-separable in the $\eta$ variable. However, as long as $x$ is small, $\eta$ will be large at the boundary $y=1$ . We therefore expect that a similarity solution with boundaries at infinity will be adequate, at least at small $x$ . In terms of the new independent variables $x^{\prime}=x$ , $\eta =y/x^{1/3}$ , (4.1) becomes (after dropping the prime in $x^{\prime}$ )

(5.2) \begin{equation} \eta x \frac {\partial n}{\partial x} = \frac {\eta ^2}{3} \frac {\partial n}{\partial \eta } + \frac {\partial ^2 n}{\partial \eta ^2} . \end{equation}

This equation may be analysed by separation of variables,

(5.3) \begin{equation} n(x,y) = x^m H_m(\eta), \end{equation}

where $H_m$ satisfies

(5.4) \begin{equation} -m \eta H_m + \frac {\eta ^2}{3} H_m^{\prime} + H_m^{\prime\prime} = 0 . \end{equation}

The eigenvalue $m$ cannot be negative, as the solution is proportional to $x^m$ and would not be finite at $x=0$ . The option $m=0$ leads to the familiar Leveque solution $H_0^{\prime}(\eta)=C \exp (-\eta ^3/9)$ , which must also be discarded because it diverges exponentially for $\eta \lt 0$ . Other than for a constant solution also associated with the eigenvalue $m=0$ , $m$ must take positive values.

The lack of a complete set of eigenfunctions in the positive domain makes this eigenvalue problem rather different from that encountered in the original $(x, y)$ variables, suggesting that the conditions for existence of a similarity solution must be newly developed. The obvious difference between the $(x,\eta)$ and the $(x,y)$ eigenvalue problems is that the $x$ dependence of the eigenmodes has shifted from exponential to a power law. The two would be equivalent in terms of the new variable $s = \ln x$ , though this change would make the problem singular at $x=0$ . The equivalent of the diffusion solution would in this case have the also singular form $\ln x + g(\eta)$ .

The asymptotic behaviour of the solutions to (5.4) at large $|\eta |$ can be inferred by assuming it is a power law. The two first terms in (5.4) would then be of the order $\eta ^3$ larger than the third and be dominant. This results in

(5.5) \begin{equation} H_m \rightarrow C_m^{\pm } |\eta |^{3 m} \quad \mbox{for} \quad |\eta | \gg 1, \end{equation}

where the two constants $C_m^{\pm }$ in the regions of positive and negative $\eta$ may evidently differ from each other. To the next approximation one would include as a non-homogeneous term the neglected $H_m^{\prime\prime}$ term computed from (5.5). Repeating the process leads to the following asymptotic form:

(5.6) \begin{align} & H_m(\eta) \rightarrow C_m^{\pm } |\eta |^{3 m} Q_m(\eta) \quad \mbox{with} \\ & Q_m(\eta) = \sum \frac {a_k}{\eta ^{3 k}} = \left [ 1 - 9 m \left ( \frac {1-3 m}{3 \eta ^3} + \frac {4 - 19 m + 34 m^2 - 9 m^3}{2 \eta ^6} \right) + \ldots \right ]\!, \nonumber \end{align}

where we have chosen $a_0=1$ , such that $Q_m$ asymptotes to $+1$ at large values of $|\eta |$ , and the sign of $C_m^{\pm }$ coincides with the sign of $H_m(\eta)$ at both the positive and the negative asymptotes. Additional terms in series (5.6) may be obtained through the recurrence relations

(5.7) \begin{equation} a_0=1, \quad a_k=a_{k-1} \frac {3(m-k+1)(3m-3k+2)}{k} . \end{equation}

A power series with such rapidly growing coefficients ( $a_k/a_{k-1} \rightarrow 9 k$ ) has zero radius of convergence. Accordingly (5.6) is at best an asymptotic series restricted to large values of $|\eta |$ , where only a finite number of terms should be kept. This series is nevertheless most useful, even though the solutions to (5.4) may be written in terms of Kummer’s confluent hypergeometric function $K_1$ :

(5.8) \begin{equation} H_m(\eta) = D_1 \eta K_1 \left ( \frac {1}{3}-m, \frac {4}{3}, -\frac {\eta ^3}{9} \right) + D_2 K_1 \left ( -m, \frac {2}{3}, -\frac {\eta ^3}{9} \right)\! . \end{equation}

Both independent solutions in (5.8) diverge as $\exp (-\eta ^3 /9)$ at large negative $\eta$ . However, the following linear combination cancels the exponential divergence, resulting in the desired solution with power law asymptotes at both large positive and negative $\eta$ :

(5.9) \begin{equation} H_m(\eta) = \eta K_1 \left ( \frac {1}{3}-m, \frac {4}{3}, -\frac {\eta ^3}{9} \right) + K_1 \left ( -m, \frac {2}{3}, -\frac {\eta ^3}{9} \right) \, \frac {3^{2/3} \varGamma (4/3) \varGamma (-m)}{\varGamma (2/3) \varGamma (1/3-m)} . \end{equation}

The numerical representation of $H_m(\eta)$ via (5.9) fails at large negative $\eta$ because it involves the difference between two astronomically large quantities. The asymptotic description (5.6) is accordingly more useful in that domain. An alternative accurate method to determine $H_m(\eta)$ is to integrate numerically (5.4) starting at a large negative $\eta$ with an initial condition given by the asymptotic expansion (5.6). The integration in the opposite direction is unstable in the negative domain, as it inevitably picks up the mode diverging as $\exp (-\eta ^3/9)$ .

Figure 3. Representation of the $H_m(\eta)$ functions for the values of $m$ indicated in the legend (bottom to top).

Figure 3 shows the $H_m(\eta)$ functions defined in (5.8) for various values of $m$ . The graphs extend down to $\eta =-7$ , the smallest $\eta$ for which Mathematica reliably calculates $H_m(\eta)$ . Of particular interest in figure 3 is the eigenfunction corresponding to $m=1/6$ , which asymptotes to zero at large positive $\eta$ . Except for a unit vertical displacement and the multiplication by a constant, this eigenfunction looks very much like what is seen in figure 1(b).

Also of relevance in cases when the initial condition $n(0,y)_{y\gt 0}$ is an analytic function of $y$ are the eigenfunctions $H_{k/3}(\eta)$ with $k$ a positive integer. Two out of three of these functions, characterized by indexes $k$ or $k+1/3$ , are polynomials. These polynomials are special cases where the coefficients in the asymptotic expansion (5.6) vanish beyond a certain power of $y^{-3}$ . Note that solution (5.8) in terms of Kummer’s function fails when $m$ is an integer. In such cases the asymptotic series (5.6) gives an exact solution when extended up to the point where $k=m$ , since all subsequent $a_k$ values are exactly null according to (5.7). A few of the polynomial eigenfunctions are given in table 1.

Table 1. Selected normalized $h_m(\eta)$ eigenfunctions (defined in (5.10)–(5.11)) given by polynomials.

5.1.1. Determination of $C_m^{\pm }$

The $C_m^{\pm }$ coefficients in (5.5) play a key role in the following analysis. They are defined as

(5.10) \begin{equation} C_m^{\pm } = \lim _{z\rightarrow \pm \infty } \frac {H_m(z)}{|z|^{3 m}} . \end{equation}

Here $C_m^-$ is particularly relevant as it is used to normalize the eigenfunctions $H_m$ via

(5.11) \begin{equation} h_m(\eta) = \frac {H_m(\eta)}{C_m^-} . \end{equation}

As a result of this normalization, the terms relevant for the calculation become

(5.12) \begin{equation} c_m^{\pm } = \frac {C_m^{\pm }}{C_m^-} = \lim _{z \rightarrow \pm \infty } \frac {h_m(z)}{|z|^{3 m}} . \end{equation}

An analytic expression giving the $m$ dependence of the $c_m^+$ coefficients will be given below as (5.16). Here $c_m^- =1$ for all $m$ values as a result of the normalization of $h_m$ .

The determination of $C_m^-$ requires the accurate computation of the ratio $H_m(\eta)/|\eta |^{3 m}$ as $\eta \rightarrow -\infty$ , which cannot be done by Mathematica based on Kummer’s function (5.9). Accordingly, we exploit the fact that the asymptotic description retaining $m$ or more terms provides an accurate value of $H_m(\eta)$ for $|\eta |\gt 6$ , such that the ratio $H_m(\eta)/Q_m(\eta)$ is indistinguishable from a constant for $|\eta |\gt 6$ . The choice $\eta =6$ is relatively arbitrary. It is a compromise of an $\eta$ large enough for the asymptotic result (5.6) for $Q_m$ to be accurate, and small enough for Mathematica’s numerical calculation of $H_m(\eta)$ to also be unproblematic. We therefore calculate this constant ratio at $\eta =-6$ . The rescaled function $H_m(\eta) Q_m(-6)/H_m(-6)$ agrees precisely with $Q_m(\eta)$ when $\eta \lt -6$ . Therefore, we compute $C_m^-$ as

(5.13) \begin{equation} C_m^- = \frac {H_m(-6)}{6^{3 m} Q(-6)}, \end{equation}

where we have used the fact that, by its definition, $Q_m(\infty)=1$ . This provides the following operational definition of the normalized eigenfunctions:

(5.14) \begin{equation} h_m(\eta) = \frac {H_m(\eta) \, 6^{3 m} Q_m(-6)}{H_m(-6)} . \end{equation}

The same process repeated at positive $\eta$ yields

(5.15) \begin{equation} C_m^+ = \frac {H_m(6)}{6^{3 m} Q_m(6)} . \end{equation}

With this calculation method we compute $c_m^+$ as a function of $m$ , and discover that it is an oscillating function numerically indistinguishable from the following analytical function:

(5.16) \begin{equation} c_m^+ =-2 \sin [\pi (m-1/6)] . \end{equation}

Note that the application of these results to increasing values of $m$ requires extending the asymptotic expansion (5.6) of $Q_m$ at least to order $k=m$ .

A noteworthy property of the $x^m H_m(\eta)$ functions defining the solution in (5.3) is that their large $|\eta |$ (small $x$ ) asymptotes depends only on $y$ as $C_m^{\pm } |y|^{3 m}$ .

The special property that $C_{1/6}^+=0$ seen in figure 3 may appear to be a coincidence. However, this particular $m$ may also be viewed as an eigenvalue satisfying the boundary condition $H_m(+\infty)=0$ . In light of the periodic dependence (5.16) of $c_m^+$ with $m$ , a whole sequence of such eigenvalues is $m_k=k+1/6$ , $k$ being a positive or negative integer. For the positive $m_k$ sequence, $H_m$ undergoes an increasing number of oscillations in the region $\eta \lt 0$ .

5.2. The relevance of eigenvalues $m \ne 1/6$

We have seen that the asymptotic behaviour as $x \rightarrow 0$ of the mode $x_m H(m,\eta)$ is $C_m^{\pm } |y|^{3 m}$ . Except for the special values of $m$ when these functions are polynomials (table 1), their limit as $x \rightarrow 0$ exhibits a slope discontinuity at $y=0$ when $0\lt m\lt 1/3$ . This may be seen in figure 4, with the smallest slope jump shown at $m=0.3$ (bottom curve), and the largest one at $m=0.05$ (top curve). For $m\gt 1/3$ discontinuities still arise but only on the second or higher $y$ derivatives. We shall see that all these discontinuities (except for $m=\{1/6, 7/6, \ldots \}$ ) must be excited through an initial condition at $y\gt 0$ with contributions of the form $y^{3 m}$ . We shall also see that the discontinuities associated with the $m=\{1/6, 7/6, \ldots \}$ modes do (apparently) always arise spontaneously, without being actively excited through an initial condition of the form $y^{3 m}$ . This point has already been seen with the initial condition $n(0,y)_{y\gt 0}=1$ . The step in the slope encountered when the initial condition at $y\gt 0$ is a constant is accordingly not the only possible discontinuity, nor the most perverse.

Figure 4. Depiction of the limit as $x \rightarrow 0$ of the eigenfunctions $x^m H_m(y/x^{1/3}) \rightarrow C_m^{\pm } |y|^{3 m}$ , with $m$ indicated in the legend (bottom to top). All profiles have been normalized at $y=-1$ . The ratio $C^-/C^+$ is calculable via (5.16).

5.3. Structure of the similarity solution

Through a series of examples with initial conditions of increasing complexity we discover a simple structure of the similarity solution, to be conjectured here with greatest generality, and to be empirically illustrated in the various examples in § 5.4 and in Appendix D.

1. (i) For sufficiently small $x$ , and for the full range $-1\lt y\lt 1$ , the similarity solution has two contributions. The first, referred to as homogeneous, has the form $\sum t_j x^{j+1/6} h(j+1/6,y/x^{1/3})$ , with $j=0,1,2,\ldots$ , with coefficients $t_j$ not determinable via local considerations.

2. (ii) For initial conditions of the form $n(0,y)_{y\gt 0}= \sum r_k y^{3 m(k)}$ , where the exponents $m(k)$ do not need to be integers, the second contribution to the similarity solution has the form $\sum r_k^{m(k)} h[m(k),y/x^{1/3}]/c_{m(k)}^+$ , with $c_m^+=-2 \sin [\pi (m-1/6)]$ . The coefficients $t_k$ are determinable from the computed global solution $n(x,y)$ via the statement

   (5.17) \begin{equation} \sum t_i |y|^{3 i} = \frac {n(0,y)_{y\lt 0} - \sum r_k |y|^{3 m(k)}/c_{m(k)}^+}{|y|^{1/2}} . \end{equation}
   There is apparently assurance that the right-hand side in this equality may be represented by a series in powers of $|y|^3$ , fixing uniquely the $t_k$ coefficients. This series often includes only a handful of terms.

3. (iii) For an initial condition $n(0,y)_{y\gt 0}$ with an essential singularity such that all its derivatives at $y=0^+$ vanish, the $t_i$ coefficients are determinable as in (ii) after taking all the $r_k$ to be null.

4. (iv) In cases when the initial condition is piecewise continuous, the relevant region fixing the $r_k$ coefficients in (5.17) is the one immediately to the right of the critical point.

5. (v) It therefore appears that the $r_k$ coefficients relevant to define the $t_k$ terms in (5.17) are determined by the continuation into the negative region of the initial condition on the right region by the rules of change from right to left of the $x^m h_m(y/x^{1/3})$ eigenmodes as $x \rightarrow 0$ .

5.4. Effect of the initial condition $n(0,y)_{y\gt 0}$ on the similarity solution

All the calculations in this section are for $L=0.1$ and use 500 eigenfunctions.

5.4.1. $n(0,y)_{y\gt 0}=1$

This step function case has already been considered in figure 1. A striking observation is the several coincidences between figure 3 and the $m=1/6$ curve in figure 4. Their asymptote at negative $y$ depends only on $y$ , and their asymptote at $y\gt 0$ is null. This situation suggests the possibility that, except for the addition of a constant, the entry ( $x=0$ ) discontinuity seen in figures 1 and 2 will have the following form involving certain constant coefficients $t_k$ :

(5.18) \begin{equation} n(x,y) = \sum t_k x^{k+1/6} h(k+1/6,y/x^{1/3}) \end{equation}

with the limit

(5.19) \begin{equation} \lim _{x \rightarrow 0} n(x,y) = \sum t_k |y|^{3k+1/2} . \end{equation}

We have confirmed for the 500 term eigenfunction sum for $n(0,y)_{y\gt 0}=1$ ; $n(L,y)_{y\lt 0}=0$ , $L=0.1$ , that its $y\lt 0$ domain at $x=0$ is precisely represented by (5.19) as $n(0,y)_{y\lt 0}=1 - 1.149 |y|^{1/2} + 0.229 |y|^{7/2} - 0.08 |y|^{13/2})$ for the whole negative interval $-1\lt y\lt 0$ (figure 5 a). Furthermore, representation (5.18) $n(x,y)= 1 - 1.149 x^{1/6} h(1/6,y/x^{1/3}) + 0.229 x^{7/6} h(7/6,y/x^{1/3}) - 0.08 x^{13/6} h(13/6,y/x^{1/3})$ describes with notable precision the downstream evolution of the solution up to $x=0.01$ , except for the decay down to zero near the right boundary $y=1$ (figure 5 b). This sharp decay may be described locally in terms of the familiar similarity variable $(1-y)/x^{1/2}$ , and will not be pursued here.

Figure 5. ( $a$ , $b$ ) Comparison between the 500-term eigenmode sum with $L=0.1$ (dashed) and the similarity solution (continuous black lines). ( $a$ ) Full range comparison keeping three modes of the similarity solution such as to closely fit $n(0,y)_{y\lt 0}$ in the full negative range as $1 - 11.149 |y|^{1/2} + 0.229 |y|^{7/2} - 0.08 |y|^{13/2}$ . ( $b$ ) Magnified depiction of the singular region revealing slight inaccuracies in the truncated sum at $x=10^{-7}$ . ( $c$ , $d$ ) Transient finite difference solution. Panel ( $c$ ) shows the collapse of 24 different profiles at different $x$ values in terms of the similarity variable $\eta$ . Panel ( $d$ ) extrapolates the data in ( $c$ ) to $x=1.65 \times 10^{-9}$ to describe the corner region.

Figure 5(b) focuses on the small $y$ region of figure 5(a). It shows that the 500 term eigenfunction series fits the similarity solution down to $x=3 \times 10^{-7}$ , with a slight disagreement at $x=3 \times 10^{-7}$ . Although much better that the 112 term truncated sum, the 500 term sum is still unable to capture the singularity. It is in fact far from doing so, as the similarity solution shows that the discontinuous $x \rightarrow 0$ asymptote is not sharply displayed even at $x=10^{-11}$ . Figure 5(c) represents 24 different profiles computed with the most refined transient finite difference mesh described in § 3 in the range $1.34 \times 10^{-7} \leq x \leq 6.52 \times 10^{-6}$ . When plotted in terms of $\eta$ , they are essentially indistinguishable. Figure 5(d) extrapolates the results of figure 5(c) to $x=1.65 \times 10^{-9}$ to describe the corner region, which follows closely the $y^{1/2}$ prediction of the similarity law.

In conclusion, in the case of a uniform entry profile $n(0,y)$ we confirm the numerical findings by Stein & Bernstein (1976) and Fisch & Kruskal (1980) that the entry slope is discontinuous at the critical point. In addition, we find analytically that the initial singularity at $y\lt 0$ is proportional to $|y|^{1/2}$ . The slope is not just discontinuous, but jumps from 0 at $y=0^+$ to infinity at $y=0^-$ . This singular slope is what gave the appearance of a discontinuity in one of our numerical descriptions of this region (dashed line in figure 1 b).

5.4.2. $n(0,y)_{y\gt 0}=4y(1-y^2)$

For simplicity, our next choice of initial conditions is such that the two eigenfunctions degenerating into $n(0,y)_{y\gt 0}$ are polynomials, and therefore present no discontinuity at $x \rightarrow 0$ :

(5.20) \begin{align} 4 x^{1/3} h_{1/3}(y/x^{1/3}) - 4 x h_1(y/x^{1/3})=4y-4(y^3+6x) . \end{align}

The situation is in this respect analogous to that in the prior example with a constant initial condition corresponding to the mode $m=0$ ( $x^0 h_0(y/x^{1/3})=1$ ). The choice $n(0,y)_{y\gt 0}=4y(1-y^2)$ has the additional advantage of satisfying the boundary condition at $y=1$ , avoiding the steep decay from $n=1$ to zero seen for $x\gt 0$ at the right-hand boundary (dashed lines in figure 5 a).

Figure 6. Concentration profiles for initial and final conditions $n(0,y)_{y\gt 0}=y-y^3$ , $n(L,y)_{y\lt 0}=0$ and $L=0.1$ . The red dashed lines represent the 500 eigenmode sum for $x=\{10^{-3}, 3 \times 10^{-4}, 10^{-4}, 3 \times 10^{-5}, 10^{-5}, 3 \times 10^{-6}, 10^{-6}, 3 \times 10^{-7}, 10^{-7}\}$ . The black curves are the similarity solution, with the lowest included line corresponding to $x=10^{-14}$ . The global calculation is excellent down to $x=10^{-7}$ , while the similarity solution remains excellent even up to $x=10^{-3}$ .

The results of the global calculation keeping 500 eigenmodes of each sign are shown by dashed lines (red in the web edition) in figure 6. Like in the previous example with $n(0,y)_{y\gt 0}=1$ , an initial slope discontinuity arises at $y=0$ . And, like in the previous case, the initial profile in the negative region combines the modes, $m=1/6, 7/6, \ldots$ , with those injected through the initial condition

(5.21) \begin{equation} n(0,y)_{y\lt 0}=4 y - 4(y^3 + 6 x) + \sum t_k x^{k+1/6} h(k+1/6,y/x^{1/3}) . \end{equation}

In view of the fact that the singular modes with $m=\{1/6, 7/6, \ldots \}$ appear to arise spontaneously under all conditions, without being actively inserted through the initial and final conditions, we shall subsequently refer to them as the homogeneous modes.

By fitting the numerically obtained $n(0,y)_{y\lt 0} - 4 y (1-y^2)$ to the expected dependence $\sum t_k |y|^{3k+1/2}$ we find a precise matching with the choice $t_k=\{2.849, -2.93, 0.0850, -0.004, \ldots \}$ , with higher-order terms making negligible contributions. This fit to the global concentration profile at $x=0$ defines the similarity solution in the full domain as

(5.22) \begin{gather} n(x,y) = 4(y-y^3-6x) + 2.849 x^{1/6} h_{1/6}(y/x^{1/3}) - 2.93 x^{7/6} h_{7/6}(y/x^{1/3}) \\ + 0.085 x^{13/6} h_{13/6}(y/x^{1/3}) - 0.004 x^{19/6} h_{19/6}(y/x^{1/3}) . \nonumber \end{gather}

The agreement between the global (red dashed) and the similarity solution (black) seen in figure 6 is as good at all the $x$ values shown as at $x=0$ . One can see that the similarity solution loses precision at $x\gt 10^{-3}$ while the 500th term eigenfunction sum is slightly off at $x=10^{-7}$ and becomes inadequate at even smaller $x$ . A profound change in $n(x,y)$ is seen in the singular region in response to slight $x$ variations. In contrast, hardly anything changes elsewhere, up to the largest $x$ shown of $10^{-3}$ . The explanation for these trends lies in the fact that the first mode of the similarity solution involves the factor $x^{1/6}$ , taking values of order unity even at rather small values of $x$ . In contrast the eigenmode $y$ is strictly independent of $x$ , while the mode $y^3 + 6x$ changes little in the interval $0\lt x\lt 10^{-3}$ . The higher homogeneous modes evolve far more slowly with $x$ , at best as $x^{7/6}$ , explaining the slight $x$ dependence seen also to the left of the critical point.

5.4.3. General analytical initial condition

We now consider the case of an analytical initial condition $n(0,y)_{y\gt 0}=\sum r_k y^k$ , including powers $k=2, 5, 8$ , etc., whose corresponding eigenfunctions $x^{k/3} h_{k/3}(y/x^{1/3})$ present discontinuities as $x \rightarrow 0$ . The same procedure previously used may be followed to determine the $t_k$ coefficients, provided one accounts for the fact that the $x \rightarrow 0$ asymptotes in the positive and the negative $y$ domains differ from each other by a factor $c_m^+$ . In other words, the similarity solution must contain terms of the form

(5.23) \begin{equation} n(x,y)_{y\gt 0} = \sum r_k x^{k/3} h(k/3,y/x^{1/3})/c_{k/3}^+, \end{equation}

whose different $x \rightarrow 0$ limits are $n(0,y)_{y\gt 0} = \sum r_k y$ and $n(0,y)_{y\lt 0} = \sum r_k |y|^k/c_{k/3}$ . Accordingly, the unknown $t_k$ coefficients are now inferred by fitting the computationally obtained $n(0,y)_{y\lt 0} - \sum r_k |y|^k/c_{k/3}^+$ to the homogeneous pattern $\sum t_k |y|^{3k+1/2}$ . In other words, a more general rule for the structure of the similarity solution in the whole domain corresponding to initial condition $n(0,y)_{y\gt 0} = \sum r_k y^k$ with $k$ integer is

(5.24) \begin{equation} n(x,y) = \sum r_kx^{k/3} h(k/3,y/x^{1/3})/c_{k/3}^+ + \sum t_k x^{k+1/6} h(k+1/6,y/x^{1/3}) . \end{equation}

Note that (5.21) was the special case of (5.24) when the $c_m^+$ coefficients were all unity.

Although the rigorous conditions required for the validity of (5.24) remain to be determined, based on the examples examined, it is reasonable to hypothesize that the general form of the similarity solution in the whole $-1\lt y\lt 1$ domain is (5.24), with the coefficients $r_k$ fixed by the initial condition $n(x,y)_{y\gt 0} = \sum r_k y^k$ , the $c_m^+$ terms given by (5.16), and the homogeneous terms $t_k$ determinable by fitting to the global solution.

5.4.4. Non-analytic initial condition: $n(x,y)_{y\gt 0}=y^{3/5}-y^3$ and $y^{3/10}-y^3$

We shall now see that the solution structure (5.24) is preserved when the initial condition includes non-integer powers $y^{3 m}$ . $y^{3 m}$ is the $x \rightarrow 0$ limit of the eigenmode $x^m h_m(y/x^{1/3})$ , and therefore introduces the sharp entry slope discontinuities seen in figure 4 at $y=0$ in the range $0\lt m\lt 1/3$ . These discontinuities are nevertheless conceptually identical to the less spectacular discontinuities affecting higher-order derivatives with integer powers $2, 5, \ldots$ , and may be dealt equally through the $c_m^+$ coefficients.

Figure 7. Concentration profiles $n(x,y)$ for $n(0,y)_{y\gt 0}=y^{3/5}-y^3$ ( $a$ ), and $y^{3/10}-y^3$ (b). $L=0.1$ . The black and red dashed lines in ( $a$ ) and ( $b$ ) are, respectively, the 500 mode eigenfunction sum and the similarity law. Panel ( $c$ ) is a magnified version of ( $b$ ) with inverted colours.

Figure 7 shows ion concentration profiles $n(x,y)$ for $n(0,y)_{y\gt 0}=y^{3/5}-y^3$ (figure 7 $a$ ), and $y^{3/10}-y^3$ (figure 7 $b$ ), computed by retaining 500 eigenmodes with $L=0.1$ (red dashed lines). The $t_k$ coefficients defining the similarity solution (black) are obtained via (5.24) as $t_k=\{4.513, 0.747, 0.\}$ and $\{-2.667,0.7499, 0.01203\}$ . They yield an excellent match to the global calculation up to $x=10^{-3}$ . The rapidity with which these two rather spiky initial discontinuities relax at increasing $x$ is such that they would be almost indiscoverable in a global numerical calculation. For instance, figure 7(c) shows that the initial spike forming at $x=0$ has more than halved at $x=10^{-7}$ . Again, this stunningly fast relaxation is understandable by the small power $x^{1/10}$ governing the amplitude of the most singular mode appearing in this example.

5.4.5. Initial conditions with essential singularities and discontinuities

In Appendix D we analyse four additional cases of initial conditions. Two of them involve essential singularities, where all derivatives at $y=0^+$ are zero. Both have negative domain solutions where (5.17) is satisfied when all the $r_k$ terms vanish, as would be expected by analytical continuation of the initial condition into the negative domain. Two other examples involve initial conditions with discontinuities. (5.17) is satisfied in both when the $r_k$ coefficients are fixed by analytical continuation of the initial condition in the domain immediately to the right of the critical point. The rules formulated in § 5.3 for construction of the similarity solution are accordingly met in all the examples considered.

6. Slow convergence near the critical point

The poor convergence of the eigenfunction series observed in the vicinity of $y=0$ is attributable primarily to the fact that the main component of the singularity is the contribution degenerating into $|y|^{1/2}$ as $x \rightarrow 0$ . Under most common conditions, this term dominates as $y \rightarrow 0^-$ . Yet, it is poorly represented by a truncated series of negative eigenfunctions, all of which have zero second derivatives at $y=0$ , while the corresponding derivative of $y^{1/2}$ is highly singular. The situation is accordingly similar to the Gibbs phenomenon.

To illustrate numerically this point, we represent the function

(6.1) \begin{equation} f_0(y)= |y|^{1/2} \quad \mbox{for} \quad y\lt 0; \quad f_0(y)= 0 \quad \mbox{for} \quad y\gt 0, \end{equation}

in terms of the positive and negative eigenfunctions,

(6.2) \begin{equation} f_{0\sigma }(y) = \sum \left [ T_k F_k(y) + S_k G_k(y) \right ]\! . \end{equation}

A terminological distinction is made between the function and its eigenfunction sum representation, not only because of our truncated sum, but also because convergence is generally assured only in a mean sense. In this case $f_0(y)$ is given over the whole range, so we may take inner products of (6.2) with both eigenfunction families over the whole $y$ range to enjoy orthogonality properties, resulting in

(6.3) \begin{equation} T_k = \frac {\int _{-1}^0{ F_k(y) y |y|^{1/2} \, \textrm {d} y}}{f_{kk}} \quad \mbox{and} \quad S_k = \frac {\int _{-1}^0{ G_k(y) y |y|^{1/2} \, \textrm {d} y}}{-f_{kk}}, \end{equation}

where $f_{jk}$ is defined in (A8) and (A10).

Figure 8 compares the representation $f_{0\sigma }(y)$ (black) with the original function $f_0(y)$ (dashed) in the singular region, showing a level of disagreement in the corner region comparable to those observed for the full-range global solution. Figure 8(b) compares also the eigenfunction representation for the $x \rightarrow 0$ limit of the second homogeneous mode: $f_1(y)=|y|^{7/2}$ for $y\lt 0$ ; $f_1(y)=0$ for $y\gt 0$ , showing a close coincidence everywhere. The agreement is even better for the higher modes, confirming the intuitive notion that the poor convergence observed near the critical point is entirely due to the first homogeneous mode.

Figure 8. Panel ( $a$ ) shows the 500 eigenmode representation (6.2) of $f_0(y)$ (6.1) (black), as well as $f_0(y)$ itself (red dashed), displaying an excellent agreement except for $|y|\lt 0.03$ . Panel (b) includes also the eigenfunction representation of $f_1(y)$ (the $x \rightarrow 0$ limit of the next homogeneous eigenmode), demonstrating a close agreement everywhere between $f_1(y)$ and $f_{1\sigma }(y)$ . The poor convergence of the eigenfunction sum near $\{x,y\}=\{0,0\}$ is therefore attributable to the slow convergence of the eigenfunction sum description of $f_0(y)$ .

7. The effect of axial diffusion

We have seen very rapid evolutions of the entry concentration profile with $x$ , with drastic changes taking place in some instances in the narrow interval $0\lt x\lt 10^{-7}$ (figure 7c ). Under these circumstances, the hypothesis that the axial diffusion term $n_{xx}$ here ignored is in fact negligible is problematic. This situation is similar to that often arising in boundary layer problems in the vicinity of a sharp leading edge, as in the Blasius problem. In the similarity variables, the neglected $n_{xx}$ term makes a contribution equal to $x^{-4/3} \epsilon ^2$ times a certain function of $\eta$ , where $\epsilon$ is an inverse Peclet number. The similarity approach is thus expected to break down at small $x$ , and cannot be used to describe the region very close to the origin. Its predictions may perhaps be hoped to apply when $x$ is large compared with $\epsilon ^{3/2}$ . In this case, even in situations with relatively small $\epsilon \sim 10^{-3}$ , it would appear that the extreme spikiness of some of the entry concentration profiles predicted here would be smoothed by axial diffusion.

On the other hand, the similarity solutions found have the property that they depend only on $y$ at $x=0$ . Here $n_{xx}$ would then be exactly zero at the entry section, apparently implying that the similarity solutions obtained are self-consistent, at least in the large $\eta$ asymptote. This general property of all of our similarity solutions may be seen more concretely for the dominant singularity $n \sim x^{1/6}H_{1/6}(y/x^{1/3}$ ), for which

(7.1) \begin{equation} \frac {\epsilon ^2 n_{xx}}{n_{yy}}= \epsilon ^2 x^{-4/3} \frac { -5 H_{1/6}(\eta)/36 + \eta H_{1/6}^{\prime}(\eta)/3 + \eta ^2 H_{1/6}^{\prime\prime}(\eta)/9}{H_{1/6}^{\prime\prime}(\eta)} .\end{equation}

Interestingly, the $\eta$ dependence of this expression decays exponentially and becomes negligible for $\eta \lt -5$ , irrespective of the magnitude of $\epsilon ^2 x^{-4/3}$ . Accordingly, if this solution were to apply for $\eta \lt -5$ , the width in $y$ space of the rounded off region of the singularity would still scale as $\Delta y \sim x^{1/3}$ . In other words, the order of magnitude of the broadening of the spiky concentration profiles with $x$ would hardly be affected by axial diffusion. Of course, there is no assurance that the form of the solution will be preserved in the regions where $\eta \lt -5$ . A reliable determination of the broadening effect of axial diffusion should therefore rely on a direct attack to the full problem.

8. Conclusions

1. (i) We have implemented two numerical approaches to compute the global solution in two-way diffusion equations. One is a transient finite difference scheme. The other follows Beals’ solution construction method by separation of variables. When implemented in a Couette flow problem, both obtain the solution efficiently, except near the critical points at the entry and exit profiles. Both calculations hint at the presence of an entry discontinuity at the critical point, but fail to resolve it accurately.

2. (ii) For large enough $L$ the uncoupled solution to the half-range problem is seen to be numerically very close to the solution to the coupled full-range problem as a result of the relative smallness of the coupling matrix $\boldsymbol{\mathsf{N}}$ defined in Appendix A, such that $\boldsymbol{\mathsf{F}}^{-1} \boldsymbol{\cdot }\boldsymbol{\mathsf{N}} \ll \boldsymbol{\mathsf{I}}$ . This approximate match is explained by the fact that $\boldsymbol{\mathsf{N}}$ involves integrals whose integrands are the product of one oscillating function and another exponentially decreasing function. This applies both to integrals over the half- $y$ -range (when the integrals are small), and to the full $y$ -range (when the integrals vanish exactly).

3. (iii) Although large non-diagonal matrices need to be computed in the separation of variable method, the task is greatly simplified by a feature analogous to the orthogonality properties of proper Sturm–Liouville problems: we show in Appendix A that the needed half-range integrals defining the governing matrices are simply related to the value of the eigenfunctions and their derivatives at the critical point.

4. (iv) For the chosen Couette problem where the coefficient $u(y)$ of the $\partial n/\partial x$ term is strictly linear with $y$ , a local solution is obtained in terms of the similarity variable $\eta =y/x^{1/3}$ . This results in an eigenvalue problem with characteristics quite different from the eigenvalue problem encountered in the original variables $x$ , $y$ . In particular, there is no complete set of eigenfunctions for the positive domain, while an infinite number of complete sets of eigenfunctions exist in the negative domain.

5. (v) A particularly interesting similarity eigenfunction set degenerates in the entry section into the functions $y^{1/2+3k}$ , where $k$ is zero or a positive integer. These special eigenfunctions, referred to as homogeneous, decay to zero in the positive domain, and are the main cause of the entry slope discontinuity observed at the critical point. They apparently always contribute to the global solution through coefficients not determined by the similarity analysis, which we obtain by matching the expected form of the similarity solution to the numerically found global solution.

6. (vi) A general structure for the similarity solution is induced without proof from numerous examples. It includes the homogeneous terms with locally undetermined coefficients, as well as additional terms driven by the initial condition. The form of these additional terms is straightforward when the initial condition at $y\gt 0$ is given as the sum of integer or non-integer powers of $y$ . This situation offers a rich range of singularities at the critical points. In several other cases with initial conditions involving essential singularities or discontinuities the form of the additional terms is determined by analytical continuation of the initial condition given to the immediate right of the critical point.

7. (vii) While neither the numerical global nor the local similarity solution provide a full accurate description to the problem, the combination of both fixes the unknown coefficients in the similarity solution and provides an analytical description of the singular region.

8. (viii) Although the present analysis is restricted to the convective diffusive problem $y \partial n/\partial x=\partial ^2 n/\partial y^2$ , with $n(x,\pm 1)=0$ , the $|y|^{1/2}$ structure of the main singularity found is identical to that in the electron scattering problem analysed by Stein & Bernstein (1976) for another two-way diffusion equation with rather different boundary conditions. This local structure is in fact general for problems where the coefficient of $\partial n/\partial x$ varies linearly with $y$ near the critical point. The problem studied by Fisch & Kruskal (1980), where the coefficient of $\partial n/\partial x$ was discontinuous at the critical point, is also analysable by a different local similarity solution.

9. (ix) For brevity we have focused on the entry region, as the behaviour at the exit region is entirely analogous.

10. (x) The slow convergence of the eigenfunction representation of the solution near the critical point at the entry and exit is (under most conditions) due entirely to the first homogeneous self-similar eigenmode, whose $x \rightarrow 0$ limit $|y|^{1/2}$ is poorly represented by a truncated series of global eigenfunctions.

11. (xi) We have determined the similarity solution keeping enough terms $x^{k+1/6} h_{k+1/6}(y/x^{1/3})$ to describe the whole $y$ range at small enough $x$ . In most cases, however, keeping just the first term provides a sufficiently accurate representation of the singular region.