Analytic shock-fronted solutions to a reaction–diffusion equation with negative diffusivity

,


Introduction
Reaction-diffusion equations (RDEs) are partial differential equation (PDE) models that have countless applications, including in wound healing [35], biological invasion [24], and movement through porous media [30].The simplest reaction-diffusion models consider the evolution of one population in one spatial dimension.Mathematically, these models have the general form where u(x, t) is the density as a function of space, x, and time, t.In (1), the symbol D(u) is the diffusivity and the function R(u) is the reaction term capturing local changes in density.In standard linear Fickian diffusion, D(u) = D, where D > 0 is a constant.
However, many systems do not undergo linear diffusion, and in these scenarios, D(u) is non-constant.In general, we refer to non-constant D(u) as nonlinear diffusion.
Reaction-diffusion models with nonlinear diffusion commonly have D(u) ≥ 0 for all feasible values of u.Positive diffusivity results in the population dispersing from regions of high density to regions of low density.However, many natural processes involve aggregation, whereby members of the population move up density gradients from regions of low density to regions of high density.For example, a near-universal phenomenon in ecology is animal populations aggregating for social reasons including mating, prey detection, or predator avoidance [17].Typical models for aggregation involve members of the popula-tion responding to chemical stimuli [22].However, another way to model aggregation is to allow the diffusivity D(u) to become negative for some values of u.Negative diffusivity can arise in the continuum limit of discrete models with social aggregation [34,37].
In [34], the strength of inter-particle attraction influences the regions in which the diffusion is negative.In [37], by considering a random walk with attraction among individuals Turchin [37] derived a model with a diffusivity that could become negative.Johnston et al. [20] showed that an agent-based model containing both isolated and group-associated individuals which can have different rates of birth, death and movement gives rise to an RDE model with a region of negative diffusivity.In this work, we focus on models with D(a) = D(b) = 0, with D(u) < 0 for a < u < b, where a and b are constants such that 0 < a < b < 1, with D(u) ≥ 0 otherwise.Although it remains to be demonstrated that negative diffusivity can describe aggregation in a real system, negative diffusivity provides a way to generate aggregation using a single PDE, without requiring chemotaxis.
Allowing negative diffusivity D(u) causes ill-posedness in Cauchy initial-value problems involving the nonlinear diffusion equation u t = (D(u)u x ) x [1,2].It also causes solutions for u to have very steep gradients [37], and solutions can even become multi-valued.
Inserting a jump discontinuity (or shock) in a multi-valued solution enables a singlevalued solution to be obtained.A shock involves an instantaneous (with respect to x) jump from one density, u r , to another, u l , such that u never takes the values between u r ≥ b and u l ≤ a. Consequently, both the region of negative diffusivity u ∈ (a, b), and the region where the solution is multi-valued are avoided altogether.Inserting a shock is an established technique for obtaining single-valued solutions to hyperbolic equations of the form ρ t +c(ρ)ρ x = 0 [39].Witelski [40] expanded this technique to construct a solution for a diffusion equation with a region of negative diffusivity.Furthermore, some reactiondiffusion models with a region of negative diffusivity admit travelling wave solutions where the density profile moves at constant speed while maintaining its shape [14,23,25,26,29].
Of particular interest to this work, Li et al. [26] proved the existence of shock-fronted travelling waves for one such reaction-diffusion model.
We use a nonclassical Lie symmetry to construct a multi-valued analytic solution to the RDE (1) with a region of negative diffusivity.Lie symmetry analysis was first proposed by Sophus Lie in the late 19th century [27] and involves finding invariant quantities in an equation and using them to construct solutions (see [19,31] for more detail).The nonclassical (or Q-conditional) method is an extension to Lie symmetry analysis introduced by Bluman and Cole [4], where in addition to requiring invariance of the PDE the invariant surface condition must also be satisfied.Nonclassical symmetry analysis can sometimes uncover symmetries not found by the classical method.While boundary and initial conditions are important for determining the behaviour of a system, they are often not included in the initial analysis because calculation of the symmetries becomes overly restrictive [15].We follow this strategy here and apply the symmetry method to the governing equation only, revisiting the initial and boundary conditions after a solution has been found.
Here we use a special nonclassical symmetry found independently by Arrigo and Hill [3] and Goard and Broadbridge [16] to construct analytic solutions for an arbitrary diffusivity provided that the reaction and diffusion terms satisfy the relationship where u * is a zero of the reaction term, and A and κ are free parameters that arise from the symmetry analysis.Generally this reaction term contains singularities (simple poles) when D(u) = 0 (that is, at u = a and u = b), however inserting a shock avoids these singularities.For the applications considered here, the case that u * is a zero of D(u) as well as a zero of R(u), is not considered.This symmetry transforms the RDE (1) to the one-dimensional spatial Helmholtz equation [7], The new dependent variable Ψ(x) is related to the density u(x, t) via the Kirchhoff transformation where Φ(u) is the Kirchhoff variable, sometimes known as the diffusive flux or potential.
After solving the Helmholtz equation ( 3), inverting the Kirchhoff transformation (4) recovers the solution for u(x, t).We restrict our attention to the choice κ = −k 2 < 0, so that the solution to the Helmholtz equation ( 3) is in terms of exponential functions, In summary, we will specify a diffusivity with a region of negative diffusion and calculate the related reaction term (2).We then construct analytic solutions by inverting the transformation (4) using the exact form of Ψ(x) (5).
In Section 2, we investigate a model where D(u) is a quadratic function with two zeros on u ∈ (0, 1), such that there is a region u ∈ (a, b), where 0 < a < b < 1, with D(u) < 0.
We consider three types of analytic solutions.The first two are time-dependent solutions, which can represent both a receding front with non-constant speed and a pair of colliding waves.The third is a receding sharp-fronted travelling wave, where the solution moves leftward at constant speed.All three solution types are sharp-fronted and satisfy a Stefanlike moving boundary condition due to the requirement that the density remains nonnegative [10,11].Since the solutions with negative diffusivity are multi-valued, we create shock-fronted solutions by inserting a jump in u around values for which D(u) < 0, so that the shock conserves Φ ′ (u) = D(u) across the shock or Φ ′ (u) is continuous rule.
After developing the shock-fronted solutions, we analyse the stability of the constant and travelling wave solutions.Section 3 examines the relationship between the Φ ′ (u) is continuous rule and the more common equal area in Φ(u) rule [32,40].We provide discussion and concluding remarks in Section 4. The existence of an analytic solution obtained using nonclassical symmetry analysis underpins the analysis in each section of the paper.

