On the relativistic heat equation in one space dimension

We study the relativistic heat equation in one space dimension. We prove a local regularity result when the initial datum is locally Lipschitz in its support. We propose a numerical scheme that captures the known features of the solutions and allows for analysing further properties of their qualitative behavior.


Introduction
In this work, we explore both analytically and numerically the implications of a new strategy to study flux-dominated nonlinear diffusions in one dimension.To be more precise, we consider the so-called relativistic heat equation (RHE) x ∈ IR, t > 0.
(1.1) introduced by Rosenau in [37] and, later on, by Brenier in [14] based on optimal transportation ideas.The name of RHE comes from the fact that (1.1) converges as c → ∞ to the heat equation both formally and rigorously [20], while the flux in (1.1), understood as a conservation law, is bounded by the speed of light c whenever the solution is positive.Many other models of nonlinear degenerate parabolic equations with flux saturation as the gradient becomes unbounded have been proposed by Rosenau and his coworkers [19,37], and Bertsch and Dal Passo [12,26].Notice also [36] for the presence of flux limited diffusion equations in the context of radiation hydrodynamics.
The general class of flux limited diffusion equations and the properties of the relativistic heat equation have been studied in a series of papers [5,4,6,21].An existence and uniqueness theory of entropy solutions for the Cauchy problem associated to the quasi-linear parabolic equation was developed in [5,4].Here, the flux function is given by b(z, ξ) = ∇ ξ f (z, ξ) and f : IR × IR N → IR + is a convex function with linear growth as ξ → ∞, such that ∇ ξ f (z, ξ) ∈ C(IR × IR N ) satisfying other additional technical assumptions.In particular, the relativistic heat equation (1.1) satisfies these assumptions, and other models considered in [37].To avoid the difficulty of the lack of a-priori estimates that ensure the compactness in time of solutions of (1.2), the existence problem was approached using Crandall-Liggett's theorem [24].For that, we first considered the associated elliptic problem and we defined a notion of entropy solution for which we developed a well-posedness theory.The notion of entropy solution permits to prove a uniqueness result using Kruzhkov's doubling variables technique [30,15].This technique was suitably adapted to work with functions whose truncatures are of bounded variation [5,4], which is the natural functional setting for (1.2) and its associated elliptic equation.
The evolution of the support of solutions of the relativistic heat equation (1.1) was studied in [6].By constructing sub-and super-solutions which are fronts evolving at speed c and using a comparison principle between entropy solutions and sub-and super-solutions, it was proved in [6] that the support of solutions evolves at speed c.Moreover, the existence of solutions which have discontinuity fronts moving at the speed c was again shown using the comparison principle with sub-solutions.This implies, in particular, that the maximal regularity in time that one can expect for general solutions of (1.1) is that u ∈ BV ([τ, T ] × IR N ) for any 0 < τ < T .That this happens for a general class of initial conditions was proved in [8] and later extended in [21].This lack of regularity is at the origin of the notion of entropy solutions for this type of equations.The only regularity result for smooth initial conditions was proved in [20] and it guarantees that ∇ ln u is bounded whenever initially is.But the study of the local regularity of solutions of (1.1) is still an open question.One of the purposes of this paper is to address this problem for the Cauchy problem associated to (1.1) in one space dimension with compactly supported bounded probability densities as initial data.
Assuming that the initial data in non-negative, we can easily change variables to observe that ũ(t, x) is a solution of (1.1) if and only if u(t, x) = ũ( ν c 2 t, ν c x) is a solution of x .
(1.3) Thus, without loss of generality we may assume that ν = c = 1, and for simplicity we shall assume it in the rest unless explicitly stated.Notice also that if u(t) is a solution corresponding to u 0 , then λu(t) is a solution corresponding to λu 0 , λ > 0. Thus, without loss of generality we assume that u 0 1 = u(t) 1 = 1 for any t > 0, and reduce our evolution to probability densities.In this paragraph, the term solution refers to entropy solution for which the well-posedness theory was developed and for which a summary of its concept is reminded to the reader in the Appendix.
The local regularity of entropy solutions to (1.3) will be done by a change variables, writing (1.3) in terms of its inverse distribution function.This change of variables has its origin in using mass transport techniques to study diffusion equations [18,13].It is known [14] that equation (1.1) has the structure of a gradient flow of a certain functional (the physical entropy) with respect to some transport distance.This structure was already used to give well-posedness results to (1.1) in [35].Nonlinear diffusions have received lots of attention from optimal transport theory viewpoint starting from the seminal works [29,33].
Transport distances between probability measures in one dimension are much easier to compute since they can be written in terms of distribution functions and their generalized inverses (pseudo-inverse), the so-called Hoeffding-Fréchet Lemma [39,Section 2.2].This result led to the following change of variables based on the distribution function F associated to the probability measure u, defined as We formally consider its inverse ϕ defined on the mass variable η ∈ (0, 1) that verifies After straightforward computations assuming that all involved functions are well-defined and smooth, one obtains the equation for the inverse distribution function ϕ.This change of variables has first been used for nonlinear diffusions in [18] to show contractivity properties of transport distances for porous-medium like equations.
It is worthy to remark that an implicit Euler discretization of (1.4) is equivalent to the variational JKO scheme whose convergence is proved in [35] for (1.1) under certain assumptions.Numerical schemes to solve the equation for the pseudo-inverse function in the case of the porous medium equation were analysed in [28].This Lagrangian approach was generalized to several dimensions in [16] in order to propose numerical schemes for equations with gradient flow structure in optimal transport theory and general quasilinear problems in divergence form.
In Section 2, we will first take advantage of this change of variables to prove the following regularity result: t ∈ (0, T ), for almost any t ∈ (0, T ), and u(t) is smooth inside its support, We emphasize that the new parts of this result with respect to the literature discussed above refer to the regularity stated on points (ii) and (iii).This result implies that sharp corners on the support of the initial data are immediately smoothed out by the evolution of the RHE.This result will be extended in Section 3, in particular, we cover the case where the initial condition u 0 vanishes at the boundary of its support.
In Section 4, we will propose an adaption of the numerical scheme in [16] based on equation (1.4) with suitable boundary conditions that fully captures the demonstrated behavior of the solutions of the RHE.Moreover, we will show different numerical tests in situations where the theory has not been developed yet.For instance, we numerically study the conditions for the formation or not of discontinuities on the bulk of the solutions for RHE and its porous medium counterparts x with m > 1 and their long-time asymptotic behaviour.Finally, we include in Appendix A some basic material to describe the notion of entropy solutions for (1.3) for the sake of completeness.