Analytic travelling wave and time-dependent solutions for quadratic diffusivity
We consider three analytic solutions for the RDE (1) with a quadratic diffusivity that is negative only for values of u between the two roots of the quadratic.This diffusivity takes the form where a, b ∈ (0, 1) and we impose a < b.Johnston et al. [20] derived a reaction-diffusion model with diffusivity (6) in the continuum limit of a discrete model where isolated and group-associated individuals had different rates of birth, death and movement.Turchin [37] also derived a diffusion equation with a negative diffusivity from a discrete model where individuals move towards regions of higher population density [37].Li et al. [25,26] then showed that travelling wave solutions exist for the RDE with quadratic diffusivity (6).
We extend these works by using a nonclassical symmetry to construct solutions with shocks, and analyse the spectral stability of travelling wave solutions.Figure 1a shows an example quadratic diffusivity of the form (6), with a = 0.2 and b = 0.4.To use the special nonclassical symmetry, we obtain the corresponding reaction term using (2).By choosing u * = 1, we fix one of the zeros of the reaction term to be at u = 1, that is, Φ(1) = 0, and so R(1) = 0. Since A is a free parameter, we choose A = −κD(0), to enforce R(0) = 0.The reaction term R(u) (2) has a third zero at u = a + b.Therefore, there are two types of reaction terms possible for this quadratic diffusivity.If a + b ≥ 1, the reaction term is typically negative for all values of u in the range of interest.Since a negative reaction term would cause a decrease in population for all densities, we do not investigate this case further.If a + b < 1, the reaction term is typically negative for lower densities and positive for higher densities.This scenario is loosely analogous to cubic reaction term commonly used to model changes in gene frequencies [6], or population dynamics due to the Allee effect [8,36]. Figure 1b shows the shape of this reaction term.The reaction term has singularities at a and b, as the solid vertical lines in Figure 1b indicate.Between the vertical lines the reaction term becomes large and positive.However, the shock solutions we construct avoid the u values in this region.We construct an analytic solution to (1) with diffusivity (6) and reaction term (2) by solving the Helmholtz equation (3) for Ψ(x), and inverting the Kirchhoff transformation (4) to find u(x, t).Since the diffusivity is quadratic, the Kirchhoff variable Φ(u) will be cubic in u(x, t).The Kirchhoff transformation (4) then determines the solution implicitly.
That is, where c 1 and c 2 are constants of integration whose values will be governed by the boundary conditions, and where we have chosen κ = −k 2 < 0. Some sections of the solution for u(x, t) are multi-valued, because Φ(u) is a cubic function.Since A = −κD(0) > 0, we require Ψ(x) < 0 so that the solution does not approach infinity as t → ∞.As a consequence, we require Φ(u) < 0 for u ∈ [0, 1), so that we must have b < (a + 2)/3.
Figure 1c shows Φ(u) which is increasing except where u ∈ (a, b) where the diffusivity is negative.Consequently the solution is multi-valued for values of u surrounding this region.

Constructing single-valued solutions by inserting shock discontinuities
With quadratic diffusivity (6), the solutions (7) obtained using the nonclassical symmetry solutions described here are multi-valued with three branches in a region centred about the points where the diffusivity is zero, i.e.where u = a and u = b.That is, the To recover a single-valued solution, we insert a shock discontinuity.Due to the form of transformation (4), the flux potential Φ(u) always conserved across the shock so that in principle, a shock could be inserted anywhere in the multi-valued region.As a result of the symmetry Φ(u) x is always conserved across the shock.To specify a second condition, we notice that equation ( 1) can be rewritten in the form and consequently we impose that Φ ′ (u) = D(u) is also continuous across the shock.That is, the shock position is determined by the following two conditions where u l and u r are the left (lower) and right (upper) endpoints of the shock.For quadratic diffusivity ( 6), these points are where we also assume that b < a(2 + √ 3) so that u l > 0. Due to the relationship between the diffusivity and the reaction term (2), R(u) is also preserved across the shock.In addition, both u x and u t are preserved across the shock.Figure 2 shows the behaviour of both the diffusivity and the reaction across the shock.In Figure 2a, the dashed horizontal line shows that the solution jumps over both the zeros of the diffusivity, u = a, b and the region where in Figure 2b the solution jumps over the locations where the reaction term is singular, and the region between the vertical lines where the reaction term is positive.Figure 2c shows the flux potential around the shock, demonstrating that Φ(u) is conserved across the shock.

Examples of analytic time-dependent solutions
We obtain different analytic solutions by adjusting the constants c 1 and c 2 in (7).For example, setting c 1 < 0 and c 2 = −c 1 gives rise to a receding solution that satisfies u(0, t) = 1 for all t.Although the receding solution u(x, t) gets steeper as the front approaches x = 0, the shock persists because the underlying analytic solution remains multi-valued.This solution is shown in Figure 3a.Alternatively, if we remove the relationship between c 1 and c 2 and instead merely impose c 1 < 0 and c 2 < 0, we obtain the  Given that the density u(x, t) typically represents a biological or chemical concentration, we usually require u(x, t) ≥ 0 in feasible solutions to reaction-diffusion equations.The nonclassical symmetry solutions presented in Figure 3 do not have u ≥ 0 for all values of x.However, we can obtain solutions with non-negative u by restricting the domain [10].
For example, the time-dependent receding solution in Figure 3a where β is the Stefan parameter.We can calculate both the speed of the moving boundary and the flux at the moving boundary exactly using our analytic solutions.Doing so reveals that our solutions do not satisfy the regular Stefan condition (11), but instead satisfy an alternate condition whereby the speed is inversely proportional to the gradient, Both the time-dependent and colliding wave solutions in Figure 3 satisfy this condition (in the case of colliding waves, both boundaries satisfy the condition).Although the physical or biological interpretation of this Stefan-like condition (12) remains to be determined, moving boundary conditions with boundary speed proportional to the reciprocal of the density gradient have been reported in the literature [28].
One consequence of our Stefan-like condition (12) is that L ′ (t) < 0 in travelling wave and time-dependent solution.This is because the gradient u ′ (x) is negative at x = L(t) (see Figures 3a and 4a), and κ and Φ(0) are both negative constants so that κΦ(0) is positive.This gives rise to receding density profiles.Although receding fronts are uncommon in reaction-diffusion models, the receding behaviour is similar to that seen in recent works that adapt the classical Stefan problem to solve reaction-diffusion (Fisher-KPP) equations on moving boundaries [9,11,12,13].Due to the inverse relationship between the gradient and the speed of the boundary (12), steeper density gradients lead to slow-moving fronts, whereas shallower gradients increase front speed.Large flux at the boundary (D(u)u x )| x=L(t) thus corresponds to slow fronts, and small flux corresponds to fast fronts.