Regularity of Solutions
As proved in [5], there exists a unique entropy solution of the Cauchy problem for (1.3) for any u 0 ∈ L 1 (IR)∩L ∞ (IR), u 0 ≥ 0, see Appendix for the full notion of solution.Moreover if u 0 has compact support in IR and is locally bounded away from zero in any interior point of its support, then supp(u(t)) = supp(u 0 ) ⊕ B(0, t) [6].The rest of this Section is devoted to the proof of the regularity statements (ii) and (iii).
Let us recall that the entropy condition on the jump set of u can be expressed by saying that the profile of u is vertical at those points.Since the support of u(t) is (a − t, b + t), and u(t) ≥ κ(t) > 0 in (a − t, b + t) for any t > 0 [6], there is a jump at the points x = a − t, b + t and we have [21] Let us consider the change of variables discussed in the introduction and define the function ϕ(t, η) by the relation Proceeding formally, assuming that the function is smooth inside its support and differentiating with respect to η we obtain u(t, x)ϕ η = 1, for x = ϕ(t, η).
Differentiating with respect to t we have Taking into account the boundary conditions (2.1) [21], one has Then the equation satisfied by ϕ is

Regularity result in mass variables
Now, let us consider the change of variables v = ϕ η .The equation satisfied by v is where we have written x instead of η.This will done through this subsection for convenience.
Proof.To prove this claim, we consider the following approximated Cauchy problem where ǫ > 0. The proof is divided in several Steps.In Steps 1 to 3 we prove some formal estimates that are also useful to state the existence of solutions of (2.5)-(2.6) in Step 4. For simplicity we write Let us observe that a(z, ξ)ξ ≥ |ξ| − z 2 . (2.7) Step 1. L p bounds on v for p ∈ [1, ∞).Let us first consider the evolution of the L 1 norm.For that we integrate (2.5) on (0, 1).We have and thus, where the inequality holds in one dimension.Using (2.7) we have Using this recurrence relation, by Gronwall's inequality, we obtain that and that where the constant C(T, p) does not depend on ǫ.
Step 2. L ∞ bounds above and below on v independent of ǫ.Let us construct a supersolution to the Cauchy problem (2.5)-(2.6).Let V (t, x) = B(t) − ǫ 2/3 + x(1 − x) with B smooth and increasing.Take B(0) such that We compute where where C does not depend on ǫ ∈ (0, 1].Take B ′ (t) ≥ C, for instance B(t) = B(0) + Ct.Let us prove that given T > 0, for ǫ > 0 small enough V (t, x) satisfies for ǫ > 0 small enough, and analogously at x = 1.Since V (t, x) is a supersolution for the Cauchy problem (2.5)-(2.6),by the classical comparison principle we get v ≤ V in [0, T ] × [0, 1], and thus there exists M > 0 depending only on u 0 and T such that Let us finally observe that v ≥ α 1 .Indeed, v = α 1 is a subsolution for the Cauchy problem (2.5)-(2.6)and by the comparison principle in its weak version, we deduce that v ≥ α 1 .
Step 3. L p bounds on v x independent of ǫ.Putting together the estimates in Step 2 and (2.9), we deduce that Step 4. Existence of smooth solutions for the Cauchy problem (2.5)-(2.6).The existence of solutions of (2.5)-(2.6)follows from classical results in [31] and [32,Theorem 13.24].We note that thanks to the a priori bounds stated above we could use the flux so that the assumptions of the existence theorems in [31] and [32,Theorem 13.24] hold.Finally, observe that we need to assume a compatibility condition on v 0 so that v 0 satisfies (2.6).If v 0 does not satisfy (2.6), we modify it to define a function v 0,ǫ ∈ W 1,∞ (0, 1) satisfying (2.6).This modification is only done in a neighborhood of x ∈ ∂(0, 1) which vanishes as ǫ → 0+, so that v 0,ǫ is locally Lipschitz inside (0, 1) with bounds independent of ǫ.Finally, we observe that this modification can be done in such a way that sup Although we omit the details of the construction, let us check that (2.10) is compatible with (2.6).
For that, notice that we can take ). Indeed substituting this expression in (2.6), we have

An asymptotic expansion shows
), and thus (2.10) is compatible with (2.6).Let v ǫ be the solution of the Cauchy problem (2.5)-(2.6).Then v ǫ has first derivatives Holder continuous up to the boundary and for g = v ǫxx , v ǫt , we have sup for some α, δ > 0, where P is the parabolic boundary of (0, 1)×(0, T ), that is [0, 1]×{0}∪{0, 1}×(0, T ), and d(•, P) denotes the distance to P. On the other hand, by the interior regularity theorem [31, Chapter V, Theorem 3.1], the solution is infinitely smooth in the interior of the domain.At this point the smoothness bounds depend on ǫ.
Step 5. A local Lipschitz bound on v ǫ uniform on ǫ.For simplicity of notation, let us write v instead of v ǫ .Let w = |v x | 2 φ 2 where φ ≥ 0 is smooth with compact support [a, b] ⊂ (0, 1).This Step is a consequence of the following inequality where A, B are smooth functions, C = (12 + ǫ 2 ), and 0 x , where P is a polynomial in v of degree 3. Assume for the moment that the last term permits to write (2.11) as wt ≤ A(t, x) wxx +B(t, x) wx .Then, using the maximum principle, this implies Let us now prove the claim (2.11).We first compute We also compute w Differentiating (2.5) with respect to x and multiplying by φ 2 we obtain Now, we get where X = a ξξ v xx v 2 x φφ x .Similarly, we obtain Finally, let us compute the term Putting all together, we get the desired claim (2.11) where P is a polynomial of degree 3 in v. Now, we have to show that ǫ v 2 x (t) ∞ ∈ L ∞ (0, T ).Let us first exploit the boundary condition in (2.6).Multiplying it by v x and using (2.7), we get and thus we get that ǫv 2 x ≤ v 2 on ∂(0, 1).Moreover, using Step 2 we finally deduce that Taking φ = 1 in (2.12), we obtain that together with (2.13) and the maximum principle, imply that for some constant C that depends on the bound (2.10), and thus independent of ǫ.
Summarizing, now the term ) with bounds independent of ǫ.Then the argument given above shows that there are local Lipschitz bounds on v ǫ uniform in ǫ.
Step 6. Interior regularity of higher order derivatives uniform in ǫ.Thanks to the smoothness results stated in Step 4 and the local uniform bounds on the gradient in Step 5, the classical interior regularity results in [31, Chapter V, Theorem 3.1] shows uniform (in ǫ) interior bounds for any space and time derivative of v ǫ .
Step 7. Passing to the limit as ǫ → 0 + .Letting ǫ → 0 + is not completely obvious due to the boundary condition (2.4).Another difficulty stems from the fact that we do not know if v ǫt are Radon measures with uniform bounds in ǫ.This means that the notion of normal boundary trace has to be considered in a weak sense as considered in [2] (see also [3,Section 5.6] or [9]).Thus, we only sketch the proof of this result.Let us first prove that the interior regularity bounds on v ǫ permit to pass to the limit and obtain a solution v of x in D ′ ((0, T ) × (0, 1)) . Let Estimate (2.14) implies that a ǫ are uniformly bounded independently of ǫ.Then by extracting a subsequence, we may assume that a ǫ ⇀ a ∈ L ∞ ((0, T ) × (0, 1)) weakly * .On the other hand, the interior regularity bounds on v ǫ ensure that a = vη √ v 4 +(vx) 2 .By passing to the limit as ǫ → 0, we have ) by ϕ and integrate by parts, we obtain This is a weak form of the boundary condition (2.4).The correct notion of weak trace is much more technical and is described in [3].Using Lemma 5.7 in [9] one can directly obtain that v satisfies (2.4) in this generalized sense.Since we do not need this result here, we skip the details that would need several technical definitions to be fully explained.
Remark 2.2.Note that we can apply Step 5 to the smooth solution obtained in Theorem 2.1 to the Cauchy problem (2.3)- (2.4).In this case f ∞ ≤ P (v, φ, φ x )|φ x | ∞ and we obtain a local Lipschitz bound for v(t, x) which only depends on local uniform bounds of v(t, x) and on the local Lipschitz bound of v 0 (x).
Remark 2.3.In Section 2.2 we will give sufficient conditions on u 0 that imply that v t is a Radon measure.In that case, the notion of weak trace a • ν can be found in [21,23].
Remark 2.4.We could define the notion on entropy solutions of equation (2.3) with boundary condition (2.4) and prove that the solution constructed is indeed an entropy solution of it.We will not pursue this here.