Analysis of a receding travelling wave solution
While the time-dependent solutions in Section 2 involve non-zero c 1 and c 2 , we can obtain a travelling wave solution by choosing c 1 < 0, c 2 = 0.The implicit solution (7) for Φ(u) then becomes We can write Equation ( 13) in terms of the travelling wave coordinate z = x − ct, where c = −A/k = −kD(0).The travelling wave solution then has an implicit form, expressed in terms of z as As described above, the travelling wave solution (14) obtained using the nonclassical symmetry is multi-valued around the region where the diffusivity takes negative values.
Figure 4 shows a shock inserted so that the diffusion is constant across the shock, that is, according to conditions (9).The start and end of the shock are marked by stars (red) while the horizontal lines (red) show the zeros of the diffusivity; between these lines, the diffusivity is negative.The arrow shows that this travelling wave solution moves leftward (c < 0), and is a receding front.This is a consequence of requiring R(0) = 0, which means that A = −κD(0) > 0. The shock avoids both the region of negative diffusivity and the singularities in the reaction term, and both the diffusivity and reaction terms are continuous across the shock.The behaviour of D(u) and R(u) across the shock corresponds to the horizontal dashed lines in Figure 2a and Figure 2b respectively.
Like the time-dependent solutions in Figure 3, the travelling wave solution in Figure 4 also satisfies the Stefan-like condition (12) at the location where u = 0. Choosing c 1 = Φ(0) means that the sharp front in u coincides with z = 0 (that is, u(0) = 0).Figure 5a shows For the colliding wave solution, the dash-dot curves (green) in Figure 5 are multi-valued because of the two contact lines; as the two boundaries approach zero the speed of the moving boundaries starts to increase and the flux at both boundaries approaches zero.
The time-dependent and travelling wave solutions are both left-moving (receding) solutions.Li et al. [26] derive a necessary condition for the existence of left-moving travelling waves satisfying the usual far-field boundary conditions (that is, u → 0, 1 and u z → 0 as z → ±∞ respectively).For a sharp-fronted wave, the condition of Li et al. [26] changes to  where g(u) = D(u)u z .Consequently, while a model with a reaction term that is positive in (0, u a ) cannot support a smooth-fronted receding travelling wave (that is, a solution where u z → 0 as z → ∞), a compactly-supported sharp-fronted receding travelling wave (for which g(0) ̸ = 0) can be supported.
To analyse the receding travelling wave solution further, we apply the transformation z = x − ct to the RDE (1) to obtain the ordinary differential equation The implicit solution ( 14) with c 2 = 0 is a travelling wave solution u(z) that satisfies (17).
Following Li et al. [25], we write (17) as the system The left-hand sides of (18) vanish when D(u) = 0, giving rise to a wall of singularities at u = a and u = b [18,25,33,38], i.e., solutions of (18) will, in general, cease to exist when reaching these walls of singularities due to a finite-time blow-up.Possible exceptions are points known as holes in the walls [18,25,33,38] where the left-hand and right-hand sides of (18) both vanish.However, the conditions D(u) = 0 and q = 0 do not correspond to holes in the wall; although the right-hand side of (18a) vanishes when q = 0, the right-hand side of (18b) does not vanish when D(u) = 0 since the product is non-zero at the walls of singularities, i.e., A u=a,b u * D(u ′ ) du ′ ̸ = 0. Thus, smooth solutions between u = 0 and u = 1 are not possible in system (18) and shocks are a necessary ingredient.
Nevertheless, we can use our analytic multi-valued solution to calculate solution trajectories in the phase space.Figure 6 shows a phase plane portrait of ( 18) with the trajectory of our multi-valued solution.The phase plane also shows why a shock is required: between the walls the diffusivity is negative and the direction of the flow switches.As a consequence, our solution cannot be represented by a single trajectory and instead is a combination of two solution branches.The phase portrait in Figure 6 also demonstrates that the sharp-fronted receding travelling wave has non-zero flux at the moving boundary (where u = 0).

Stability analysis
We analyse the spectral stability of the constant and travelling wave solutions to (1), relative to perturbations in an appropriately chosen space.To that end, we first derive the linearised equation, linearised about a solution ū, which we write in the form where p is a perturbation to the solution ū and L(ū) is the linearised operator.To derive the linearised equation, we consider solutions to (1) of the form u = ū + εp(x, t), where For a general solution of (1), the right hand side of (19), and in particular the operator L(ū) will depend on both x and t.However, for the special solutions we consider, L will either be a constant coefficient operator, or depend on a single variable only.