Getting an entropy solution of (1.3) from (2.3)
Here, we use several notations and definitions that are introduced in the Appendix to which we refer for details.In this Section, we come back to the notation v(t, η) instead of v(t, x), η ∈ (0, 1).Recall that by passing to the limit as ǫ → 0 + we have found a solution v of for any T > 0. Thus, let v(t, η) be the solution of (2.15) constructed in Theorem 2.1 which satisfies [a(t, η) • ν] = 1 for η = 0, 1 and a.e. for t ∈ (0, T ) in a weak sense.As we shall see, we do not need this here, we only need a weaker form of the boundary condition as expressed in (2.17) below.
In the next Lemma we construct an entropy solution of (1.1) from a solution v(t, η) of (2.15).To prepare its statement, let , By (2.8), we have ) and any t > 0.
The statement (ii) in Theorem 1.1 follows from next Proposition.
Proof.(i) Since v is bounded and bounded away from zero from Step 2 in Theorem 2.1, then u is bounded and bounded away from zero in its support.The smoothness properties of v prove that u ∈ C([0, T ], L 1 (IR)), u(0) = u 0 , and u(t) is smooth inside its support.By Step 3 from Theorem 2.1, we have that u(t) ∈ W 1,1 (a − t, b + t) for almost any t ∈ (0, T ).This implies that u(t) ∈ BV (IR) for almost any t ∈ (0, T ).From the change of variables (2.16) we have that (ii) For simplicity, let us write Q T = (0, T ) × IR, and Ω(t) = (a − t, b + t).Since we have that u ∈ L 1 loc,w (0, T ; BV (IR)).We have denoted by u i (t) the trace of u| Ω(t) on ∂Ω(t).Note that it coincides with u + (t).Let us prove that where (2.18) was used.Thus (2.19) holds.
(iii) To prove that u is an entropy solution of (1. 3), we have to prove that holds for any any T, S ∈ T + and any φ ∈ D((0, T ) × IR), φ(t, x) = η(t)ρ(x).As in [6, Proposition 1], we have (h S (u(t), DT (u(t))) where R(r) = r, r ∈ IR.Thus, by (2.21) and (2.22), we get On the other hand, it is easy to prove that Adding (2.23) and (2.24), we obtain To simplify the subsequent notation let us denote p(u) = T (u)S(u) = J ′ (u) and J(u) = J T S (u).Let us now prove that The main technical difficulty comes from the fact that we do not know that u t = z x is a Radon measure.We circumvent this difficulty by using instead discrete derivatives.Let us denote Then, we can obtain which is a discrete version of (2.26).Note that we have used the inequality ∆ + τ u(t)p(u(t+τ )) ≥ ∆ + τ J(u) which is a consequence of the convexity of J.By letting τ → 0 + , we need to show that and This will result in (2.26).The limit (2.27) follows since u(t) ∈ BV (IR) a.e. in t, u ∈ L 1 w (0, T ; BV (IR)) (hence u x (t) ∈ L 1 (0, T )) and the trace functions u(t, a − t), u(t, b + t) are integrable in [0, T ].The second limit (2.28) follows easily.To prove (2.29), for any τ > 0 let ψ τ (t, x) := 1 τ t+τ t φ(s, x)p(u(s, x)) ds, Observe also that Then, as τ → 0 We have proved (2.29).Finally we observe that from (2.25) and (2.26) we obtain (2.20).
Remark 2.6.In a similar way, using this time and ∆ + τ (J(u))(t) ≥ p(u(t))∆ + τ (u)(t) one can prove that the opposite inequality in (2.26) holds, and we have equality.Note also that equality holds also in the entropy conditions (2.20).
With some additional regularity on the initial condition, one has that u t is a Radon measure in (0, T ) × IR.Indeed, the following proposition follows immediately from the results in [10,21].
If u is the entropy solution of (1.3) with initial data u 0 , then u t is a Radon measure in (0, T ) × IR.
From Proposition 2.7 and the results in [21], it follows that [z • ν Ω(t) ] = −u i (t) on ∂Ω(t) for almost any t ∈ (0, T ).This permits also to define the notion of normal trace of a(v, v η ) in the sense of [21,23].

Regularity for touching-down initial data
Let us start by getting local estimates.
Let v 0δ (η) be the functions obtained by the change of variables (2.2) (with t = 0).Let u δ (t, x) be the entropy solution of (1.1) with u δ (0) = u 0δ .By Theorem 1.1 we know that each u δ (t, x) is smooth inside (a, b).Let us note that the local bounds on u δ and its derivatives do not depend on δ.It suffices to observe that this is true for the associated functions v δ (t, η) which are solutions of (2.3), (2.4), with initial data v δ (0, η) = v 0δ (η).Note that the bounds in Steps 1, 2, 3 in the proof of Theorem 2.1 are independent of δ.By Remark 2.2, the Lipschitz bound in Step 5 depends only on the local Lipschitz bounds of v 0δ (η) and are, thus, uniform in δ.Step 6 proves uniform (in δ) interior bounds for any space and time derivative of v δ (t, η).By passing to the limit as δ → 0+ we conclude that u(t, x) is smooth inside its support and (i) and (ii) in Theorem 1.1 hold.
) with uniform local Lipschitz bounds in (a, b).Let v 0δ (η) be the functions obtained by the change of variables (2.2) (with t = 0).Let u δ (t, x) be the entropy solution of (1.1) with u δ (0) = u 0δ .By Theorem 1.1 we know that each u δ (t, x) is smooth inside (a, b).Let us note that the local bounds on u δ and its derivatives do not depend on δ.Again, it suffices to observe that this is true for the associated functions v δ (t, η) which are solutions of (2.3), (2.4), with initial data v δ (0, η) = v 0δ (η).
The L p bounds follow from Step 1 in the proof of Theorem 2.1 for p ∈ [1, ∞) and they only depend on the L p bound of v 0δ .Actually, we have that depends on the integrability of 1 u 0δ (x) at the boundary points.But multiplying (2.3) by v p δ φ where p ∈ [1, ∞) and φ is a positive smooth test function with compact support in (0, 1) we obtain Thus we derive local L p bounds for v δ which are independent of δ.We also obtain local bounds on the total variation of v p δ which are independent of δ.To obtain a local L ∞ bound independent of δ we observe that this follows from the identity v δ (t, η) = 1 u δ (t,x) , where x = ϕ δ (t, η) is given by (2.2), since we know that u δ (t, x) is locally bounded away from zero in its support [6].Thus Steps 1, 2, 3 hold in their local versions.By Remark 2.2, the Lipschitz bound in Step 5 depends only on the uniform local bounds on v δ (t, η) and on the local Lipschitz bounds of v 0δ (η) and are, thus, uniform in δ.Step 6 proves uniform (in δ) interior bounds for any space and time derivative of v δ (t, η).By passing to the limit as δ → 0+ we conclude that u(t, x) is smooth inside its support and (i) and (ii) in Theorem 1.1 hold.The last assertion is a consequence of the comparison principle using Lemma 3.4 below.Remark 3.3.Note that the last assertion implies that if the initial profile is not vertical at the boundary at t = 0 it remains non-vertical for any t > 0.Moreover, during the proof we have observed that if u 0 has a vertical profile with 1 u 0 ∈ L p (a, b), then 1 u(t,x) ∈ L p (a − t, b + t) for any t > 0. Thus in that case u(t, x) has a vertical profile at the boundary of its support.
Due to translational invariance of (1.3), we state our next Lemma in an interval symmetric around zero.
Proof.Computing the derivatives, we get

and
( where , and then Thus, the claim x holds if and only if Indeed, the above inequality is implied by 2R ≥ 4αx 2 /Q and 4αx 2 ≤ (2R)2α|x| ≤ 2RQ.Now, we choose A such that that is, ) is a supersolution of (1.3).

Numerical experiments and heuristics
In this section, we will propose a numerical scheme for more general equations than the RHE (1.1).We deal with the Cauchy problem for the generic porous media relativistic heat equation (RHEm) [22] given by with initial data u 0 a probability density with compact support.In order to propose the numerical scheme, we make use of the change of variables to Lagrangian coordinates.As in the introduction, let us denote by F the distribution function associated to the probability density u and ϕ(t, η) its inverse or generalized inverse, defined by Here, we have preferred to shift the mass variable to the interval (− 1 2 , 1 2 ) to simplify the notations about boundary conditions.In this way, we simply have the relation For simplicity, most of the numerical tests have been chosen for even initial data.Observe that this change of variables is a weak diffeomorphism in case of connected compactly supported smooth u, say on the interval (−A(t), A(t)) in which case lim Straightforward computations show that the equation satisfied by ϕ in (− 1 2 , 1 2 ) is while at the boundary, formally, by (4.2) and (4.4), we have to impose Moreover, thanks to the vertical contact angle property (see (2.1) for the RHE and [22] for the RHEm), we have that lim The purpose of this section is two-fold.On one hand, we heuristically observe some qualitative properties from the Lagrangian viewpoint.On the other hand, these properties are confirmed by numerical experiments with the use of an adaptation of the algorithm proposed in [16] for general equations in continuity form for the 2-dimensional case.