Linear stability analysis of the constant solutions
We suppose first that our solution ū is a constant solution, either 0 or 1.Then, (19) becomes a constant-coefficient second-order PDE, which is solvable via separation of variables.Supposing that p(x, t) takes the form p(x, t) = e λt p(x) we obtain the ordinary differential equation We choose the domain of our perturbations to be p ∈ L 2 (R), which provides the boundary conditions lim x→±∞ p = 0, and enables us to solve (20) via the Fourier transform.We are led to the dispersion relation, where α is the Fourier variable in space.If Re(λ) < 0 then, to first order, the perturbation will decay to the constant solution ū, indicating that u = 0, or u = 1 is linearly stable to L 2 (R) perturbations.In contrast, if Re(λ) > 0 the perturbation to the constant solution ū will grow with time.

The derivative of the reaction term (with
At u = 0 we have Since κ < 0, Φ(0) < 0, and D(0) > 0, the derivative R ′ (0) has the same sign as D ′ (0) = −a − b, which is always negative as 0 < a < b.Therefore, small perturbations around the constant solution u = 0 will decay, and the constant solution u = 0 is (linearly) stable.
Likewise, we can consider the stability of the constant solution u = 1.Since D(1) > 0, stability depends on the value of R ′ (1), which is Since −κD(1) > 0 for quadratic diffusivity, the sign of R ′ (1) depends on the sign of (D(0)/D(1) − 1) .This term is positive if D(1) < D(0), negative if D(0) < D (1), and is is zero when a + b = 1.Quadratic diffusivity with a + b > 1 thus admits long-wavelength instabilities for sufficiently small α.The solutions presented in Figure 4 are for a + b < 1, for which so u = 1 is stable relative to perturbations in L 2 (R).

Stability analysis of the receding travelling wave solution
Having considered constant ū, we now address the linearisation of (1) in the moving frame (z, t) = (x − ct, t), about the receding travelling wave solution described in Section 2.3.
Again, we subject the solution to perturbations of the form p(z)e λt , where p(z) will be in an appropriate space.Our equation for the perturbation is The right hand side of ( 25) depends only on z.We require that p(z) be a solution to (25) which is square integrable on (−∞, 0) which also vanishes at z = 0 (i.e. in L 2 0 ((−∞, 0))).
To uniquely solve the linearised ODE (25), we need to impose conditions across the shock.
To this end, we ask that p be C 1 across the shock.We note that this means that the coefficient of the highest derivative of D(ū) will be continuous across the shock, but the lower order coefficients may not be.We remark that since the shock occurs at a single point, the set of functions which are C 1 across the shock is still a densely defined subset of L 2 0 ((−∞, 0)), so imposing this condition will not affect the spectrum of L(ū).
For our discussion on spectral stability of the receding travelling wave solution, we use the language and notation outlined in Chapters 2 and 3 of the textbook by Kapitula and Promislow [21].The set of λ ∈ C such that L(ū) − λ does not have a bounded inverse on L 2 0 ((−∞, 0)), will be called the spectrum of L(ū), and denoted σ(L(ū)).We note that the spectrum of a closed linear operator on a Hilbert space is a (topologically) closed set, so if we can find an open set which is in the spectrum, then its closure must be as well.We will call the wave ū spectrally stable, if λ ∈ σ(L(ū)) means that Re(λ) < 0.
Remark.Sometimes waves are referred to as spectrally stable if they have spectrum in the left-half plane, and perhaps a simple eigenvalue associated with translation invariance at the origin.As our perturbations are in L 2 0 ((−∞, 0)), the eigenvalue associated with translation invariance λ = 0 will not be an eigenvalue of the linearised operator.
Proof.We begin with an analysis of the so-called far-field operator of (25).Letting z → −∞, we have that solutions to (25) will asymptotically behave like solutions to the constant coefficient ordinary differential equation where we recall that c = −kD(0).The characteristic equation of ODE ( 26) is with characteristic exponents given as the roots

1) .
We are interested in the sign of the real parts of r ± for a given value of λ.Supposing that r = iα for α ∈ R, we have the dispersion relation Equation ( 28) is a parametrised curve for λ in the complex plane, parametrised by α, and divides the complex plane into two disjoint regions.We denote the region to the right of the curve by Ω, while the region to the left of the curve we will call B, see Figure 7.For λ ∈ Ω we have that Re(r + ) > 0, while Re(r − ) < 0. However, for λ ∈ B we have that both Re(r ± ) > 0, meaning that L(ū) − λ is not invertible here.As both r ± have positive real parts, we have a spanning set of solutions to (25) which decay to 0 as z → −∞.Then the solution satisfying the right-hand boundary condition would necessarily be a linear combination of these solutions, and would hence lie in the kernel of L(ū) − λ.Owing to our remark above about the spectrum being closed, we conclude that the spectrum at least contains the closure of the set B.
We have that D(1) > 0, so the maximum real value of λ in (28) depends on R ′ (1), which we know is negative provided that the conditions in the statement of Theorem 2.1 are satisfied.Thus this part of the spectrum at least is contained entirely in the left half plane.Remark.In fact, what the above calculation partially shows is that the so-called Fred-holm index of L(ū) − λ is not zero for λ ∈ B (it is in fact 2), but as we can show lack of invertibility of L(ū) − λ for λ ∈ B directly, we avoid the diversion into functional analysis terminology.For a more complete description of why the previous calculation computes the Fredholm index of L(ū) − λ, we refer the reader to [21].
The previous computation shows that for λ ∈ Ω we can compute the Green's function of (25), and thus invert L(ū) − λ, provided we do not have an eigenvalue/eigenfunction to (25).We first show that any eigenvalues must be real, and then show that for any a and b satisfying the conditions in the statement of Theorem 2.1 the eigenvalues must be negative.
To see that eigenvalues of ( 25) with λ ∈ Ω must be real, we multiply (25) by the integrating factor so that equation ( 25) is put into Sturm-Liouville form where the three z-dependent functions are Noting that S 3 is a positive, bounded, continuous function (even across the shock because u always takes values where D(u) > 0 and D(u) is continuous across the shock (9)), a calculation shows that for λ ∈ Ω, if λ is an eigenvalue of L(ū), then it will be an eigenvalue of 1 S 3 (z) M in the weighted Hilbert space (L 2 0 ((−∞, 0))) S 3 with the (usual) weighted norm < v, v > S 3 := vvS 3 dz.
Lastly to see that any eigenvalues λ ∈ Ω must be negative, we note that if λ is real, the corresponding eigenfunction will also be real.We multiply ( 25) by p and integrate, where z s is the location of the shock.Integrating the speed terms results in where p(z − s ) is the value of p(z) as z → z s from below and p(z + s ) is the value of p(z) as z → z s from above.For the diffusion term, using integration by parts twice we arrive at where We have D(ū) z = D ′ (ū)ū z , and for quadratic diffusion D ′ (ū) > 0 when ū > (a+b)/2 which is in the region to the left of the shock z ∈ (−∞, z s ) and D ′ (ū) < 0 when ū < (a + b)/2 which is in the region to the right of the shock z ∈ (z s , 0].Also ūz < 0, therefore will vanish if p is continuous at the shock, and because D(ū) is conserved across the shock the term [D(ū)pp z ] z=z − s − [D(ū)pp z ] z=z + s will vanish if p z is also continuous across the shock.Alternatively this entire term will vanish if p also vanishes at the shock.
With S < 0 the right-hand side of Equation (35) will always be negative since D(ū) is always positive (since the shock avoids the regions of negative diffusivity) and D(ū) zz /2 + R ′ (ū) is assumed to be negative.We then need to have λ < 0 for the left-hand side of (35) to be negative as well.Hence, any real eigenvalues of L(ū) in Ω will be negative.This completes the proof of Theorem 2.1.
We have proven that receding travelling waves with quadratic diffusivity are spectrally stable whenever D(ū) zz /2 + R ′ (ū) < 0. We now explore numerically the values of a and b for which the stability analysis applies. Figure 8