Numerical Method
Equations (1.3) and (4.1) have been numerically treated in [34,38] using the connection between nonlinear diffusions and Hamilton-Jacobi equations and numerical methods for conservation laws and in [10] using an appropiate WENO scheme.Here, we propose a completely different approach based on the optimal transportation viewpoint.As we already mentioned in the introduction an explicit Euler discretization of the equation satisfied by the generalized inverse (4.5) coincides with the variational scheme introduced in [29,33].Moreover, the theoretical result proven in [35] shows that this scheme applied to (1.3) is convergent for initial data compactly supported smooth in their support and bounded below and above.Therefore, we plan to use a similar algorithm for Eq.(4.1).This Lagrangian formulation in 1D for nonlocal and nonlinear diffusion problems was numerically analysed in [28,13].These Lagrangian coordinates ideas were generalized to several dimensions in [16].
The advantages of this method are the adaptation of the mesh to the mass distribution of the solution in an automatic way, the immediate positivity of the solutions, and the decay of the natural Liapunov functional of the equations.We refer to [16] for more details and discussions on these issues.
Here, we propose an adaptation of the algorithm in [16].First of all, the discretization in the mass variable has been treated by finite difference approximations of the derivatives of the unknown ϕ.We consider a partition {η i } i=1:N of the spatial interval [− 1 2 , 1 2 ] and we let ∆ i := η i+1 − η i .Note that, due to (4.2), first derivatives at the points corresponding to the nodes η 2 and η N −1 have to be taken from the inside of the domain.In order to avoid higher errors in the approximation of the derivative at the boundaries, we decide to approximate ϕ η as with η(t) to be specified.The derivative of the term 1 ϕη is computed in the other direction for better stability properties of the approximation of (ϕ η ) −1 η .At the boundary we just impose (4.6).As explained in [16], the point η(t) has to be taken as the global maximum for u, which can be tracked at any time step.In all examples computed, initial data are taken to be radially symmetric and decreasing from the point x = η = 0.In all of them, the global maximum stays at x = η = 0. Therefore, we choose to take an even number of points N in the discretization and to take a symmetric partition Let us point out that the spatial partition is never uniform since the change to Lagrangian coordinates produces the accumulation of nodes near the global maximum.We instead want to follow some particular features of these type of equations such as propagation of fronts with a vertical contact angle or formation of singularities.Therefore, the partitions will be chosen accordingly in order to accumulate more points around the points ± 1 2 and other points of interest.The time derivative is evaluated through a simple explicit Euler scheme with the CFL condition proposed in [16]; i.e: with α CF L > 2, for the porous-medium equation which is the large-time limit behaviour of (4.1), see [22] and subsection 5.3.All our simulations are done with α CF L = 8.Although the CFL analysis in [16] applies only to equations written in variational form that includes (4.1)only for m = 1, all numerical tests seem not to be affected by the chosen CFL condition.Finally, we point out that u(t, ϕ(t, η 1 )) = u(t, ϕ(t, η N )) = 0.Because of this fact, in all the plots which follow, the first and last nodes are never plotted.
We define next w(t, η) := ψ(t, η)ψ η (t, η) = u x (t, ϕ(t, η)).The analysis above also shows that, in case ψ(t 0 , ± 1 2 ∓ ) = 0, then |w(t, ± 1  2 )| ≤ |w(0, ± 1 2 )|.On the other hand, in the bulk, w verifies the following equation Thus, if w 0 is initially bounded, w remains bounded in [− 1 2 , 1  2 ] as proved in Section 3. Observe that at a point η 0 of maximum of w, we have (η 0 ) ≤ 0, implying the claim.We show a numerical experiment with u 0 (x) = (1 − |x|) + as initial datum which does not satisfy the conditions of Theorem 1.1.We take N = 1000 for the simulations.We point out that since the initial datum is 0 at the extremes of the support, we need a lot of nodes in the discretization near them since due to the change of variables (4.3), then ϕ η (± 1 2 ∓ ) = ∓∞ and we want the numerical scheme to be able to capture this feature.We report in Fig. 1 the precise evolution of the support showing the smoothing effect at x = 0, the boundedness of the derivative all over the support including the boundaries, and the expansion of the boundary at precise unit speed as expected by the theory in Theorem 1.1 and the heuristic arguments above.Let us now take m > 1.In case u(±A(0) ∓ ) = 0 (i.e.ψ(0, ± 1 2 ∓ ) = 0), then (4.8) implies that the support of the solution does not move at all whenever u(±A(t)) = 0.The solution will become positive at the tip of the support u(±A(t)) > 0 with t > t 0 > 0 if and only if ψ t (t 0 , ± 1 2 ∓ ) ∈ (0, +∞] with u(±A(t 0 ) ∓ ) = ψ(t 0 , ± 1 2 ∓ ) = 0.In case u x (±A(t 0 )) = 0, we can use (4.11) to approximate terms 2 in the expression of (4.9) to get As a consequence, ψ t (t 0 , ± (u m+1 (x)) x (t 0 , ±A(t)) > 0 .
In such case, the solution becomes positive at ±A(t) and then, according to (4.8), its support starts to increase.We note that this waiting time phenomenon is similar to that of the classical porous medium equation but the condition for the support to start moving is completely different to the one obtained  in [11].Supposing a potential growth of ψ, i.e. ψ(t 0 , η) we obtain that ψ t (t 0 ± 1 2 ∓ ) = +∞ if and only if p < 1 m+1 .We point out that this behavior has already been numerically obtained in [10].In Fig. 2, we show this waiting time phenomenon for m = 1.5.One can observe that initially the support does not move since the behavior near the boundary is ψ(0, η) ≃ C 0, 1  2 − |η| 1 2 , then the derivative at the boundary builds up until the behavior at the boundary reaches the critical value producing the lift-off of the boundary point.More interesting is the case m = 3 which we show in Fig. 3. There, a discontinuity in the bulk appears before the support starts to move.

Formation of discontinuities in the bulk
In view of the first example in the last section, one may think that discontinuities may appear only as a consequence of the waiting time phenomenon; i.e. particles tend to dissipate but their support does not move, which may create the discontinuities.In this section we heuristically study that it is possible to create discontinuities inside the bulk even if the solution are far away from zero as seen in Fig. 3.
First we treat the case m = 1.In case of an upwards jump discontinuity or a vertical angle at a point η 0 ∈] − 1 2 , 1 2 [ such that ψ η (η 0 ) ± = +∞ , then we also have ϕ t (η 0 ) = −1.Since |ϕ t | ≤ 1, then ϕ t (η 0 ) = −1 implies that ϕ t is nonincreasing to the left and nondecreasing to the right of η 0 , i.e., ((ϕ η ) − ) t ≤ 0 and ((ϕ η ) + ) t ≥ 0. This shows that (ψ(η 0 ) − ) t ≥ 0 while (ψ(η 0 ) + ) t ≤ 0, which implies that the size of the discontinuity reduces in for an upwards jump discontinuity or that no discontinuity is created if initially there is a vertical angle.This last phenomenon is not true if m > 1 in the case of a vertical angle at a point η 0 ∈] − 1 2 , 1 2 [ such that ψ η (η 0 ) ± = +∞.From the equation (4.9) for ψ as in previous subsection, we deduce that ψ t (η 0 ) = (m − 1)ψ m ψ η (η 0 ), and thus, a discontinuity is created.Once we have a discontinuity at η 0 the evolution is theoretically unknown.In order to show this behavior we have taken two types of initial datum with N = 1000: We imposed a high concentration of nodes around the vertical angles or discontinuities (i.e.x = ± 1 2 ).In Fig. 4 we observe the evolution of the solutions corresponding to the initial datum u 0 , demonstrating the above heuristics.In Fig. 5, we see how an initially discontinuous initial datum ũ0 is smoothed during the evolution both for m = 1 (as heuristically deduced before) and for m = 2.We observe that the smoothing of the discontinuity is slower with m > 1 than that of m = 1.

Asymptotic behavior
In this Section, guided by heuristics, we numerically observe the asymptotic behavior of solutions to (4.1) and the rate of convergence towards their asymptotic steady state, for which no result is available in the literature.Performing the classical self-similar change of variables [17] that translates porous medium equation onto nonlinear Fokker-Planck equations given by v(x, t) = e t u(e t x, k(e with k = 1 m+1 , then equation (4.1) transforms into Therefore, formally, when t → ∞ solutions of (4.14) should converge to a stationary solution of v t = div xv + v m−1 ∇v , i.e., to a Gaussian V (x) for m = 1 or to the corresponding Barenblatt solution V m (x) when m > 1 given by where C m is uniquely determined by the conservation of mass.In the original variables, then solutions should converge to the corresponding self-similar profiles obtained from V and V m via the change of variables (4.13) except time translations.To be precise, the self-similar solutions are given by where Cm is determined as above.In the following computations we have taken u ] , N = 100.We plot the evolution of the initial datum for different values of m and an estimate of the difference of u − U m in the L 1 -norm.More precisely, we took u Some comments are in order.First of all, in Figure 6, we note that for m = 1, while time is small, the numerical solution satisfies both the linear propagation of the support property, as well as the vertical contact angle property.However, for larger times, these two conditions are lost during the computation.This is due to the fact that we took a fixed number of nodes (N = 100), and as time increases, this number of nodes is clearly insufficient.We have observed that by increasing the number of nodes (for instance to N = 1000) the time in which the numerical solution is more accurate increases.We can also see in Figure 6, that, in spite of this, the numerical solution tends to a Gaussian with an algebraic rate of convergence that seems to be 1  2 the one of the heat equation.However, it is exactly by the same reason as before that when time increases, the rate of convergence degenerates.For this reason, we have included in Figure 6 the L 1 -convergence rate with N = 1000.Instead, when m > 1, the support of the solution does not propagate so fast and we can observe in Figures 7 and 8 how the vertical contact angle property is preserved even for large times.Moreover, in Figures 7 and 8 we can see how the numerical solution tends to U m for m = 2 and m = 10.In both cases, the rate of convergence is algebraic and, numerically, it is surprisingly seen that it might correspond to 1  3 in the first case and to 1  11 in the second one; i.e.: the same convergence rate as for the porous medium equation, see [17,33]   when ν → ∞.If φ ∈ C c (IR N ), we write φ = φ + − φ − with φ + = max(φ, 0), φ − = − min(φ, 0), and we define g(u, DT (u)), φ := g(u, DT (u)), φ + − g(u, DT (u)), φ − .
Recall that, if g(z, ζ) is continuous in (z, ζ), convex in ζ for any z ∈ IR, and φ ∈ C 1 (IR N ) + has compact support, then g(u, DT (u)), φ is lower semi-continuous in T BV + (IR N ) with respect to L 1 (IR N )-convergence [27].This property is used to prove existence of solutions of (A.1).

Figure 2 :
Figure 2: Evolution of u 0 in case m = 1.5 at different times.Left: Before the discontinuity at the tip of the support appears.Right: Evolution of the discontinuity front after.

Figure 3 :
Figure 3: Evolution of u 0 in case m = 3 at different times.Top left: Before a discontinuity on the bulk appears.Top right: After the discontinuity front forms till it reaches the tip of the support.Bottom: After the discontinuity front starts to move.

Figure 4 :
Figure 4: Evolution of solutions corresponding to u 0 .Top left: Initial datum u 0 .Top right: Evolution for m = 1 at small times.Bottom left: Evolution for m = 1 for larger times.Bottom right: Evolution for m = 4.

Figure 5 :
Figure 5: Evolution of solutions for ũ0 .Top: Evolution for m = 1 at different times.Bottom: Evolution for m = 2 at different times.

Figure 6 :
Figure 6: Top and Bottom left: Evolution of u 0 in case m = 1 at different times with N = 100.Bottom right: log-log plot of the estimate u(t) − U (t) 1 with N = 1000.

Figure 7 :
Figure 7: Top and Bottom left: Evolution of u 0 in case m = 2 at different times.Bottom right: log-log plot of the estimate u(t) − U 2 (t) 1 . .

4. 4
Convergence toward the homogeneous relativistic heat equationWe finally show numerically how solutions to (1.1), converge to solutions of the homogeneous relativistic heat equation   u t = u u x |u x | x in IR N × [0, T ] u 0 (x) = u 0 in IR N

Figure 8 :
Figure 8: Top: Evolution of u 0 in case m = 10 at different times.Bottom left: log-log plot of the estimate u(t) − U 10 (t) 1 .Bottom right: zoom of the final time intervalwhen the kinematic viscosity ν → +∞ as already proved in[7].In Fig.9we estimate the evolution in time of the difference in the L 1 -norm for solutions corresponding to the initial data u0 = χ [− 1 2 , 1 2 ]for different values of ν with respect to the explicit solution u hom , given byu hom (x, t) = 1 1 + 2t χ [− 1 2 −t, 1 2 +t]

Figure 9 :
Figure 9: Top left: Numerical solution at t = 1 for different values of ν.Top Right: Numerical solution at t = 100 for different values of ν.Bottom: Evolution of the L 1 -difference with respect to u hom .