The equal area rule and non-symmetric diffusivity
As described in Section 2.1, our shock solutions are constructed by connecting two branches of an analytic multi-valued solution; here we choose the shock location according to the conditions described in equation ( 9) and as shown in Figure 9.In practice, the shock could be located anywhere in the multi-valued region with two conditions required to locate the shock.The first condition in (9) requires the diffusive flux to be continuous across the shock and is a constraint used almost universally (see for example [26,40]).In previous works [32,40,41], the second condition has been chosen to correspond to the equal area in Φ(u) rule, ur u l Φ(u) − Φ(u l ) du = 0 (37) as illustrated in Figure 9b.The horizontal line indicates the value of Φ(u) at the shock chosen so that the two shaded regions have equal area.Pego [32] used the equal-area rule to describe solutions of the nonlinear Cahn-Hilliard equation for phase separation, and Witelski [40] showed how the equal-area rule arises when constructing a shock solution to the diffusion equation using a non-local regularisation where ε ≪ 1 is a small parameter.Li et al. [26] showed that the equal-area rule also arises as a shock solution for reaction-diffusion equations using the same type of regularisation (that is, equation (38) with a reaction term added).Other types of regularisations, including composite regularisations, have also been used to construct shock solutions [26,5].
In Section 2 the location of the shock was determined by the requirement that the derivative of the flux potential Φ ′ (u) must also be continuous, that is, conditions (9) must be satisfied across the shock.As a consequence, the diffusivity is constant across the shock (in the case of the symmetry solution presented here, relationship (2) means that the reaction term is also constant across the shock, but this may not be true in general).In the case of quadratic diffusivity above, the condition on the first derivative recovers the more commonly imposed equal area rule (37).Here we show that this is true for any diffusivity that has two zeros and is symmetric about the midpoint of its zeros.
is odd.Also, since Φ ′ (ū) = D(ū) is even, the requirement that Φ ′ ( ūl ) = Φ ′ ( ūr ) is satisfied by ūl = − ūr , because Φ is a vertical translation of the odd function f (ū) and the requirement that Φ( ūl ) = Φ(− ūl ) means that Φ( ūl ) = Φ| ū=0 .The left-hand side of the equal area condition (37) then becomes where f (ū) dū = F (ū), and F (ū) is an even function since f (ū) is odd.That is, the equal area rule is equivalent to the requirement that Φ ′ (u) be continuous, providing the diffusivity is symmetric about the midpoint of its two zeros.In contrast, if the diffusivity is not symmetric about the midpoint of its two zeros, the requirement that Φ ′ (u) is continuous does not recover the familiar equal area rule.For example, consider a quartic diffusivity given by   to resolve regions where u(x, t) is multi-valued.We place shocks such that the flux potential Φ(u) and its derivative Φ ′ (u) are continuous across the shock.For diffusivity that is symmetric about the midpoint of its zeros (e.g. the quadratic (6)), these conditions are equivalent to the more commonly-used equal-area rule.With symmetric diffusivity, the continuity and equal-area requirements recover the shock with longest possible length.For a non-symmetric diffusivity, the largest jump coincides with the continuity condition only.Since another way to write equation ( 1) is equation ( 8), the requirement that Φ ′ (u) is continuous has stronger foundation than imposing equal-area.Our planned future investigation of the entropy jump across the shock and regularisation to the PDE might yield further understanding about the link between negative diffusivity and physical phenomena.
The solutions presented here avoid the values of u(x, t) for which the diffusivity becomes negative.It is interesting then that the solutions do not display the usual characteristics of a standard diffusion model.This may be a consequence of the loss of mass at the moving boundary, although it is more likely that the negative values of D(u) impact the solution despite the solution not taking on the relevant values of u.While the travelling wave solutions presented are exact solutions that satisfy the PDE, they are not weak solutions in the usual sense of averaging through the PDE with an arbitrary compactly-supported test function.We refer the reader to Whitham [39] for further information about weak solutions that contain shocks.
In addition to presenting the analytic solutions, we perform spectral stability analysis for constant and travelling wave solutions.We show that the constant solution u = 0 is stable, but that u = 1 admits long-wavelength instability for quadratic diffusivity with a + b > 1.We also prove spectral stability of travelling wave solutions to (1) with quadratic diffusivity (6), for some combinations of a and b.However, open problems remain in the analysis of analytic solutions to reaction-diffusion equations with negative diffusivity.For non-symmetric diffusivity, the continuity of flux Φ(u) and its derivative Φ ′ (u) conditions for defining shocks do not yield an equal-area rule.The regularisation to the nonlinear reaction-diffusion equation (1) corresponding to the requirement for Φ ′ (u) to be continuous remains unknown.Investigation of this and other regularisations will be the subject of future work, as will a study of the jump in entropy across the shock.

Figure 2 :Figure 3 :
Figure 2: Nonlinear (a) diffusivity, (b) reaction function, and (c) flux potential around the shock.Here, κ = −1, A = 0.08, u * = 1, a = 0.2 and b = 0.4.The vertical lines (red) show the location of the zeros of the diffusivity.Horizontal dashed lines indicate the transition across the shock, and the stars (red) indicate the end points of the shock.The reaction term is zero at u = 0, u = a + b = 0.6, u = 1.

Figure 4 :
Figure 4: Travelling wave solution.Curves (blue) show the solution at different times and the arrow indicates how the solution moves with time.Stars (red) indicated the beginning and end of the vertical shock.Horizontal lines (red) show where the diffusivity is zero.(b) The dashed curve shows the multi-valued part of the solution at a single time.Here κ = −1, A = 0.08, u * = 1, a = 0.2, b = 0.4, c 1 = Φ(0) and c 2 = 0.

Figure 5 :
Figure 5: Comparison of the boundary position L(t) and the flux −D(u)u x at x = L(t) for the receding time-dependent, colliding wave, and travelling wave solutions.The dash-dot curves (green) are multi-valued because it corresponds to the colliding wave solution and has two boundaries.Solutions obtained using parameter values κ = −1, A = 0.08, u * = 1, a = 0.2, b = 0.4.For the time-dependent solution, we set c 1 = −c 2 = −0.0153,for the colliding wave solution we have c 1 = c 2 = −0.0153,and for the travelling wave solution we use c 1 = −0.0153,c 2 = 0.

Figure 6 :
Figure 6: Phase portrait of(18).The arrows (green) indicate the direction field of the system.The dashed curve (magenta) is the solution trajectory of our analytic solution, which emanates from (1, 0) and terminates at (0, −q 0 ) satisfying the boundary condition(12).The vertical lines (red) are walls of singularities.The solid horizontal line and the solid curve (black) are the nullclines of the system.

Figure 7 :
Figure 7: Schematic of the dispersion relation (28) in the complex plane dividing the plane into two disjoint regions, denoted B and Ω.

Figure 8 :
Figure 8: Values of a and b for which a + b < 1 and 0 < a < b < 1. Circles (blue) indicate the pairs that satisfy the conditions for spectral stability outlined in Theorem 2.1, and the condition b < a(2+ √ 3), which guarantees u l > 0 (10).Crosses (red) indicate combinations of a and b that do not satisfy both conditions.Our analysis proves spectral stability for (a, b) pairs denoted by blue circles, but makes no statement about the stability of (a, b) pairs denoted by red crosses.

Figure 9 :
Figure 9: (a) Stars (red) show the shock position for the solutions shown in Figures 3 and 4 and the dashed line (blue) shows that Φ ′ (u) = D(u) is conserved across the shock.(b) The equal area rule for Φ(u).The curve (blue) is Φ(u) and the horizontal line (red) is the value of Φ(u) at the shocks, Φ(u l ) = Φ(u r ).The equal area rule in Φ(u) requires that the areas in the two lobes cut off by the horizontal line (the shaded regions) are equal.
)where a < b, a, b ∈ (0, 1), c arbitrary and d > 0 are constants.This diffusivity has two zeros (D(a) = D(b) = 0), and is not symmetric about their midpoint.Figure10(a) shows this diffusivity, where the parameters have been chosen to mimic the diffusivity examined in Section 2. Figure10(b) shows a possible corresponding (through (2)) reaction-term chosen to have similar characteristics to that discussed in Section 2; vertical red lines indicate the location of the zeros of the diffusivity.In both panels, the horizontal blue dashed line indicates the transition across the shock as chosen by enforcing that Φ ′ (u) is continuous.A direct consequence of this choice is that both D(u) and R(u) are continuous across the shock.The green dash-dot lines indicate the transition across the shock if the equal area rule is imposed.Notice that location of the shock is different, and also that neither the diffusivity nor reaction is continuous across a shock chosen in this way.An implicit solution to equation (1) with (40) and (2) can be constructed using the nonclassical symmetry described in Section 1. Figures10c and 10dshow a receding travelling wave solution.The dashed line in panel (d) indicates the difference in shock location as chosen by the requirement that Φ ′ (u) is continuous and the equal area rule.

Figure 10 :
Figure 10: Analytic solution to reaction-diffusion equation (1) with non-symmetric diffusivity (40).Parameter values are κ = −1, A = 0.0448, a = 0.2, b = 0.4, c = 0.6, and d = 0.2.(a) Diffusivity D(u) (40).Inset shows the transition across the shock: horizontal dashed line (light blue) indicates the position chosen if Φ ′ (u) is continuous, dash-dot line (green) shows the position chosen by the equal area rule.(b) Corresponding reaction term R(u), given by (2).(c) Multi-valued travelling wave solution.Curves (blue) show the solution at different times, the arrow indicates the direction of travel.Horizontal lines (red) show where the diffusivity is zero.(d) Enlargement of the multi-valued part of the solution at a single time.Stars (red) show the end points of two possible shock positions.The dashed line (left, light blue) corresponds to Φ ′ (u) is continuous.The dash-dot line (right, green) corresponds to the equal-area rule.