by Lemma 5.4 it satisfies the Oleĭnik E-condition (4.17) at any discontinuity
of $M(t,x)$ for a.e. $t$, and so by (5.2) and Lemma 5.2 we conclude that
$U(t)-V(t)$ is contained in the normal cone $N_{X(t)}\mathscr{K}$ for a.e.
$t$. Then it follows that $V(t)=\operatorname{P}_{T_{X(t)}}U(t)$ for a.e. $t$.
Therefore $X(t)$ solves the differential inclusion (3.21a) and is almost
everywhere differentiable again due to (4.24). Then, as argued in the proof of
[7, Theorem 3.5], see also [10, Remark 3.9], at every point of
differentiability we must have
$\dot{X}(t,m)=V(t,m)=\operatorname{P}_{\mathscr{H}_{X(t)}}U(t)$. It remains to
show that $\lim\limits_{t\to 0+}X(t)=X^{0}$ in $L^{2}(\Omega)$, but this
follows from the identity (4.20) and the time-continuity (4.24) for $p=2$. ∎

######  Remark 5.8.

In the proof of Theorem 5.6 we used the Kružkov, or integral, formulation to
show that we obtain an entropy solution. Of course, we could have
alternatively used Lemma 5.4 and only verified the Oleĭnik E-condition at the
discontinuities.

######  Remark 5.9 (Connection to a uniqueness result for a related system).

As mentioned in the introduction, another Euler-type system on the line is
studied in [32], namely

$\partial_{t}\rho+\partial_{x}(\rho f^{\prime}(v))=0,\qquad\partial_{t}(\rho
v)+\partial_{x}(\rho vf^{\prime}(v))=\rho
h\left(\int_{-\infty}^{x}\rho\mathop{}\\!\mathrm{d}y\right),$

where $h\colon\mathbb{R}\to\mathbb{R}$ is a given Lipschitz function, and
$f\colon\mathbb{R}\to\mathbb{R}$ is a smooth, convex function satisfying
$\lim_{\lvert v\rvert\to+\infty}f(v)/\lvert v\rvert=+\infty$. In particular,
for $f(v)=\frac{1}{2}v^{2}$ we recover a version of (1.1). In their study they
make use of a generalization of the Hopf–Lax–Oleĭnik formula, which is similar
to the sticky Lagrangian solutions in Section 3 in the sense that there is a
representation formula, see also [46]. An interesting observation is that
their uniqueness results require $h$ to be nondecreasing, which means that its
primitive is concave. If one, to make sense of the above system in the case of
measures with Dirac parts, replaced
$h(\int_{-\infty}^{x}\rho\mathop{}\\!\mathrm{d}y)$ with its Vol’pert average
(4.7), the Lagrangian representation satisfying (3.16) would be
$F[X](m)=h(m)$, independent of $X$ as in the Euler–Poisson case. Furthermore,
since $h$ is nondecreasing, using Lemma 3.1 one can then verify that this
functional is sticking in the sense of Definition 3.5, meaning the unique
Lagrangian solutions are globally sticky as in Definition 3.9.

######  Remark 5.10 (Connection to a uniqueness result for the Euler–Poisson
system).

Note that in the particular case of (1.1) with $f[\rho]\equiv 0$, the “other”
Oleĭnik condition (1.7) has been established independently for both entropy
solutions [36] and Lagrangian solutions [35]. As mentioned in the
introduction, [26] establishes uniqueness of entropy solutions for the
attractive Euler–Poisson system (1.3) with $\alpha=1$, $\beta=0$, and initial
data $\rho_{0}\in\mathscr{P}_{\mathrm{c}}(\mathbb{R})$ and $v_{0}\in
L^{\infty}(\mathbb{R},\rho_{0})$.

Observe that
$\mathscr{P}_{\mathrm{c}}(\mathbb{R})\subset\mathscr{P}_{2}(\mathbb{R})$,
while $v_{0}\in L^{\infty}(\mathbb{R},\rho_{0})$ for
$\rho_{0}\in\mathscr{P}_{\mathrm{c}}(\mathbb{R})$ also belongs to
$L^{2}(\mathbb{R},\rho_{0})$, such that Lagrangian solutions also cover such
initial data. Of course, showing that $\rho$ and $v$ remain in these more
restrictive spaces requires more analysis, but should rely on the finite speed
of propagation due to bounded velocities, or equivalently, the flux function
(4.9) here being Lipschitz.

Compared to our result, we merely argue that the Lagrangian solutions of the
attractive Euler–Poisson system studied here are also entropy solutions in the
sense of [26]. Indeed, the weak continuity of the kinetic energy follows from
the right-continuity of $V$, cf. Lemma 3.7, which then carries over to $\lVert
V\rVert^{2}_{L^{2}(\Omega)}$. It remains to show (1.7), and to that end, as
noted in [35], it is sufficient that the map $m\mapsto X(t,m)-tV(t,m)$ is
nondecreasing. Using the fact that (1.3) for $\alpha\geq 0$ has globally
sticky solutions given by the projection formulas (3.22) and (3.23), together
with (3.13) we find
$X(t)-tV(t)=\operatorname{P}_{\mathscr{H}_{X(t)}}(Y(t)-tU(t))$. Since
$\operatorname{P}_{\mathscr{H}_{X}}$ preserves monotonicity, it is sufficient
to show that $Y-tU$ is nondecreasing, and we easily compute
$Y(t)-tU(t)=X^{0}(m)+\frac{1}{2}t^{2}\alpha(m-\frac{1}{2})$ which in our
attractive case $\alpha>0$ clearly is increasing.

##  6\. Particle dynamics for the Euler–Poisson equation

To illustrate the equivalence between the Lagrangian and entropy solutions, we
will in this section return to the particle dynamics for the Euler–Poisson
system (1.3). The motivation for considering particle dynamics comes from its
central role in proving the equivalence of entropy solutions and gradient
flows in [4]. In addition, limits of particle dynamics play key roles in both
the gradient flow formulation of [7] and the scalar conservation law-point of
view in [36, 37].

Particle solutions correspond to piecewise constant distribution functions of
the corresponding conservation law, and this is where we will make use of the
connection between the $L^{2}$-projections and shock admissibility conditions
from the previous section. Indeed, a piecewise constant function will be a
piecewise smooth solution of the conservation law, for which the Kružkov
entropy condition (4.14) is equivalent to the Oleĭnik E-condition (4.17) for
the shocks.

The benefit of considering the specific case of the Euler–Poisson system is,
as seen in Section 3.4, that the forcing term $F[X](m)$ turns out to be
independent of the position $X$. In particular, the prescribed velocity $U$
can be computed for all times directly from the initial data, and will be a
linear function of $t$ since trajectories undergo constant acceleration
between collisions. It is also of interest to study the distinction between
the repulsive ($\alpha<0$) and non-repulsive ($\alpha\geq 0$) cases, since the
global behavior of solutions is quite different in the two cases, as
emphasized in [7].

Although the particle dynamics and their connection to the conservation law
(4.4) for $\alpha=0$ and $\alpha\geq 0$ are respectively treated in [8] and
[36], we recall the key points here for completeness and ease of comparison.

###  6.1. Discrete setup

Recalling the setting of Section 2, we have $n$ initially distinct particles
located at positions $x_{i}$ with velocities $v_{i}$ and masses $m_{i}$. From
their masses we can define a partition of $\Omega=(0,1)$ given by
$0=\theta_{0}<\theta_{1}<\dots<\theta_{n-1}<\theta_{n}=1$ where

(6.1) $\theta_{0}=0,\qquad\theta_{i}=\sum_{j=1}^{i}m_{j},\quad
i\in\\{1,\dots,n\\}.$

Based on (3.5) we can define the functions

(6.2)
$X(m)\coloneqq\sum_{i=1}^{n}x_{i}\chi_{[\theta_{i-1},\theta_{i})}(m),\qquad
V(m)\coloneqq\sum_{i=1}^{n}v_{i}\chi_{[\theta_{i-1},\theta_{i})}(m),$

the finite-dimensional Hilbert space

$\mathscr{H}_{\mathbf{m}}\coloneqq\left\\{X(m)=\sum_{i=1}^{n}x_{i}\chi_{[\theta_{i-1},\theta_{i})}(m)\colon\mathbf{x}=(x_{1},\dots,x_{n})\in\mathbb{R}^{n}\right\\}\subset
L^{2}(\Omega),$

and the convex cone

$\mathscr{K}_{\mathbf{m}}\coloneqq\left\\{X(m)=\sum_{i=1}^{n}x_{i}\chi_{[\theta_{i-1},\theta_{i})}(m)\colon\mathbf{x}=(x_{1},\dots,x_{n})\in\mathbb{K}^{n}\right\\}\subset\mathscr{K}\subset
L^{2}(\Omega),$

where we recall $\mathbb{K}^{n}$ from (2.4), so that
$X\in\mathscr{K}_{\mathbf{m}}\subset\mathscr{K}$ and $V\in\mathscr{H}_{X}$.
Here $X$ and $V$ provide us with particle solutions as in (2.1), since

$\rho=X_{\\#}\mathfrak{m}=\sum_{i=1}^{n}m_{i}\delta_{x_{i}},\qquad V=v\circ
X,\qquad\rho v=\sum_{i=1}^{n}m_{i}v_{i}\delta_{x_{i}}.$

Letting $H(x)$ be the right-continuous Heaviside function, that is,

$H(x)=\begin{cases}1,&x\geq 0\\\ 0,&x<0,\end{cases}$

we can integrate the measures above to obtain discrete counterparts of (4.1),
namely

$M(x)=\sum_{i=1}^{n}m_{i}H(x-x_{i}),\qquad
Q(x)=\sum_{i=1}^{n}m_{i}v_{i}H(x-x_{i}).$

Note that $M$ above only takes values in the set $\\{\theta_{i}\\}_{i=1}^{n}$
defined in (6.1). We recall that in order to define the scalar conservation
law (4.4) we needed to find $\mathcal{V}^{0}$ such that
$Q=\mathcal{V}^{0}(M)$, which turned out to be the primitive of $V$. For the
“acceleration”-part $\mathcal{A}(M)$ we need to compute the square of $M$, and
to do this we conveniently rewrite it as a sum of step functions with disjoint
supports,

$M(x)=\sum_{i=1}^{n}\theta_{i}\chi_{[x_{i},x_{i+1})}(x),$

where we for further convenience define $x_{n+1}\equiv\infty$. Then it is
clear that

$M(x)^{2}=\sum_{i=1}^{n}\theta_{i}^{2}\chi_{[x_{i},x_{i+1})}(x)=\sum_{i=1}^{n}(\theta_{i}^{2}-\theta_{i-1}^{2})H(x-x_{i})=\sum_{i=1}^{n}(\theta_{i}+\theta_{i-1})m_{i}H(x-x_{i}).$

This shows that $-\frac{\alpha}{2}(M(t,x)^{2}-M(t,x))$ can be written as
$\int_{0}^{M(t,x)}A_{n}(m)\mathop{}\\!\mathrm{d}m$ for the piecewise constant
function

(6.3) $A_{n}(m)=\sum_{i=1}^{n}a_{i}\chi_{[\theta_{i-1},\theta_{i})}(m),\qquad
a_{i}\coloneqq-\frac{\alpha}{2}\left(\theta_{i}+\theta_{i-1}-1\right).$

Note that
$\mathcal{A}_{n}(m)=\int_{0}^{m}A_{n}(\omega)\mathop{}\\!\mathrm{d}\omega$ is
exactly the piecewise linear interpolation of the quadratic function
$\mathcal{A}(m)=-\frac{\alpha}{2}(m^{2}-m)$ on $[0,1]$ with breakpoints
$\\{\theta_{i}\\}_{i=1}^{n}$. Moreover, (6.3) is exactly the
$L^{2}$-projection of $F[X](m)=-\frac{\alpha}{2}(2m-1)$ onto
$\mathscr{H}_{\mathbf{m}}$ as defined in [7, Section 5], namely

$F_{\mathbf{m}}[X]\coloneqq\operatorname{P}_{\mathscr{H}_{\mathbf{m}}}F[X]=\sum_{i=1}^{n}a_{\mathbf{m},i}\chi_{[\theta_{i-1},\theta_{i})}(m),\quad
a_{\mathbf{m},i}=\frac{1}{m_{i}}\int_{\theta_{i-1}}^{\theta_{i}}F[X](m)\mathop{}\\!\mathrm{d}m=-\frac{\alpha}{2}\left(\theta_{i}+\theta_{i-1}-1\right),$

where we recognize $a_{\mathbf{m},i}=A_{n}(\theta_{i-1})=a_{i}$, the
acceleration for our particles in (6.3). This projected force satisfies
(3.16); indeed, recall from the derivation of $\mathcal{A}$ in the continuous
flux function (4.9) that $f[\rho]=-\frac{\alpha}{2}(M(x)+M(x-)-1)\rho$, and so
we compute using (6.2) and (6.3)

$\displaystyle\int_{\mathbb{R}}\varphi(x)f[\rho](\mathop{}\\!\mathrm{d}x)$
$\displaystyle=-\frac{\alpha}{2}\int_{\mathbb{R}}\varphi(x)(M(x)+M(x-)-1)\rho\mathop{}\\!\mathrm{d}x$
$\displaystyle=-\frac{\alpha}{2}\sum_{i=1}^{n}\varphi(x_{i})(\theta_{i}+\theta_{i-1}-1)(\theta_{i}-\theta_{i-1})$
$\displaystyle=\sum_{i=1}^{n}\int_{\theta_{i-1}}^{\theta_{i}}\varphi(X(m))A_{n}(m)\mathop{}\\!\mathrm{d}m=\int_{0}^{1}\varphi(X(m))A_{n}(m)\mathop{}\\!\mathrm{d}m.$

To summarize, from either point of view the acceleration of the $i$th particle
is $a_{i}$ before collision, and we have justified the initial dynamics (2.2),
which can be restated as the system

$\dot{x}_{i}=v_{i},\qquad\dot{v}_{i}=-\frac{\alpha}{2}\sum_{j=1,j\neq
i}^{n}m_{j}\operatorname*{sgn}(x_{i}-x_{j}),$

being the particle analogue of (1.3).

At this stage, with the initial particle evolution defined, we still need to
derive the corresponding conservation law (4.4) for this discrete setting.
From now on, to emphasize that we are considering particle solutions, we will
add the subscript $n$ to both the Lagrangian variables $X$ and $V$, the
cumulative quantities $M$ and $Q$, and the primitive $\mathcal{V}^{0}$.
Moreover, as the position $t\mapsto x_{i}(t)$ and velocity $t\mapsto v_{i}(t)$
of the $i$th particle in their definitions depend on time, we have that
$X_{n}(t,m)$, $V_{n}(t,m)$, $M_{n}(t,m)$, $Q_{n}(t,m)$ depend on both time $t$
and the material coordinate $m$.

Recall from Section 2 that the prescribed velocities
$u_{i}(t)=v_{i}^{0}+ta_{i}$ played a key role in defining the dynamics at
collision time, and so we define

(6.4)
$U_{n}(t,m)\coloneqq\sum_{i=1}^{n}u_{i}(t)\chi_{[\theta_{i-1},\theta_{i})}(m)=V_{n}^{0}(m)+tA_{n}(m)$

inspired by the definition of $V_{n}$. This is also the discrete counterpart
of (3.29), and from the same considerations as in Section 4 we find that the
flux function for the scalar conservation law is the primitive
$\mathcal{U}_{n}$ of (6.4), i.e.,

(6.5)
$\mathcal{U}_{n}(t,m)=\int_{0}^{m}U_{n}(t,\omega)\mathop{}\\!\mathrm{d}\omega=\mathcal{V}_{n}^{0}(m)+t\mathcal{A}_{n}(m),$

which is exactly the piecewise linear and continuous interpolation of the flux
function (4.9) with breakpoints $\\{\theta_{i}\\}_{i=1}^{n}$. In summary, the
scalar conservation law for the particle dynamics is

(6.6) $\partial_{t}M_{n}+\partial_{x}\mathcal{U}_{n}(t,M_{n})=0,$

with flux function $\mathcal{U}_{n}$ defined in (6.5).

From the point of view of Lagrangian solutions, we have from our definition of
the prescribed velocity that $\dot{U}_{n}=F_{\mathbf{m}}[X_{n}]$, and so, for
$X_{n}$ to be a Lagrangian solution we require $U_{n}-\dot{X}_{n}\in
N_{X_{n}}\mathscr{K}_{\mathbf{m}}$.

###  6.2. Dynamics prior to first collision

Here we show that the initial particle dynamics correspond to both the
Lagrangian solutions and the entropy solutions of (6.6).

Before any collision has happened, the particles $\mathbf{x}$ are located
within the interior of the cone $\mathbb{K}^{n}$, and consequently the
projection of the prescribed velocity onto the tangent cone is simply the
identity operator, that is,
$\operatorname{P}_{T_{\mathbf{x}}\mathbb{K}^{n}}\mathbf{u}=\mathbf{u}$.
Therefore the initial particle trajectories $\mathbf{x}(t)$ are uniquely
determined as the “free” trajectories
$\mathbf{y}(t)=\mathbf{x}^{0}+t\mathbf{v}^{0}+\frac{1}{2}t^{2}\mathbf{a}$. In
turn, $X_{n}$ is uniquely defined as $X_{n}=Y_{n}$, where

$Y_{n}(t,m)\coloneqq\sum_{i=1}^{n}y_{i}(t)\chi_{[\theta_{i-1},\theta_{i})}(m),\qquad
y_{i}(t)=x_{i}^{0}+tv_{i}^{0}+\frac{1}{2}t^{2}a_{i}.$

Note that $A_{n}=F_{\mathbf{m}}[X]$ is trivially Lipschitz-continuous as a
function
$F_{\mathbf{m}}[X]\colon\mathscr{K}_{\mathbf{m}}\to\mathscr{H}_{\mathbf{m}}$
since it does not depend on $X_{n}$, and so $(X_{n},V_{n})=(Y_{n},U_{n})$ is
the unique Lagrangian solution from Theorem 3.8.

From the point of view of conservation laws, we see that the particles,
through their representation by the piecewise constant mass distribution
$M_{n}$ and the piecewise linear and continuous flux function
$\mathcal{U}_{n}$, correspond exactly to Dafermos’ [18], or the front tracking
[9, 25], method for scalar conservation laws, at least initially. This
correspondence is exactly what was used to provide “sticky particle” solutions
of the Euler system in [8]. Since $M_{n}(t,x)$ is piecewise constant, it will
be a weak solution of (6.6), the discrete counterpart of the scalar
conservation law (4.4), provided it satisfies the Rankine–Hugoniot condition
(4.16) along the discontinuities $x_{i}(t)$. Moreover, it will be an entropy
solution if it additionally satisfies the Oleĭnik E-condition (4.17) along the
discontinuities, and in particular, it will be the unique $BV$ solution of
this conservation law for $t\in[0,T]$. This is indeed the case, since along
the $i$th trajectory we have a left and right state, $\theta_{i-1}$ and
$\theta_{i}$ respectively, such that $[[M_{n}]]=\theta_{i}-\theta_{i-1}=m_{i}$
and $[[\mathcal{U}_{n}(t,M_{n})]]=m_{i}(v_{i}^{0}+ta_{i})=m_{i}u_{i}(t)$.
Therefore, the Rankine–Hugoniot condition (4.16) is satisfied since

(6.7)
$\dot{x}_{i}(t)=u_{i}(t)=\frac{[[\mathcal{U}_{n}(t,M_{n}(t,x_{i}(t)))]]}{[[M_{n}(t,x_{i}(t))]]}=\frac{\mathcal{U}_{n}(t,\theta_{i})-\mathcal{U}_{n}(t,\theta_{i-1})}{\theta_{i}-\theta_{i-1}},$

and because $M_{n}$ cannot take any other values between $\theta_{i-1}$ and
$\theta_{i}$, the Oleĭnik E-condition (4.17) is trivially satisfied. Note that
since $U_{n}-\dot{X}_{n}=0\in N_{X_{n}}\mathscr{K}_{\mathbf{m}}$, this is also
the unique Lagrangian solution.

###  6.3. Collisions in the non-repulsive and repulsive case

Now that we have seen that the two notions of solution are equivalent for the
initial particle dynamics, it remains to show that this is still true when
particles meet. We shall treat the non-repulsive ($\alpha\geq 0$) and
repulsive ($\alpha<0$) cases separately, as there are major differences in the
collision dynamics in these two cases; in fact, as emphasized in [7],
$\alpha\geq 0$ corresponds to a sticking functional $F[X]$, while for
$\alpha<0$ the functional does not satisfy the sticking condition of
Definition 3.5.

To begin with we will establish some convenient notation and terminology.
Since the particles are ordered according to their positions, we see that the
particles on the left must necessarily have velocities larger than or equal to
the particles on the right before touching, or else they would not have met.
To be precise, we define the quantities

$\mathcal{J}_{i}(\tau)\coloneqq\\{j\in\\{1,\dots,n\\}\colon
x_{j}(\tau)=x_{i}(\tau)\\},\quad
i_{*}(\tau)\coloneqq\min\mathcal{J}_{i}(\tau),\quad
i^{*}(\tau)\coloneqq\max\mathcal{J}_{i}(\tau).$

That is, $\mathcal{J}_{i}(\tau)$ is the set indices of particles located at
$x=x_{i}(\tau)$ at time $t=\tau$, while $i_{*}(\tau)$ and $i^{*}(\tau)$ are
respectively the smallest and greatest index in this set. From our previous
reasoning, at a time of collision $t=\tau$ we must necessarily have the chain
of inequalities

(6.8) $v_{i_{*}(\tau)}(\tau-)\geq v_{i_{*}(\tau)+1}(\tau-)\geq\cdots\geq
v_{i^{*}(\tau)-1}(\tau-)\geq v_{i^{*}(\tau)}(\tau-).$

These are exactly the relations that give rise to the barycentric lemma used
in [8, 36] to connect the momentum conservation of particles to the Oleĭnik
E-condition for the scalar balance law.

For any $\alpha\in\mathbb{R}$, the $i$th particle initially moves along its
free trajectory $y_{i}(t)=x_{i}^{0}+tv_{i}^{0}+\frac{1}{2}t^{2}a_{i}$ as long
as one remains in the interior of the cone $\mathbb{K}^{n}$. Now assume that
some subset of particles $\mathcal{J}_{i}(\tau)\ni i$ collide at time $\tau$.
Recalling the set (3.7) and the definition of $X_{n}$ in (6.2), we see that
$(\theta_{i_{*}(\tau)-1},\theta_{i^{*}(\tau)})$ is a maximal interval of
$\Omega_{X_{n}(\tau)}$. For this set of particles, the Rankine–Hugoniot
condition reads

(6.9)
$v_{i}(t)=\frac{\mathcal{U}_{n}(t,\theta_{i^{*}(\tau)})-\mathcal{U}_{n}(t,\theta_{i_{*}(\tau)-1})}{\theta_{i^{*}(\tau)}-\theta_{i_{*}(\tau)-1}}=\frac{\sum_{j\in\mathcal{J}_{i}(\tau)}m_{j}u_{j}(t)}{\sum_{j\in\mathcal{J}_{i}(\tau)}m_{j}}.$

On the other hand, suppose $k\in(\theta_{l-1},\theta_{l}]$ for some
$l\in\mathcal{J}_{i}(\tau)$, then the expressions appearing in the Oleĭnik
E-condition (4.17) are given by

$\frac{\mathcal{U}_{n}(t,k)-\mathcal{U}_{n}(t,\theta_{i_{*}(\tau)-1})}{k-\theta_{i_{*}(\tau)-1}}=\frac{\sum_{j=i_{*}(\tau)}^{l-1}m_{j}u_{j}(t)+(k-\theta_{l-1})u_{l}(t)}{\sum_{j=i_{*}(\tau)}^{l-1}m_{j}+k-\theta_{l-1}}$

and

$\frac{\mathcal{U}_{n}(t,\theta_{i^{*}(\tau)})-\mathcal{U}_{n}(t,k)}{\theta_{i^{*}(\tau)}-k}=\frac{\sum_{j=l+1}^{i^{*}(\tau)}m_{j}u_{j}(t)+(\theta_{l}-k)u_{l}(t)}{\sum_{j=l+1}^{i^{*}(\tau)}m_{j}+\theta_{l}-k}.$

We observe that these expressions attain their extremal values at either
$k=\theta_{l-1}$ or $k=\theta_{l}$ depending on the value of $u_{l}$, and so
it is clear that it is sufficient for us to verify (4.17) for $k=\theta_{l}$
for all $l\in\mathcal{J}_{i}(\tau)$. Moreover, we recall Lemma 5.2, which says
that $U_{n}-\operatorname{P}_{\mathscr{H}_{X_{n}}}U_{n}\in
N_{X_{n}}\mathscr{K}_{\mathbf{m}}$ is equivalent to the Oleĭnik E-condition.
From this it follows that, for these particle dynamics, the Lagrangian and
entropy solutions coincide.

Since the particles in $\mathcal{J}_{i}(\tau)$ collided in the first place, we
know that their free velocities $u_{j}(\tau)=v_{j}(\tau-)$ for
$j\in\mathcal{J}_{i}(\tau)$ satisfy the ordering (6.8), meaning
$\theta\mapsto\mathcal{U}_{n}(t,\theta)$ is concave on
$(\theta_{i_{*}(\tau)-1},\theta_{i^{*}(\tau)})$ at the collision time
$t=\tau$. This implies that the Oleĭnik E-condition is satisfied at $t=\tau$,
but what happens for $t>\tau$ is highly dependent on the sign of $\alpha$,
i.e., if we are in the repulsive or non-repulsive case. Note that if one or
more of the inequalities in (6.8) are strict, there is a loss of kinetic
energy at impact due to Jensen’s inequality and the convexity of $v\mapsto
v^{2}$.

####  6.3.1. The non-repulsive case

The important point in this case is that $\mathcal{A}_{n}(\cdot)$ in (6.5) is
strictly concave on $[0,1]$ for $\alpha>0$, or vanishes for $\alpha=0$. We
already established that $\theta\mapsto\mathcal{U}_{n}(\tau,\theta)$ is
concave on $(\theta_{i_{*}(\tau)-1},\theta_{i^{*}(\tau)})$, and so, since
$\theta\mapsto(t-\tau)\mathcal{A}_{n}(\theta)$ is concave for $t>\tau$, it
follows that
$\mathcal{U}_{n}(t,\theta)=\mathcal{U}_{n}(\tau,\theta)+(t-\tau)\mathcal{A}_{n}(\theta)$
remains concave on this interval. Therefore the Oleĭnik E-condition remains
satisfied, and the entropic solution is to let the particles continue as a
single amassed particle with velocity given by the Rankine–Hugoniot condition
(6.9).

Equivalently, since $\alpha\geq 0$ we have that the accelerations $a_{i}$ from
(6.3) are decreasing in $i$, which combined with (6.8) for
$u_{j}(\tau)=v_{j}(\tau-)$, $j\in\mathcal{J}_{i}(\tau)$ means that this chain
of inequalities remains valid for $u_{j}(t)$, $t>\tau$. If $\alpha>0$, the
inequalities even become strict, so that $\mathbf{u}(t)\notin
T_{\mathbf{x}(t)}\mathbb{K}^{n}$ for all $t>\tau$. Therefore, from here on the
particle velocities will be given by the projection
$\mathbf{v}(t)=\operatorname{P}_{T_{\mathbf{x}(t)}\mathbb{K}^{n}}\mathbf{u}(t)$.
The same is true when $\alpha=0$, as there is no acceleration and the
particles will continue moving with the mass-averaged velocity given by
momentum conservation, or the Rankine–Hugoniot condition, at impact.

Therefore, in the non-repulsive case, once particles collide, they stick
together for all time, exhibiting globally sticky behavior. This is in
agreement with $F_{\mathbf{m}}[X]=A_{n}$ for $\alpha\geq 0$ being a sticking
functional in the sense of Definition 3.5. In particular, for an attractive
potential $\alpha>0$ there is some time $T$ such that for $t>T$ all particles
have been lumped together, moving along the trajectory of the center of mass,
which is unaffected by the internal forces and collision dynamics. Hence we
recover exactly the sticky dynamics already detailed in [36] and [7]. As we
will see next, the behavior for a repulsive potential, i.e., $\alpha<0$, is
less straightforward.

####  6.3.2. The repulsive case

In this case $A_{n}(\theta)$ is strictly convex since $\alpha<0$, which will
change the dynamics, as the flux function will not remain concave on the jump
set for all time after a collision.

For instance, suppose that at the time of impact $\tau$, (6.8) holds with all
inequalities replaced by equality, that is, the particles meet tangentially.
Then $\mathcal{U}(\tau,\theta)$ is linear for
$\theta\in(\theta_{i_{*}(\tau)-1},\theta_{i^{*}(\tau)})$, and for $t>\tau$ the
flux function $\mathcal{U}_{n}(t,\theta)$ will immediately turn strictly
convex on $(\theta_{i_{*}(\tau)-1},\theta_{i^{*}(\tau)})$, violating the
Oleĭnik E-condition there, except on the “trivial” intervals
$(\theta_{i-1},\theta_{i})$, $i\in\mathcal{J}_{i}(\tau)$ corresponding to
single particles. The entropy solution is therefore to have the particles
immediately split apart and follow the free trajectories defined by (6.7), and
one can think of the particles as only grazing one another at $t=\tau$.

On the other hand, if $\mathcal{U}_{n}(\tau,\cdot)$ is strictly concave on
$(\theta_{i_{*}(\tau)-1},\theta_{i^{*}(\tau)})$, it will take some time before
the Oleĭnik E-condition is violated. This coincides with the lower convex
envelope $\theta\mapsto\mathcal{U}_{n,\smallsmile}(t,\theta)$ on this interval
changing from a chord connecting the endpoints of the graph of
$\mathcal{U}_{n}(t,\theta)$ to breaking into two or more affine segments
connected at breakpoints $\theta_{j}$ for $j\in\mathcal{J}_{i}(\tau)$. This
happens when there is either a left subgroup
$\\{i_{*}(\tau),\dots,j\\}\subsetneq\mathcal{J}_{i}(\tau)$ where the average
of their prescribed velocities $u_{j}(t)$ dip below the (average) velocity
$v_{i}(t)$ of the amassed particles $\mathcal{J}_{i}(\tau)$, or alternatively,
there is a right subgroup
$\\{j,\dots,i^{*}(\tau)\\}\subsetneq\mathcal{J}_{i}(\tau)$ such that their
averaged prescribed velocities surpass the velocity of the amassed particle.
The entropy solution is then to have the amassed particle break into two or
more lumped particles corresponding to the aforementioned breakpoints. Since
this “fission” happens when the subgroups of particles have the same velocity
as that of the amassed particle, the kinetic energy does not change as it
breaks apart. This solution can then be thought of as “physical” in the sense
that the energy is nonincreasing.

Since we now know that we can have amassed particles in the repulsive case as
well, if only for a limited time, we observe that similar results hold for
collisions between aggregated particles as well. Two amassed particles can
collide and stick together for some limited time before breaking apart into
new subgroups. All of this is governed by the flux function
$\mathcal{U}_{n}(t,\theta)$. For this repulsive case $\alpha<0$, we see that
there must be some time $T$ such that the flux function is strictly convex for
$t>T$ and there are no more collisions between particles, which then move
freely with velocities $v_{i}(t)=u_{i}(t)$.

In the introduction of [7], the Lagrangian selection principle is motivated by
an inequality involving the collision of two particles. Here we will give a
concrete such example, to illustrate the ideas.

######  Example 6.1 (Two repulsive particles).

Consider the case $\alpha=-2$, $\beta=0$ and two particles where
$m_{1}=m_{2}=\frac{1}{2}$, $x_{1}^{0}=0$, $v_{1}^{0}=2$, $x_{2}^{0}=1$ and
$v_{2}^{0}=0$ such that $a_{1}=-\frac{1}{2}$ and $a_{2}=\frac{1}{2}$. Clearly,
the free trajectories $y_{i}$ and prescribed velocities $u_{i}$ for $i=1,2$
are given by

$y_{1}(t)=2t-\frac{t^{2}}{4},\quad u_{1}(t)=2-\frac{t}{2},\qquad
y_{2}(t)=1+\frac{t^{2}}{4},\quad u_{2}(t)=\frac{t}{2}.$

The undisturbed trajectories intersect at two points, given by the times
$t=2\pm\sqrt{2}$. One way of obtaining a generalized Lagrangian solution is to
project the undisturbed trajectories onto the cone $\mathbb{K}^{2}$, and
writing $\tilde{\mathbf{x}}=\operatorname{P}_{\mathbb{K}^{2}}\mathbf{y}$ we
then obtain

$\tilde{x}_{1}(t)=\begin{cases}2t-\frac{t^{2}}{4},&\lvert
t-2\rvert\geq\sqrt{2},\\\ \frac{1}{2}+t,&\lvert
t-2\rvert<\sqrt{2},\end{cases}\qquad\tilde{x}_{2}(t)=\begin{cases}1+\frac{t^{2}}{4},&\lvert
t-2\rvert\geq\sqrt{2},\\\ \frac{1}{2}+t,&\lvert
t-2\rvert<\sqrt{2}.\end{cases}$

That is, this generalized Lagrangian solution is given by letting the
particles stick together as long as $\mathbf{y}\notin\mathbb{K}^{2}$, i.e.,
$y_{1}(t)>y_{2}(t)$, and otherwise they follow the free trajectories. Note
that these trajectories correspond to the weak solution of the repulsive
Euler–Poisson system found in the similar example [46, Example 3]. On the
other hand, the Lagrangian, or entropic, trajectories are given by letting the
particles stick together as long as the prescribed velocity $\mathbf{u}\notin
T_{\mathbf{x}(t)}\mathbb{K}^{2}$, i.e., $u_{1}(t)>u_{2}(t)$. This means that
the “right” time to split is at $t=2$, and the corresponding trajectories are
given by

$x_{1}(t)=\begin{cases}2t-\frac{t^{2}}{4},&t<2-\sqrt{2},\\\
\frac{1}{2}+t,&2-\sqrt{2}\leq t<2,\\\ \frac{1}{2}+t-\frac{(t-2)^{2}}{4},&t\geq
2,\end{cases}\qquad x_{2}(t)=\begin{cases}1+\frac{t^{2}}{4},&t<2-\sqrt{2},\\\
\frac{1}{2}+t,&2-\sqrt{2}\leq t<2,\\\ \frac{1}{2}+t+\frac{(t-2)^{2}}{4},&t\geq
2.\end{cases}$

(a)

(b)

Figure 1. Free trajectories $y_{i}$ compared with (a) projected trajectories
$\tilde{x}_{i}=(\operatorname{P}_{\mathbb{K}^{2}}\mathbf{y})_{i}$ and (b)
Lagrangian, or entropy, trajectories $x_{i}(t)$ for $i\in\\{1,2\\}$.

Figure 1 highlights the differences between the two solutions. In both cases,
when the particles collide at time $t=2-\sqrt{2}$, the velocity is
discontinuous, as one expects from the conservation of momentum during a
collision, and Jensen’s inequality shows that the kinetic energy decreases on
impact. However, there is an important difference in the splitting of the
particles. For $x_{i}$ the velocity is continuous at $t=2$ when the particles
split, which ensures that the kinetic energy
$\frac{1}{2}(\dot{x}_{1}^{2}+\dot{x}_{2}^{2})$ is conserved at the time of
splitting. On the other hand, for $\tilde{x}_{i}$, the velocity is
discontinuous also at the splitting time $t=2+\sqrt{2}$, which by Jensen’s
inequality shows that the kinetic energy increases at this time. This would
then suggest that the “physically relevant” way of splitting the particles is
given by the entropic solution. Of course, another solution for which the
energy is nonincreasing is to let the particles remain stuck after collision,
but this would not take into account the repulsive forces of the problem, and
violates the Oleĭnik E-condition. Indeed, for this amassed particle, the
inequalities (4.17) are equivalent to $u_{2}(t)\leq\dot{x}_{i}\equiv 1\leq
u_{1}(t)$, which is violated for $t>2$.

######  Remark 6.2 (Numerical algorithm).

Since one can use empirical measures for approximations, cf. [4, Section 4],
one could approximate general initial data for the Euler–Poisson system (1.3)
and evolve it according to the particle dynamics described in this section,
i.e., with particles following parabolic trajectories to approximate
solutions. In practice, one could then implement a front-tracking algorithm in
the spirit of [34], where one considers time-varying velocity fields leading
to nonaffine front trajectories. One could also imagine doing the same with
the more general system (1.1), using an ODE solver to compute the trajectories
between collisions.

##  7\. Repulsive Poisson Forces

In this section we want to use the equivalence we have established to study
some more involved examples. The first of these expands upon an example from
[7], and can be seen as a more involved, continuous version of Example 6.1,
although we here consider symmetric data which simplifies the analysis and
formulas. In the second example and its subcases we use the equivalence
framework on the repulsive Euler–Poisson system with a quadratic confinement
potential studied in [14] with the aim of providing a well-defined notion of
solution beyond singularity formation.

###  7.1. Entropic and non-entropic solutions to the repulsive Euler–Poisson
system

Let us consider the example found in [7, Section 6.2.2] where the authors
consider the Euler–Poisson system (1.3) with $\alpha=-2$ and $\beta=0$, and
initial data $X^{0}(m)=m-\frac{1}{2}$,
$V^{0}(m)=-\operatorname*{sgn}\left(m-\frac{1}{2}\right)$ for the differential
inclusion (3.21). Note that this corresponds to the initial data
$\rho_{0}=\chi_{[-\frac{1}{2},\frac{1}{2}]}$ and
$v_{0}(x)=-\operatorname*{sgn}(x)$ for (1.3), which immediately leads to mass
concentration, or collision of trajectories, at the origin. In the example,
they use the representation formula (3.22) for sticky Lagrangian solutions to
compute a trajectory which fails to satisfy the conditions for a Lagrangian
solution.

We will use this example as a starting point to highlight the equivalence of
Lagrangian solutions, as defined in Section 3, and entropy solutions, as
defined in Section 4, for the repulsive Euler–Poisson system. In particular,
we will compute the Lagrangian solution $X(t,m)$ and its generalized inverse
$M(t,x)$ and show that $M$ is an entropy solution. Furthermore, we will show
that the inverse $\tilde{M}$ of the trajectory $\tilde{X}$ computed from
(3.22), although it is a weak solution, does not satisfy the entropy
condition.

####  7.1.1. Computing the (non-)Lagrangian trajectories

According to (3.29), the the prescribed velocity is

$U(t,m)=-\operatorname*{sgn}\left(m-\frac{1}{2}\right)+t\left(2m-1\right),$

from which we find the free trajectory

$Y(t,m)=(1+t^{2})\left(m-\frac{1}{2}\right)-t\operatorname*{sgn}\left(m-\frac{1}{2}\right).$

To obtain an admissible, i.e., nondecreasing in $m$, trajectory we can simply
project $Y$ onto the cone $\mathscr{K}$, as indicated by the projection
formula (3.22); to this end we introduce the primitive

$\mathcal{Y}(t,m)=\frac{1}{2}(1+t^{2})\left[\left(m-\frac{1}{2}\right)^{2}-\frac{1}{4}\right]+t\left[\frac{1}{2}-\left\lvert
m-\frac{1}{2}\right\rvert\right],$

and then according to Lemma 3.11 we have that the projection of $Y$ onto
$\mathscr{K}$ is
$\tilde{X}=\operatorname{P}_{\mathscr{K}}Y=\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}m}\mathcal{Y}^{**}$.
To find the lower convex envelope $\mathcal{Y}^{**}$ of $\mathcal{Y}$, we note
that $\mathcal{Y}(0,m)$ is convex, while for fixed $t>0$, $\mathcal{Y}(t,m)$
has critical points at $m-\frac{1}{2}=\pm\frac{t}{1+t^{2}}$ and a singular
point at $m=\frac{1}{2}$. The critical points are local minima, while the
singular point is a local maximum, and from this we deduce that the lower
convex envelope is

$\mathcal{Y}^{**}(t,m)=\begin{cases}\mathcal{Y}(t,m),&\lvert
m-\frac{1}{2}\rvert\geq\frac{t}{1+t^{2}},\\\
\mathcal{Y}\left(t,\frac{t}{1+t^{2}}\right),&\lvert
m-\frac{1}{2}\rvert<\frac{t}{1+t^{2}},\end{cases}$

from which we find the corresponding trajectory

$\displaystyle\tilde{X}(t,m)$
$\displaystyle=\begin{cases}(1+t^{2})(m-\frac{1}{2})-t\operatorname*{sgn}(m-\frac{1}{2}),&\lvert
m-\frac{1}{2}\rvert\geq\frac{t}{1+t^{2}},\\\ 0,&\lvert
m-\frac{1}{2}\rvert<\frac{t}{1+t^{2}}.\end{cases}$
$\displaystyle=\operatorname*{sgn}\left(m-\frac{1}{2}\right)\max\left\\{\left(1+t^{2}\right)\left\lvert
m-\frac{1}{2}\right\rvert-t,0\right\\}$

For the Lagrangian, or entropy, solution we will have to instead project the
prescribed velocity $U$ onto the tangent cone $T_{X}\mathscr{K}$,
characterized by Lemma 3.2. Note for a given $t$, $U(t,m)$ is decreasing for
$\lvert m-\frac{1}{2}\rvert<\frac{1}{2t}$ and nondecreasing otherwise. In
particular, we see that it is decreasing on $[0,1]$ for $t<1$. On the other
hand, from the trajectories $Y$ we see that initially we have
$\Omega_{X(t)}=\left\\{m\colon\lvert m\rvert<\frac{t}{1+t^{2}}\right\\}$, and
so for $t<1$ we have

$\operatorname{P}_{T_{X(t)}\mathscr{K}}U(t,m)=\begin{cases}U(t,m),&\lvert
m-\frac{1}{2}\rvert\geq\frac{t}{1+t^{2}}\\\ 0,&\lvert
m-\frac{1}{2}\rvert<\frac{t}{1+t^{2}}.\end{cases}$

Combining the above, we see that at $t=1$, all mass is concentrated at the
origin with zero velocity. For $0\leq t<1$, $\Omega_{X(t)}$ was contained in
the set where $U$ is decreasing. For $t>1$ we have the opposite situation, and
therefore

$\operatorname{P}_{T_{X(t)}\mathscr{K}}U(t,m)=\begin{cases}U(t,m),&\lvert
m-\frac{1}{2}\rvert\geq\frac{1}{2t},\\\ 0,&\lvert
m-\frac{1}{2}\rvert<\frac{1}{2t}\end{cases}$

in this case. Integrating in time we obtain the Lagrangian solution

(7.1)
$X(t,m)=\begin{cases}(1+t^{2})(m-\frac{1}{2})-t\operatorname*{sgn}(m-\frac{1}{2}),&\frac{t}{1+t^{2}}\leq\left\lvert
m-\frac{1}{2}\right\rvert\leq\min\left\\{\frac{1}{2t},\frac{1}{2}\right\\},\\\
0,&\left\lvert
m-\frac{1}{2}\right\rvert<\min\left\\{\frac{t}{1+t^{2}},\frac{1}{2t}\right\\},\\\
\left(m-\frac{1}{2}\right)\left(t-\frac{1}{\lvert
2m-1\rvert}\right)^{2},&\left\lvert
m-\frac{1}{2}\right\rvert\geq\min\left\\{\frac{1}{2t},\frac{1}{2}\right\\}.\end{cases}$

The generalized Lagrangian solution $\tilde{X}$ and the Lagrangian solution
$X$ are compared in Figure 2, where we observe that the trajectories emanating
from mass diffusion in $X$ are tangential, in contrast to $\tilde{X}$. Note
that in the original example of [7], they use the characterization of Lemma
3.1 to show that $U-\dot{\tilde{X}}\notin N_{\tilde{X}}\mathscr{K}$, hence
$\tilde{X}$ is not a Lagrangian solution. In particular, they provide a
counterexample to the equivalent condition

(7.2) $\frac{\partial}{\partial
t}\left(\mathcal{Y}(t,m)-\mathcal{Y}^{**}(t,m)\right)\geq 0,\qquad
m\in\Omega_{\tilde{X}(t)}.$

This is then again equivalent to the Oleĭnik E-condition (4.17).

(a)

(b)

Figure 2. The free trajectories $Y(t,\theta_{i})$ for
$\theta_{i}=\frac{i}{30}$, $i\in\\{0,1,\dots,30\\}$, together with the (a)
generalized Lagrangian solution $\tilde{X}(t,\theta_{i})$ and (b) Lagrangian
solution $X(t,\theta_{i})$. Note that $\tilde{X}(t,m)=X(t,m)$ for $t\in[0,1]$.

(a)

(b)

Figure 3. The (a) non-entropic $\tilde{M}(t,\cdot)$ and (b) entropic
$M(t,\cdot)$ weak solution for various times $t$. These are the respective
generalized inverses of $\tilde{X}(t,\cdot)$ and $X(t,\cdot)$ from Figure 2.
Note that $\tilde{M}(t,\cdot)=M(t,\cdot)$ for $t\in[0,1]$.

####  7.1.2. Entropy admissibility

Next we will consider this problem from the point of view of the scalar
conservation law (4.4), which in this case has the flux function

$\mathcal{U}(t,m)=t\left[\left(m-\frac{1}{2}\right)^{2}-\frac{1}{4}\right]+\frac{1}{2}-\left\lvert
m-\frac{1}{2}\right\rvert.$

Some computations show that the corresponding right-continuous inverses of
$\tilde{X}$ and $X$ are respectively

$\tilde{M}(t,x)=\begin{cases}\frac{1}{2}\left(1+\operatorname*{sgn}(x)\right),&\lvert
x\rvert>\frac{1}{2}(t-1)^{2},\\\
\frac{1}{2}+\frac{x+t\operatorname*{sgn}^{\\!+}(x)}{1+t^{2}},&\lvert
x\rvert\leq\frac{1}{2}(t-1)^{2}\end{cases}$

and

$M(t,x)=\begin{cases}\tilde{M}(t,x),&t<1,\\\
\begin{cases}\frac{1}{2}\left(1+\operatorname*{sgn}(x)\right),&\lvert
x\rvert>\frac{1}{2}(t-1)^{2},\\\
\frac{1}{2}+\operatorname*{sgn}^{\\!+}(x)\left(\frac{\sqrt{\lvert
x\rvert}+\sqrt{2t+\lvert x\rvert}}{2t}\right)^{2},&\lvert
x\rvert\leq\frac{1}{2}(t-1)^{2}\end{cases},&t\geq 1,\end{cases}$

where we have introduced the right-continuous sign function
$\operatorname*{sgn}^{\\!+}(x)$, i.e., $\operatorname*{sgn}(x)$ with the
altered convention $\operatorname*{sgn}^{\\!+}(0)=1$. Further computations,
rather cumbersome in the case of $M$, show that the generalized inverses are
piecewise smooth solutions of (4.4), with a jump discontinuity, or stationary
shock, at the origin for $t>0$. Hence, to be weak and entropy solutions it is
sufficient that they satisfy respectively the Rankine–Hugoniot condition and
Oleĭnik E-condition at the shock. The symmetry of the flux function around
$m=\frac{1}{2}$ yields
$[[\mathcal{U}(t,\tilde{M}(t,0))]]=0=[[\mathcal{U}(t,M(t,0))]]$, so we see
that both solutions satisfy the Rankine–Hugoniot condition at the stationary
shock.

In order to check the entropy conditions we have to study the convexity
properties of the flux function $\mathcal{U}(t,m)$. To this end we note that
$\mathcal{U}(0,m)$ is concave, and for a given $t>0$, $\mathcal{U}(t,m)$ has
critical points at $t\lvert 2m-1\rvert=1$ and a singular point at
$m=\frac{1}{2}$. Again, the critical points are local minima and the singular
point is a local maximum, showing that $\mathcal{U}(t,m)$ is convex for
$\lvert m-\frac{1}{2}\rvert\geq\frac{1}{2t}$ and concave for $\lvert
m-\frac{1}{2}\rvert\leq\frac{1}{2t}$.

Now we consider the sets in $[0,1]$ where $\tilde{X}(t,\cdot)$ and
$X(t,\cdot)$ are constant, namely

$\Omega_{\tilde{X}(t)}=\left\\{m\colon\left\lvert
m-\frac{1}{2}\right\rvert<\frac{t}{1+t^{2}}\right\\},\qquad\Omega_{X(t)}=\left\\{m\colon\left\lvert
m-\frac{1}{2}\right\rvert<\min\left\\{\frac{t}{1+t^{2}},\frac{1}{2t}\right\\}\right\\}.$

Observe that $\mathcal{U}(t,m)$ is always concave on the set $\Omega_{X(t)}$,
and since both this set and the flux function $\mathcal{U}(t,m)$ is symmetric
around $m=\frac{1}{2}$, its lower convex envelope
$\mathcal{U}_{\smallsmile}(t,m)$ on this set will be a horizontal chord
connecting the endpoints. In consequence, we have
$\mathcal{U}(t,m)>\mathcal{U}_{\smallsmile}(t,m)$ for $m\in\Omega_{X(t)}$,
which in turn implies the Oleĭnik E-condition (4.17); compare this to the
condition (7.2). On the other hand, on the set
$\Omega_{\tilde{X}(t)}\setminus\Omega_{X(t)}$, which becomes nonempty only for
$t>1$, $\mathcal{U}(t,m)$ is convex. Therefore, the lower convex envelope on
the set $\Omega_{\tilde{X}(t)}$ will have a negative slope for
$-\frac{t}{1+t^{2}}<m<-\frac{1}{2t}$ and a positive slope for
$\frac{1}{2t}<m<\frac{t}{1+t^{2}}$, and this shows that the stationary shock
of $\tilde{M}$ violates (4.17). These two solutions are compared in Figure 3,
where we observe, as seen from the formulas, that the shock decreases faster
for the entropy solution $M$.

####  7.1.3. Energy of trajectories

Recall that the energy for the system is the sum of the kinetic and potential
energy, which for our Euler–Poisson system with $\beta=0$ reads

$\frac{1}{2}\int_{\mathbb{R}}v(t,x)^{2}\mathop{}\\!\mathrm{d}\rho(t,x)+\frac{\alpha}{4}\iint_{\mathbb{R}^{2}}\lvert
x-y\rvert\mathop{}\\!\mathrm{d}\rho(t,x)\mathop{}\\!\mathrm{d}\rho(t,y).$

Using the push-forward relation and $\alpha=-2$ this becomes

$\frac{1}{2}\int_{0}^{1}V(t,m)^{2}\mathop{}\\!\mathrm{d}m-\frac{1}{2}\iint_{[0,1]^{2}}\left\lvert
X(t,m)-X(t,\omega)\right\rvert\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}\omega.$

For the initial configuration we can easily compute

$\int_{0}^{1}(V^{0})^{2}(m)\mathop{}\\!\mathrm{d}m=1,\qquad\iint_{[0,1]^{2}}\left\lvert
X^{0}(m)-X^{0}(\omega)\right\rvert\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}\omega=\frac{1}{3},$

meaning the initial energy is $\frac{1}{3}$. At $t=1$, all mass is
concentrated at the origin with zero velocity, and so both the kinetic and
potential energy are zero, i.e., the initial amount of energy has dissipated.
Now, given the formula (7.1) for the entropic solution for $t\geq 1$, it is
not too difficult to compute the kinetic energy, but considerably more
involved to compute the potential energy, and find the identity

$\int_{0}^{1}V(t,m)^{2}\mathop{}\\!\mathrm{d}m=\iint_{[0,1]^{2}}\left\lvert
X(t,m)-X(t,\omega)\right\rvert\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}\omega=\frac{(t-1)^{3}}{3t},\qquad
t\geq 1.$

Hence, as the system tends to decrease its potential energy, for $t\geq 1$
there is no more dissipation of energy for the entropic solution $(X,V)$, as
the kinetic and potential energy exactly balance one another. On the other
hand, since the evolution of $\tilde{X}$ is “mirrored” around $t=1$, it
follows that the energy will increase again and at $t=2$ it will regain its
initial energy $\frac{1}{3}$. Therefore, as in the discrete case, we see that
the Lagrangian, or entropy, solution provides a more physically desirable
behavior in the sense that energy is nonincreasing.

######  Remark 7.1 (The globally sticky solution).

In the above we could also have considered a globally sticky solution where
for $0\leq t<1$ the evolution is identical to $X$ and $\tilde{X}$, while for
$t\geq 1$ all mass is concentrated in the stationary shock at the origin. This
would also be a generalized Lagrangian solution, as the zero element is always
contained in the tangent cone $T_{X}\mathscr{K}$, as well as a weak solution
of the scalar conservation law since the Rankine–Hugoniot condition is
satisfied for the generalized inverse. In addition, there would clearly be no
more dissipation of energy. However, from analogous arguments as before, it
fails to satisfy the entropy condition, as $\mathcal{U}(t,m)$ is no longer
concave on the whole of $[0,1]$ for $t>1$.

####  7.1.4. General energy considerations

Consider the repulsive Euler–Poisson system (1.3) with no external forcing,
i.e. $\alpha<0$ and $\beta=0$. Once more we write the energy in Lagrangian
variables to obtain

$E(t)\coloneqq\frac{1}{2}\int_{0}^{1}V(t,m)^{2}\mathop{}\\!\mathrm{d}m+\frac{\alpha}{4}\iint_{[0,1]^{2}}\left\lvert
X(t,m)-X(t,\omega)\right\rvert\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}\omega.$

Since $V=\operatorname{P}_{T_{X(t)}\mathscr{K}}U$ we have $\lVert
V\rVert_{L^{2}(0,1)}\leq\lVert U\rVert_{L^{2}(0,1)}$, we can instead use an
expression involving $U$ to bound the energy, and the advantage of this is
that we have a differential equation for $U$. That is, for

$\hat{E}(t)\coloneqq\frac{1}{2}\int_{0}^{1}U(t,m)^{2}\mathop{}\\!\mathrm{d}m+\frac{\alpha}{4}\iint_{[0,1]^{2}}\left\lvert
X(t,m)-X(t,\omega)\right\rvert\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}\omega$

we find

$\displaystyle\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\hat{E}(t)$
$\displaystyle=\frac{-\alpha}{2}\int_{0}^{1}U(t,m)(2m-1)\mathop{}\\!\mathrm{d}m+\frac{\alpha}{4}\iint_{[0,1]^{2}}\operatorname*{sgn}(X(t,m)-X(t,\omega))(V(t,m)-V(t,\omega))\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}\omega$
$\displaystyle=\frac{-\alpha}{2}\int_{0}^{1}U(t,m)(2m-1)\mathop{}\\!\mathrm{d}m+\frac{\alpha}{2}\int_{0}^{1}V(t,m)\int_{0}^{1}\operatorname*{sgn}(X(t,m)-X(t,\omega))\mathop{}\\!\mathrm{d}\omega\mathop{}\\!\mathrm{d}m.$

Now we note that

$\int_{0}^{1}\operatorname*{sgn}(X(t,m)-X(t,\omega))\mathop{}\\!\mathrm{d}\omega=\begin{cases}\displaystyle\int_{0}^{m}\mathop{}\\!\mathrm{d}\omega-\int_{m}^{1}\mathop{}\\!\mathrm{d}\omega=2m-1,&\text{for
a.e.\ }m\notin\Omega_{X(t)}\\\\[5.69054pt]
\displaystyle\int_{0}^{m_{l}}\mathop{}\\!\mathrm{d}\omega-\int_{m_{r}}^{1}\mathop{}\\!\mathrm{d}\omega=m_{l}+m_{r}-1,&m\in(m_{l},m_{r})\subset\Omega_{X(t)},\end{cases}$

that is, for a.e. $m\in\Omega$ we have

(7.3)
$\int_{0}^{1}\operatorname*{sgn}(X(t,m)-X(t,\omega))\mathop{}\\!\mathrm{d}\omega=\operatorname{P}_{\mathscr{H}_{X(t)}}(2m-1).$

Rewriting a bit we then have

$\displaystyle\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\hat{E}(t)=-\frac{\alpha}{2}\int_{0}^{1}(U(t,m)-V(t,m))(2m-1)\mathop{}\\!\mathrm{d}m-\frac{\alpha}{2}\int_{0}^{1}V(t,m)\left(2m-1-\operatorname{P}_{\mathscr{H}_{X(t)}}(2m-1)\right)\mathop{}\\!\mathrm{d}m.$

In the above, the final term vanishes due to $V\in\mathscr{H}_{X(t)}$ and
$2m-1-\operatorname{P}_{\mathscr{H}_{X(t)}}(2m-1)\in\mathscr{H}_{X(t)}^{\perp}$.
On the other hand, integrating by parts, we have

$\int_{0}^{1}(U(t,m)-V(t,m))(2m-1)\mathop{}\\!\mathrm{d}m=-2\int_{0}^{1}\int_{0}^{m}(U(t,\omega)-V(t,\omega))\mathop{}\\!\mathrm{d}\omega\mathop{}\\!\mathrm{d}m\leq
0,$

where the final inequality comes from $U-V\in N_{X(t)}\mathscr{K}$ and the
characterization provided by Lemma 3.1. Since $\alpha<0$ we then see that
$\hat{E}(t)$ is decreasing in time. Hence,
$E(t)\leq\hat{E}(t)\leq\hat{E}(0)=E(0)$. Note that a different strategy is
needed for the attractive case, as the acceleration will make the prescribed
velocity $\lVert U\rVert_{L^{2}(0,1)}$ increase through collisions.

###  7.2. Repulsive Poisson and quadratic confinement forces

As we have seen, internal forces of the Euler–Poisson system (1.3) are either
repulsive or attractive depending on the sign of $\alpha$. The dynamics become
even more complex when both attractive and repulsive forces are present: for
instance, one could combine repulsive Poisson forces with a quadratic
confinement potential, like in [14]. Following these lines, we modify the
damped Euler–Poisson system (4.10) by taking $\alpha=-\lambda^{2}<0$,
$\gamma\geq 0$ and $\Phi(t,x)=\phi\ast\rho(t,x)$ with $\phi(x)=\lvert
x\rvert-\frac{1}{2}x^{2}$; this gives a repulsive Poisson potential with
quadratic confinement and linear damping as described by (1.6), namely

$\partial_{t}\rho+\partial_{x}\left(\rho v\right)=0,\qquad\partial_{t}(\rho
v)+\partial_{x}\left(\rho v^{2}\right)=-\gamma\rho
v+\lambda^{2}\rho\int_{\mathbb{R}}\left(\operatorname*{sgn}(x-y)-(x-y)\right)\mathop{}\\!\mathrm{d}\rho(t,y)-\beta\rho,$

where we observe that the effect of the confinement potential is to drive
$\rho$ toward the center of mass with a force that increases linearly with
distance. As mentioned in Example 4.4, this force distribution $f[\rho,v]$ is
not directly covered by (1.1) since it also involves the velocity $v$.
However, the same argument as in [7, Theorem 3.6] can be adapted for this case
as well. Indeed, one can still consider the operator $A(X,U)=(\partial
I_{\mathscr{K}}(X)-U,F[X,U])$ as a Lipschitz perturbation of $\partial
I_{\mathscr{K}}(X)$, as long as $F[X,U]$ is Lipschitz in both $X$ and $U$.

####  7.2.1. Derivation of the prescribed velocity and flux function in the
smooth case

Note that, using the a priori knowledge that $\rho$ is a probability measure,
we can, as for the Euler–Poisson system, rewrite
$(\operatorname*{sgn}\ast\rho)\rho=\partial_{x}(M^{2}-M)$ as before. On the
other hand, under the assumption that $\rho(t)\in\mathscr{P}_{2}(\mathbb{R})$,
the term $(\mathrm{Id}\ast\rho)\rho$ is a well-defined measure since
$(\mathrm{Id}\ast\rho)$ is continuous for a given $t$. Indeed, for such $\rho$
we should have $\lvert(\mathrm{Id}\ast\rho)(t,x)\rvert\leq 1+\lvert
x\rvert+\overline{x^{2}}(t)$, where $\overline{x^{2}}(t)$ is the second moment
at time $t$. Moreover,
$(\mathrm{Id}\ast\rho)(t,x)=x+(\mathrm{Id}\ast\rho)(t,0)=x-\bar{x}(t)$ and
$(\mathrm{Id}\ast\rho)(t,x)-(\mathrm{Id}\ast\rho)(t,y)=x-y$, where we have
introduced the center of mass, or first moment, $\bar{x}(t)$. That is, for
fixed $t$ it is a locally bounded, linear function in $x$, with Lipschitz-
constant 1. Integrating $(\mathrm{Id}\ast\rho)\rho$ we obtain

$\int_{\mathbb{R}}\mathop{}\\!\mathrm{d}\rho(t,y)\int_{(-\infty,x]}y\mathop{}\\!\mathrm{d}\rho(t,y)-\int_{\mathbb{R}}y\mathop{}\\!\mathrm{d}\rho(t,y)\int_{(-\infty,x]}\mathop{}\\!\mathrm{d}\rho(t,y)=\int_{(-\infty,x]}y\mathop{}\\!\mathrm{d}\rho(t,y)-\bar{x}(t)M(t,x).$

Therefore we see that this term cannot directly be expressed in terms of the
distribution function $M$, contrary to the primitive of
$(\operatorname*{sgn}\ast\rho)\rho$.

However, we can still apply some a priori knowledge to the trajectory of the
center of mass. This should not be affected by the internal forces, i.e., the
forces coming from $\phi^{\prime}\ast\rho$, but only the external damping.
Multiplying the continuity equation with $x$ and integrating we obtain

$\dot{\bar{x}}(t)=\partial_{t}\int_{\mathbb{R}}x\mathop{}\\!\mathrm{d}\rho(t,x)=\int_{\mathbb{R}}v\mathop{}\\!\mathrm{d}\rho(t,x)\eqqcolon\bar{v}(t),$

that is, the velocity of the center of mass is given by the total momentum.
Furthermore, integrating the momentum equation, the internal forces vanish due
to symmetry, or, if one prefers, Newton’s third law of motion, and we are left
with

$\dot{\bar{v}}(t)=-\gamma\bar{v}(t)-\beta.$

Integrating these equations we then obtain

(7.4) $\bar{x}(t)=\bar{x}(0)+\bar{v}(0)\frac{1-\mathrm{e}^{-\gamma
t}}{\gamma}+\beta\frac{1-\gamma t-\mathrm{e}^{-\gamma
t}}{\gamma^{2}}\qquad\bar{v}(t)=\mathrm{e}^{-\gamma
t}\bar{v}(0)-\beta\frac{1-\mathrm{e}^{-\gamma t}}{\gamma},$

where the case $\gamma=0$ corresponds to the limit $\gamma\to 0$. These
relations should remain valid no matter how one resolves the breakdown of
classical solutions, i.e., concentration of mass.

Now we want to derive the flux function for our problem, by following the
arguments leading up to (4.3). We then suggestively introduce

$\displaystyle u(t,X(t,m))$ $\displaystyle=v_{0}\circ
X^{0}(m)+\int_{0}^{t}F[X](m)\mathop{}\\!\mathrm{d}s$
$\displaystyle=V^{0}(m)+\int_{0}^{t}\left[-\gamma
u(s,X(s,m))+\lambda^{2}\left(2m-1-X(s,m)+\bar{x}(s)\right)\right]\mathop{}\\!\mathrm{d}s,$

where $F[X]$ is in accordance with (3.16) for our right-hand side. This
relation together with $\dot{X}(t,m)=u(t,X(t,m))=U(t,m)$ should be satisfied
by the solution as long as it remains smooth; alternatively we can write

$\dot{X}(t,m)=U(t,m),\qquad\dot{U}(t,m)=-\gamma
U(t,m)+\lambda^{2}\left(2m-1+\bar{x}(t)-X(t,m)\right),$

which should be compared to the smooth trajectories of the numerical example
[7, Equation (1.40)]. Differentiating we obtain the second-order inhomogeneous
ODE

$\ddot{U}+\gamma\dot{U}+\lambda^{2}U=\lambda^{2}\bar{v}(t),\quad$

with initial conditions

$U(0,m)=V^{0}(m),\qquad\dot{U}(0,m)=-\gamma
V^{0}(m)+\lambda^{2}(2m-1+\bar{x}(0)-X^{0}(m)).$

Note that since the total momentum, or mean velocity, $\bar{v}$ from (7.4)
satisfies $\ddot{\bar{v}}+\gamma\dot{\bar{v}}=0$, the above can be restated as
a homogeneous ODE for the deviation $\tilde{U}=U-\bar{v}$. Let us from here on
write $\gamma=2\kappa$ to simplify expressions. We are interested in the case
where classical solutions oscillate, and so we take $\lambda>\kappa\geq 0$.
Solving the ODE we obtain the following expression for the smooth prescribed
velocity,

(7.5) $\displaystyle U(t,m)$
$\displaystyle=\bar{v}(t)+\left(V^{0}(m)-\bar{v}(0)\right)\mathrm{e}^{-\kappa
t}\cos\left(\sigma t\right)$
$\displaystyle\quad+\left[\beta-\kappa\left(V^{0}(m)-\bar{v}(0)\right)-\lambda^{2}\left(X^{0}(m)-2m+1-\bar{x}(0)\right)\right]\mathrm{e}^{-\kappa
t}\frac{\sin\left(\sigma t\right)}{\sigma},$

where we have introduced $\sigma\coloneqq\sqrt{\lambda^{2}-\kappa^{2}}$. The
corresponding flux function is then given by the primitive of this velocity,
which is of the form (4.18). Indeed, this follows since
$\rho_{0}\in\mathscr{P}_{2}(\mathbb{R})$ means $\mathrm{Id}\in
L^{2}(\mathbb{R},\rho_{0})$ which again implies $\mathrm{Id}\circ
X^{0}=X^{0}\in L^{2}(0,1)$, while $v_{0}\in L^{2}(\mathbb{R},\rho_{0})$ leads
to $V^{0}=v_{0}\circ X^{0}\in L^{2}(0,1)$. Writing

$\mathcal{X}^{0}(m)=\int_{0}^{m}X^{0}(\omega)\mathop{}\\!\mathrm{d}\omega,\qquad\mathcal{V}^{0}(m)=\int_{0}^{m}V^{0}(\omega)\mathop{}\\!\mathrm{d}\omega=\int_{0}^{m}v_{0}(X^{0}(\omega))\mathop{}\\!\mathrm{d}\omega,$

the flux function can be expressed as

$\displaystyle\mathcal{U}(t,m)$
$\displaystyle=\bar{v}(t)m+\left(\mathcal{V}^{0}(m)-\bar{v}(0)m\right)\mathrm{e}^{-\kappa
t}\cos\left(\sigma t\right)$ $\displaystyle\quad+\left[\beta
m-\kappa\left(\mathcal{V}^{0}(m)-\bar{v}(0)m\right)-\lambda^{2}\left(\mathcal{X}^{0}(m)-m^{2}+m-\bar{x}(0)m\right)\right]\mathrm{e}^{-\kappa
t}\frac{\sin\left(\sigma t\right)}{\sigma}$
$\displaystyle=m\frac{\mathop{}\\!\mathrm{d}\bar{x}}{\mathop{}\\!\mathrm{d}t}+\left(\mathcal{V}^{0}(m)-\bar{v}(0)m\right)\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\mathrm{e}^{-\kappa
t}\frac{\sin(\sigma t)}{\sigma}$ $\displaystyle\quad+\left[\beta
m-\kappa\left(\mathcal{V}^{0}(m)-\bar{v}(0)m\right)-\lambda^{2}\left(\mathcal{X}^{0}(m)-m^{2}+m-\bar{x}(0)m\right)\right]\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\mathrm{e}^{-\kappa
t}\left(\cos(\sigma t)+\frac{\kappa}{\sigma}\sin(\sigma t)\right).$

######  Remark 7.2 (Comparison with existing formulas).

Let us verify that this agrees with the formulas found in [14], where we first
note that their mass is not normalized to one. In their formulas they denote
the total (conserved) mass by $M_{0}$ and the initial momentum by $M_{1}$; we
will instead denote these by respectively $\mu_{0}$ and $\mu_{1}$ to avoid
confusion with the distribution function $M$. Moreover, the parameters for
their system are $\lambda=1$, $\gamma=2\kappa=1$ and $\beta=0$, such that
$1<4\mu_{0}$ is what they call case C, exhibiting oscillatory behavior. We
would therefore have to rescale $\rho$ to obtain the corresponding set of
equations. Writing their density as $\tilde{\rho}=\mu_{0}\rho$ for
$\rho\in\mathscr{P}(\mathbb{R})$, this amounts to the rescaling
$\lambda^{2}\mapsto\lambda^{2}\mu_{0}$ and the identities
$\sigma=\frac{1}{2}\sqrt{4\mu_{0}-1}$, $\bar{v}(0)=\mu_{1}/\mu_{0}$ in our
equations. Now, for this set of parameters we would like our expression (7.5)
to coincide with the smooth characteristic velocity given in [14, Equation
(2.8)]. However, their characteristics $\eta(t,x)$, satisfying
$\partial_{t}\eta(t,x)=u(t,\eta(t,x))$ and $\eta(0,x)=x$, are parametrized by
the initial position $x\in\mathbb{R}$ rather than the material coordinate
$m\in[0,1]$ for our $X(t,m)$. Since we are in the smooth case, we have the
identity $M(t,\eta(t,x))=M^{0}(x)$, and so evaluating (7.5) at $m=M^{0}(x)$
yields the characteristic velocity parametrized by $x$. In addition we have
the identities $X^{0}(M^{0}(x))=x$ and $V^{0}(M^{0}(x))=v_{0}(x)$ which hold
$\rho_{0}$-a.e., and so we plug into (7.5) to find

$\displaystyle U(t,M^{0}(x))$
$\displaystyle=\frac{\mu_{1}}{\mu_{0}}\mathrm{e}^{-t}+\left(v_{0}(x)-\frac{\mu_{1}}{\mu_{0}}\right)\mathrm{e}^{-t/2}\cos\left(\frac{\sqrt{4\mu_{0}-1}}{2}t\right)$
$\displaystyle\quad+\frac{2}{\sqrt{4\mu_{0}-1}}\left[\mu_{0}\left(2M^{0}(x)-1+\bar{x}(0)-x\right)-\frac{1}{2}\left(v_{0}(x)-\frac{\mu_{1}}{\mu_{0}}\right)\right]\mathrm{e}^{-t/2}\sin\left(\frac{\sqrt{4\mu_{0}-1}}{2}t\right),$

which is exactly the formula we are after, since we have the relations

$\mu_{0}M^{0}(x)=\int_{-\infty}^{x}\tilde{\rho}_{0}(y)\mathop{}\\!\mathrm{d}y,\qquad\mu_{0}\bar{x}(t)=\int_{\mathbb{R}}y\tilde{\rho}(t,y)\mathop{}\\!\mathrm{d}y.$

Hence, we are in agreement with [14] as long as the solutions remain
classical, which is as far as their analysis goes.

Next, with the Lagrangian and entropy solution framework in hand, we aim to
provide well-defined solutions beyond the breakdown of classical solutions,
and for this we give some examples.

####  7.2.2. A simple (computable) case

Let us, as in Section 7.1, take $X^{0}(m)=m-\frac{1}{2}$, but, as explained
below, we will for simplicity take $V^{0}(m)=-\left(m-\frac{1}{2}\right)$.
Then it is clear that the initial center of mass is $\bar{x}(0)=0$ and the
initial total momentum, or average velocity, is $\bar{v}(0)=0$, which
simplifies the flux function a bit. To simplify even more, we let the external
force be zero, $\beta=0$. Computing $\mathcal{X}^{0}(m)$ we find

$\mathcal{X}^{0}(m)=\int_{0}^{m}X^{0}(\omega)\mathop{}\\!\mathrm{d}\omega=\frac{m^{2}-m}{2},$

which leads to the flux function

$\displaystyle\mathcal{U}(t,m)$
$\displaystyle=\mathcal{V}^{0}(m)\mathrm{e}^{-\kappa t}\left[\cos\left(\sigma
t\right)-\frac{\kappa}{\sigma}\sin\left(\sigma
t\right)\right]+\frac{\lambda^{2}}{2}\left(m^{2}-m\right)\mathrm{e}^{-\kappa
t}\frac{\sin(\sigma t)}{\sigma}$
$\displaystyle=-\frac{1}{2}\left(m^{2}-m\right)\mathrm{e}^{-\kappa
t}\left[\cos\left(\sigma
t\right)-\frac{\kappa+\lambda^{2}}{\sigma}\sin\left(\sigma t\right)\right].$

Then the prescribed velocity $U=\frac{\partial}{\partial m}\mathcal{U}$ is

$U(t,m)=-\left(m-\frac{1}{2}\right)\mathrm{e}^{-\kappa
t}\left[\cos\left(\sigma
t\right)-\frac{\kappa+\lambda^{2}}{\sigma}\sin\left(\sigma t\right)\right]$

which we integrate to find the free trajectory

$\displaystyle Y(t,m)$ $\displaystyle=X^{0}(m)+V^{0}(m)\mathrm{e}^{-\kappa
t}\frac{\sin\left(\sigma
t\right)}{\sigma}-\frac{m^{2}-m}{2}\mathrm{e}^{-\kappa
t}\left[\cos\left(\sigma t\right)+\frac{\kappa}{\sigma}\sin\left(\sigma
t\right)\right]$
$\displaystyle=\left(m-\frac{1}{2}\right)\left[2-\mathrm{e}^{-\kappa
t}\left(\cos\left(\sigma t\right)+\frac{1+\kappa}{\sigma}\sin\left(\sigma
t\right)\right)\right].$

Since $V^{0}$ is decreasing, the flux function is initially concave, and will
periodically alternate between being concave and convex in between the times
$t$ satisfying $\cot(\sigma t)=\frac{\kappa+\lambda^{2}}{\sigma}$, which
coincide with the times when the free trajectory $Y(\cdot,m)$ has its extremal
amplitudes. Moreover, since mass concentrates when the free trajectories are
no longer increasing, we see that either all or no mass concentrates depending
on the function

$c(t)\coloneqq 2-\mathrm{e}^{-\kappa t}\left(\cos\left(\sigma
t\right)+\frac{1+\kappa}{\sigma}\sin\left(\sigma t\right)\right),$

which attains its minimum $c(\tau_{*})$ at $t=\tau_{*}$, where

$\tau_{*}=\arctan\left(\frac{\sigma}{\kappa+\lambda^{2}}\right),\qquad
c(\tau_{*})=2-\exp\left(-\frac{\kappa}{\sigma}\arctan\left(\frac{\sigma}{\kappa+\lambda^{2}}\right)\right)\frac{\sqrt{\lambda^{2}+2\kappa+1}}{\lambda}.$

If $c(\tau_{*})$ is nonnegative, there is no mass concentration, except
possibly at the instant $t=\tau_{*}$ if $c(\tau_{*})=0$. However, in this
case, the velocity of the free trajectories will all be zero, and thus
tangential to one another, and the Lagrangian solution will still be given by
the smooth curves given by the free trajectory $Y(t,m)$. On the other hand, if
$c(\tau_{*})<0$, then there is some $\tau_{0}\in(0,\tau_{*})$ such that
$c(\tau_{0})=0$ and all mass collides at the origin. We note that in the
undamped case $\kappa=0$ we explicitly find

$\tau_{0}=\arcsin\left(\frac{\lambda\left(2-\sqrt{1-3\lambda^{2}}\right)}{1+\lambda^{2}}\right),$

which also emphasizes that mass concentration can only occur for
$0<\lambda<1/\sqrt{3}$ in the undamped case. From time $t=\tau_{0}$ onward we
can no longer use the smooth free trajectory found by solving

$\dot{X}(t,m)=U(t,m)=V^{0}(m)+\int_{0}^{t}\left[-2\kappa
U(s,m)+\lambda^{2}\left(2m-1-X(s,m)\right)\right]\mathop{}\\!\mathrm{d}s,$

since $X(t,m)=0$ due to the concentration at the origin. Instead we take the
velocity $U(\tau_{0},m)$ as initial condition for the equation

$U(t,m)=U(\tau_{0},m)+\int_{\tau_{0}}^{t}\left[-2\kappa
U(s,m)+\lambda^{2}\left(2m-1\right)\right]\mathop{}\\!\mathrm{d}s$

and find

$\displaystyle U(t)$
$\displaystyle=\mathrm{e}^{-2\kappa(t-\tau_{0})}U(\tau_{0},m)+\lambda^{2}\mathrm{e}^{-2\kappa
t}(t-\tau_{0})(2m-1)$
$\displaystyle=\left(m-\frac{1}{2}\right)\mathrm{e}^{-2\kappa
t}\left[2\lambda^{2}(t-\tau_{0})-\mathrm{e}^{\kappa\tau_{0}}\left(\cos(\sigma\tau_{0})-\frac{\kappa}{\sigma}\sin(\sigma\tau_{0})\right)\right].$

Since $U(\tau_{0},m)$ is decreasing in $m$, meaning $\mathcal{U}(\tau_{0},m)$
is concave, we see that the mass will remain at the origin until $U(t,m)$
becomes increasing at times $t>\tau_{1}$ where

$\tau_{1}=\tau_{0}+\frac{\mathrm{e}^{\kappa\tau_{0}}}{2\lambda^{2}}\left(\cos(\sigma\tau_{0})-\frac{\kappa}{\sigma}\sin(\sigma\tau_{0})\right),$

and we note that $U(\tau_{1},m)=0$. From this point on we solve the original
second order, inhomogeneous equation with initial data
$X(\tau_{1},m)=U(\tau_{1},m)=0$ to find

$X(t,m)=\left(2m-1\right)\left[1-\mathrm{e}^{-\kappa(t-\tau_{1})}\cos\left(\sigma(t-\tau_{1})\right)\right].$

Note that even in the undamped case $\kappa=0$, although some kinetic energy
is lost at the first collision $t=\tau_{0}$, for the remaining collisions
happening at $t=\tau_{1}+n2\pi/\lambda$ for $n\in\mathbb{N}$, the velocities
are tangential and there is no loss of energy. On the other hand, in the
damped case we observe that $X(t,m)\to 2m-1$ as $t\to\infty$, corresponding to
the density $\rho(t,x)$ tending to $\frac{1}{2}\chi_{[-1,1]}(x)$.

In total, the Lagrangian solution is given by

$X(t,m)=\begin{cases}Y(t,m),&c(\tau_{*})\geq 0,\\\ \begin{cases}Y(t,m),&0\leq
t<\tau_{0},\\\ 0,&\tau_{0}\leq t<\tau_{1},\\\
(2m-1)\left[1-\mathrm{e}^{-\kappa(t-\tau_{1})}\cos\left(\sigma(t-\tau_{1})\right)\right],&t\geq\tau_{1}\end{cases}&c(\tau_{*})<0.\end{cases}$

As we have shown, this solution then also corresponds to the entropy solution
since it is in accordance with the characterization of the tangent cone in
Lemma 3.2. Like in the previous section, we can compare this with the
trajectory $\tilde{X}=\operatorname{P}_{\mathscr{K}}Y$ defined by the
projection formula (3.22), which in this case can be expressed as

$\tilde{X}(t,m)=\operatorname*{sgn}\left(m-\frac{1}{2}\right)\max\left\\{\operatorname*{sgn}\left(m-\frac{1}{2}\right)Y(t),0\right\\}=\left(m-\frac{1}{2}\right)\max\left\\{c(t),0\right\\}.$

See Figures 4 and 5 for illustrations of the undamped and damped cases
respectively. Note that for our initial velocity, $\lambda>0$ needs to be
rather small to have concentration, meaning the repulsive force is quite weak
and we see in Figures 4(b) and 5(b) that the mass stays concentrated for a
while before spreading out again.

(a) $\lambda=0.7$

(b) $\lambda=0.5$

Figure 4. The trajectories $X(t,\theta_{i})$,
$\tilde{X}(t,\theta_{i})=(\operatorname{P}_{\mathscr{K}}Y)(t,\theta_{i})$ and
$Y(t,\theta_{i})$ for $\theta_{i}=\frac{i}{20}$, $i\in\\{0,1,\dots,20\\}$,
$\kappa=0$ and sub- (a) and supercritical (b) values of $\lambda$.

(a) $\lambda=0.7$

(b) $\lambda=0.5$

Figure 5. The trajectories $X(t,\theta_{i})$,
$\tilde{X}(t,\theta_{i})=(\operatorname{P}_{\mathscr{K}}Y)(t,\theta_{i})$ and
$Y(t,\theta_{i})$ for $\theta_{i}=\frac{i}{20}$, $i\in\\{0,1,\dots,20\\}$,
$\kappa=0.05$, and sub- (a) and supercritical (b) values of $\lambda$.

####  7.2.3. A more involved example with particles

The choice of initial data $X^{0}$ and $V^{0}$ in the previous example
simplified the analysis, since the convexity properties of the flux function
$m\mapsto\mathcal{U}(t,m)$, or equivalently, the monotonicity properties of
the prescribed velocity $m\mapsto U(t,m)$ depended only on time $t$. The
smoothness of $V^{0}$ then meant that the inter-particle forces had to be
rather weak, i.e., $\lambda$ small enough, for concentration to happen.
Choosing instead initial data as in the repulsive Euler–Poisson case before,
namely $X^{0}(m)=m-\frac{1}{2}$ and
$V^{0}(m)=-\operatorname*{sgn}(m-\frac{1}{2})$, the discontinuous $V^{0}$
ensures that mass will immediately concentrate near the origin. However, in
this case, the convexity properties will depend on both $t$ and $m$, causing
far more complex dynamics where mass will concentrate also outside the origin.

Therefore, to simplify the analysis we will consider a particle approximation
of this problem, where the uniformly distributed mass on
$(-\frac{1}{2},\frac{1}{2})$ is first divided into four intervals
$(\frac{i-3}{4},\frac{i-2}{4})$ for $i\in\\{1,2,3,4\\}$ for then to be
concentrated at the midpoints $x_{i}^{0}=\frac{2i-5}{8}$ with
$m_{i}=\frac{1}{4}$. Moreover we define the cumulative quantity
$\theta_{i}=\frac{i}{4}$ for $i\in\\{0,\dots,5\\}$, where we note the
relations $x_{i}^{0}=\frac{1}{2}\left(\theta_{i-1}+\theta_{i}-1\right)$ and
$v_{i}^{0}=-\operatorname*{sgn}(\theta_{i-1}+\theta_{i}-1)$ for
$i\in\\{1,2,3,4\\}$. As in the previous examples, the symmetry of the initial
data ensures that the center of mass is always located at the origin,
simplifying our equations, and so the initial particle dynamics are given by

$\dot{x}_{i}=v_{i},\qquad\dot{v}_{i}=-2\kappa
v_{i}+\lambda^{2}\left(\theta_{i-1}+\theta_{i}-1-x_{i}\right).$

This we may integrate to find the “smooth” trajectories

$x_{i}(t)=x_{i}^{0}\left[2-\mathrm{e}^{-\kappa t}\left(\cos(\sigma
t)+\frac{\kappa}{\sigma}\sin(\sigma
t)\right)\right]+v_{i}^{0}\mathrm{e}^{-\kappa t}\frac{\sin(\sigma
t)}{\sigma}.$

Initially we see that the inter-particle repulsion dominates the confinement
force which attracts the particles to the center of mass at the origin. For
$\lambda>\kappa>0$ small enough, the repulsive and damping forces will be too
weak to stop the particles from colliding. By symmetry, particles $i=2,3$ will
collide first, at $t=\tau_{1}^{-}$. Now these particles will be stuck at the
origin for some time, where their prescribed velocities satisfy

$\dot{u}_{i}+2\kappa
u_{i}=\lambda^{2}\left(\theta_{i-1}+\theta_{i}-1\right),\qquad
u_{i}(\tau_{1}^{-})=v_{i}(\tau_{1}^{-}),$

which for $t\geq\tau_{1}^{-}$ integrates to

$u_{i}(t)=\mathrm{e}^{-2\kappa(t-\tau_{1}^{-})}v_{i}(\tau_{1}^{-})+\lambda^{2}\left(\theta_{i-1}+\theta_{i}-1\right)\frac{1-\mathrm{e}^{-2\kappa(t-\tau_{1}^{-})}}{2\kappa},\qquad
i\in\\{2,3\\}.$

Since the particles collided at $t=\tau_{1}^{-}$, we know that
$v_{2}(\tau_{1}^{-})>v_{3}(\tau_{1}^{-})=-v_{2}(\tau_{1}^{-})$, and so
$v_{i}(t)=(\operatorname{P}_{T_{\mathbf{x}(t)}\mathbb{K}^{n}}\mathbf{u})_{i}=0$
for $i=2,3$ until $u_{3}(t)\geq u_{2}(t)$. However, the second pair of
particles hits the origin at time $t=\tau_{1}^{+}>\tau_{1}^{-}$ before this
can happen. At this time we will have the chain of inequalities

$u_{1}(\tau_{1}^{+})>u_{2}(\tau_{1}^{+})>u_{3}(\tau_{1}^{+})>u_{4}(\tau_{1}^{+}),$

where by symmetry $u_{4}(\tau_{1}^{+})=-u_{1}(\tau_{1}^{+})$ and
$u_{3}(\tau_{1}^{+})=-u_{2}(\tau_{1}^{+})$. Analogously, the prescribed
velocity of the second particle pair for $t\geq\tau_{1}^{+}$ is

$u_{i}(t)=\mathrm{e}^{-2\kappa(t-\tau_{1}^{+})}v_{i}(\tau_{1}^{+})+\lambda^{2}\left(\theta_{i-1}+\theta_{i}-1\right)\frac{1-\mathrm{e}^{-2\kappa(t-\tau_{1}^{+})}}{2\kappa},\qquad
i\in\\{1,4\\}.$

The Oleĭnik E-condition will be violated when $u_{1}(t)$ changes sign at the
time $t=\tau_{2}^{+}$, and then
$v_{i}(t)=(\operatorname{P}_{T_{\mathbf{x}(t)}\mathbb{K}^{n}}\mathbf{u}(t))_{i}=u_{i}(t)$
for $i=1,4$, releasing them from the origin with initial data
$x_{i}(\tau_{2}^{+})=0$ and $v_{i}(\tau_{2}^{+})=0$, yielding trajectories

$x_{i}(t)=\left(\theta_{i-1}+\theta_{i}-1\right)\left[1-\mathrm{e}^{-\kappa(t-\tau_{2}^{+})}\left(\cos(\sigma(t-\tau_{2}^{+}))+\frac{\kappa}{\sigma}\sin(\sigma(t-\tau_{2}^{+}))\right)\right]$

At the later time $t=\tau_{2}^{-}>\tau_{2}^{+}$, $u_{2}(t)$ will change sign,
and particles $i=2,3$ will be released with corresponding trajectories

$x_{i}(t)=\left(\theta_{i-1}+\theta_{i}-1\right)\left[1-\mathrm{e}^{-\kappa(t-\tau_{2}^{-})}\left(\cos(\sigma(t-\tau_{2}^{-}))+\frac{\kappa}{\sigma}\sin(\sigma(t-\tau_{2}^{-}))\right)\right]$

For weak enough damping, at a time $t=\tau_{3}$ there will be a new collision
between particles $i=1,2$ and $i=3,4$ respectively. Since there was a
collision, one has that $v_{i}(\tau_{3})>v_{i+1}(\tau_{3})$ for $i=1,3$, and
the particles will stick together until the Oleĭnik E-condition no longer
holds, i.e., when $u_{i}(t)\leq u_{i+1}(t)$, which happens at $t=\tau_{4}$. In
the interval $t\in[\tau_{3},\tau_{4}]$, the pairs of particles will in fact
move along their centers of mass, which is the average of the two previous
trajectories. In the damped case, $t=\tau_{4}$ is the last time the particles
are in contact, see Figure 6. Indeed, for $t\geq\tau_{4}$, the trajectory of
the $i$th particle for $i\in\\{1,2,3,4\\}$ is given by

$\displaystyle x_{i}(t)$
$\displaystyle=x_{i}(\tau_{4})+\left(\theta_{i-1}+\theta_{i}-1-x_{i}(\tau_{4})\right)\left[1-\mathrm{e}^{-\kappa(t-\tau_{4})}\left(\cos(\sigma(t-\tau_{4}))+\frac{\kappa}{\sigma}\sin(\sigma(t-\tau_{4}))\right)\right]$
$\displaystyle\quad+v_{i}(\tau_{4})\mathrm{e}^{-\kappa(t-\tau_{4})}\frac{\sin(\sigma(t-\tau))}{\sigma}.$

As we can see, for $\kappa>0$, the particle trajectories will asymptotically
tend to the equilibrium state

$\lim\limits_{t\to\infty}x_{i}(t)\coloneqq\theta_{i-1}+\theta_{i}-1,\qquad\lim\limits_{t\to\infty}v_{i}(t)=0,$

which is the four-particle approximation of $\rho=\frac{1}{2}\chi_{(-1,1)}$
and $v\equiv 0$. In the undamped case $\kappa=0$, the trajectories of
particles 1 and 2, and similarly 3 and 4, will still periodically meet at
times $t=\tau_{4}+2\pi n/\sigma$, but tangentially, so there is no loss of
energy, see Figure 7.

Figure 6. Poisson repulsion and quadratic confinement with $\lambda=0.6$ and
damping parameter $\kappa=0.05$ for $t\in[0,6\pi/\sigma]$. Figure 7. Poisson
repulsion and quadratic confinement with $\lambda=0.6$ and no damping for
$t\in[0,6\pi/\sigma]$.

The undamped case should be compared to the numerical example considered in
[7, Section 1.4] where a time-discrete scheme is applied to the repulsive
Euler–Poisson system with neutralizing background and periodic domain; in this
case the resulting equations correspond to the repulsive Poisson and quadratic
confinement potential above. There it is observed that after an initial period
of trajectories colliding and dissipating energy, at some point the
trajectories only touch tangentially and no more energy is dissipated.

####  7.2.4. Asymptotic behavior of the damped system

In our previous examples we saw that the presence of damping drove our system
to a steady equilibrium state in the form of a uniform probability
distribution on a finite interval. This asymptotic behavior was also shown to
happen at an exponential rate in [14] for solutions which remain classical. We
will show that this is also the case for our Lagrangian, or entropic,
solutions. Let us as above consider $\lambda>\kappa>0$, this is not strictly
necessary and all the other cases can also be dealt with in the same manner.
In Lagrangian variables, the system reads

(7.6) $\dot{X}(t)+\partial I_{\mathscr{K}}(X(t))\ni
U(t),\qquad\dot{U}(t)=-2\kappa
U(t)+\lambda^{2}\left(2m-1+\int_{0}^{1}X(t,\omega)\mathop{}\\!\mathrm{d}\omega-X(t)\right).$

Observe that $\int_{0}^{1}X(t,\omega)\mathop{}\\!\mathrm{d}\omega$ is simply
the center of mass $\bar{x}(t)$ expressed in Lagrangian variables, and in our
case $\beta=0$ we find from (7.4)

(7.7)
$\bar{x}_{\infty}\coloneqq\lim\limits_{t\to\infty}\bar{x}(t)=\bar{x}(0)+\frac{\bar{v}(0)}{2\kappa},\qquad\bar{v}_{\infty}=\lim\limits_{t\to\infty}\bar{v}(t)=0.$

This suggests the steady equilibrium state with $\dot{X}=V=0$ and

(7.8)
$X=2m-1+\int_{0}^{1}X^{0}(m)\mathop{}\\!\mathrm{d}m+\frac{1}{2\kappa}\int_{0}^{1}V^{0}(m)\mathop{}\\!\mathrm{d}m=2m-1+\int_{\mathbb{R}}x\mathop{}\\!\mathrm{d}\rho_{0}(x)+\frac{1}{2\kappa}\int_{\mathbb{R}}v_{0}(x)\mathop{}\\!\mathrm{d}\rho_{0}(x),$

which, with appropriate rescaling of $\rho$ and parameter choices, agrees with
the asymptotics for global-in-time classical solutions in [14, Theorem 4.1].
We will show that the damping term indeed forces any initial data
$X^{0}\in\mathscr{K}$, $V^{0}\in L^{2}(0,1)$ to this equilibrium state
exponentially fast.

######  Theorem 7.3 (Asymptotic behavior).

Let $(\rho,v)$ be the solution of (1.6) corresponding to the Lagrangian, or
entropy, solution for initial data $\rho_{0}\in\mathscr{P}_{2}(\mathbb{R})$
and $v_{0}\in L^{2}(\mathbb{R},\rho_{0})$. Then the asymptotic solution is
given by

$\rho_{\infty}(x)\coloneqq\lim\limits_{t\to\infty}\rho(t,x)=\frac{1}{2}\chi_{(-1,1)}(x-\bar{x}_{\infty}),\qquad
v_{\infty}(x)\coloneqq\lim\limits_{t\to\infty}v(t,x)\equiv 0$

with $\bar{x}_{\infty}$ defined in (7.7). Moreover, this convergence is
exponentially fast in the $\mathscr{T}_{2}$-metric (3.25) at a rate
proportional to the damping factor $\gamma=2\kappa$.

###### Proof.

Let us start by considering the evolution of the energy, which is the sum of a
kinetic and a potential part

$\frac{1}{2}\int_{\mathbb{R}}v(t,x)^{2}\mathop{}\\!\mathrm{d}\rho(t,x)+\frac{\lambda^{2}}{2}\iint_{\mathbb{R}^{2}}\left(\frac{1}{2}(x-y)^{2}-\lvert
x-y\rvert\right)\mathop{}\\!\mathrm{d}\rho(t,y)\mathop{}\\!\mathrm{d}\rho(t,x),$

and which in Lagrangian variables reads

$E(t)\coloneqq\frac{1}{2}\int_{0}^{1}V(t,m)^{2}\mathop{}\\!\mathrm{d}m+\frac{\lambda^{2}}{2}\iint_{[0,1]^{2}}\left(\frac{1}{2}(X(t,m)-X(t,\omega))^{2}-\lvert
X(t,m)-X(t,\omega)\rvert\right)\mathop{}\\!\mathrm{d}\omega\mathop{}\\!\mathrm{d}m.$

However, there is the issue that we do not really have an evolution equation
for the velocity $V$. Fortunately, we have the relation
$V=\operatorname{P}_{T_{X}\mathscr{K}}U$, meaning $\lVert
V\rVert_{L^{2}(0,1)}\leq\lVert U\rVert_{L^{2}(0,1)}$, which allows us to bound
$E(t)$ from above by

$\hat{E}(t)\coloneqq\frac{1}{2}\int_{0}^{1}U(t,m)^{2}\mathop{}\\!\mathrm{d}m+\frac{\lambda^{2}}{2}\iint_{[0,1]^{2}}\left(\frac{1}{2}(X(t,m)-X(t,\omega))^{2}-\lvert
X(t,m)-X(t,\omega)\rvert\right)\mathop{}\\!\mathrm{d}\omega\mathop{}\\!\mathrm{d}m,$

and we have an evolution equation for $U$. Differentiating this quantity we
find

$\displaystyle\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\hat{E}(t)$
$\displaystyle=\int_{0}^{1}U(t,m)\left[-2\kappa
U(t,m)+\lambda^{2}\left(2m-1+\int_{0}^{1}X(t,\omega)\mathop{}\\!\mathrm{d}\omega-X(t,m)\right)\right]\mathop{}\\!\mathrm{d}m$
$\displaystyle\quad+\frac{\lambda^{2}}{2}\int_{0}^{1}\int_{0}^{1}\left(X(t,m)-X(t,\omega)-\operatorname*{sgn}\left(X(t,m)-X(t,\omega)\right)\right)(V(t,m)-V(t,\omega))\mathop{}\\!\mathrm{d}\omega\mathop{}\\!\mathrm{d}m$
$\displaystyle=-2\kappa\int_{0}^{1}U(t,m)^{2}\mathop{}\\!\mathrm{d}m+\lambda^{2}\int_{0}^{1}U(t,m)\left(2m-1+\int_{0}^{1}X(t,\omega)\mathop{}\\!\mathrm{d}\omega-X(t,m)\right)\mathop{}\\!\mathrm{d}m$
$\displaystyle\quad+\lambda^{2}\int_{0}^{1}V(t,m)\left(X(t,m)-\int_{0}^{1}X(t,\omega)\mathop{}\\!\mathrm{d}\omega-\int_{0}^{1}\operatorname*{sgn}\left(X(t,m)-X(t,\omega)\right)\mathop{}\\!\mathrm{d}\omega\right)\mathop{}\\!\mathrm{d}m$
$\displaystyle=-2\kappa\int_{0}^{1}U(t,m)^{2}\mathop{}\\!\mathrm{d}m-\lambda^{2}\int_{0}^{1}\left(U(t,m)-V(t,m)\right)\left(X(t,m)-\int_{0}^{1}X(t,\omega)\mathop{}\\!\mathrm{d}\omega\right)\mathop{}\\!\mathrm{d}m$
$\displaystyle\quad+\lambda^{2}\int_{0}^{1}V(t,m)\left(2m-1-\operatorname{P}_{\mathscr{H}_{X(t)}}(2m-1)\right)\mathop{}\\!\mathrm{d}m+\lambda^{2}\int_{0}^{1}\left(U(t,m)-V(t,m)\right)\left(2m-1\right)\mathop{}\\!\mathrm{d}m,$

where we again have used (7.3). In the final expression, the second and third
integral vanish due to Lemma 3.7 and (5.8); the second because of $U-V\in
N_{X}\mathscr{K}\subset\mathscr{H}_{X}^{\perp}$, while the third vanishes for
a.e. $t$ since $V(t)\in\mathscr{H}_{X(t)}$ for a.e. $t$. We note that $U-V\in
N_{X}\mathscr{K}$ while $2m-1$ is nondecreasing, which allow us to combine the
characterizations (3.6) and (3.9) to deduce that the final term is
nonpositive, yielding

$\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\hat{E}(t)\leq-2\kappa\int_{0}^{1}U(t,m)^{2}\mathop{}\\!\mathrm{d}m.$

From this we deduce $E(t)\leq\hat{E}(t)\leq\hat{E}(0)=E(0)$. Moreover, this
implies that $V$ decays exponentially: indeed,

$\displaystyle\frac{1}{2}\int_{0}^{1}U(t,m)^{2}\mathop{}\\!\mathrm{d}m\leq\hat{E}(t)\leq\hat{E}(0)-2\kappa\int_{0}^{t}\int_{0}^{1}U(s,m)^{2}\mathop{}\\!\mathrm{d}m\mathop{}\\!\mathrm{d}s$

which by Grönwall’s inequality implies

$\lVert V(t)\rVert_{L^{2}(0,1)}^{2}\leq\lVert U(t)\rVert_{L^{2}(0,1)}^{2}\leq
2\mathrm{e}^{-4\kappa t}\hat{E}(0)=2\mathrm{e}^{-4\kappa t}E(0).$

This suggests that a damped system, i.e., with $\kappa>0$, will approach an
equilibrium state at an exponential rate proportional to $\kappa$. To simplify
notation we introduce

$\displaystyle\tilde{X}(t,m)$ $\displaystyle=X(t,m)-2m+1-\bar{x}(t),$
$\displaystyle\tilde{V}(t,m)$
$\displaystyle=V(t,m)-\bar{v}(t),\quad\tilde{U}(t,m)=U(t,m)-\bar{v}(t),$

where we note that the center of mass $\bar{x}(t)=\bar{X}(t)$ and total
momentum $\bar{v}(t)=\bar{V}(t)$ satisfy

$\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\bar{X}(t)=\bar{V}(t),\qquad\bar{v}(t)=\int_{0}^{1}V(t,\omega)\mathop{}\\!\mathrm{d}\omega=\int_{0}^{1}U(t,\omega)\mathop{}\\!\mathrm{d}\omega.$

Here $(X-\tilde{X},V-\tilde{V})$ corresponds to a uniform distribution
centered at the center of mass for the system, and according to (7.7) this
converges uniformly pointwise in $\Omega$ to the proposed steady state (7.8).
From the fact that $-\tilde{X},\tilde{V}\in T_{X}\mathscr{K}$, (3.6), (5.8),
$\tilde{U}-\tilde{V}=U-V\in N_{X}\mathscr{K}$ and
$\int_{0}^{1}\tilde{X}\mathop{}\\!\mathrm{d}m=0=\int_{0}^{1}\tilde{U}\mathop{}\\!\mathrm{d}m$
we derive the following relations,

$\displaystyle\int_{0}^{1}\tilde{X}\tilde{V}\mathop{}\\!\mathrm{d}m$
$\displaystyle\leq\int_{0}^{1}\tilde{X}\tilde{U}\mathop{}\\!\mathrm{d}m=\int_{0}^{1}\tilde{X}U\mathop{}\\!\mathrm{d}m,$
$\displaystyle\int_{0}^{1}\tilde{V}\tilde{U}\mathop{}\\!\mathrm{d}m$
$\displaystyle=\int_{0}^{1}\tilde{V}^{2}\mathop{}\\!\mathrm{d}m=\int_{0}^{1}V^{2}\mathop{}\\!\mathrm{d}m-\bar{v}^{2}\leq\int_{0}^{1}U^{2}\mathop{}\\!\mathrm{d}m-\bar{v}^{2}=\int_{0}^{1}\tilde{U}^{2}\mathop{}\\!\mathrm{d}m=\int_{0}^{1}\tilde{U}U\mathop{}\\!\mathrm{d}m.$

The first identity in the second line holds for a.e. $t$ since
$\tilde{V}(t)\in\mathscr{H}_{X(t)}$ for a.e. $t$, alternatively it holds for
all times with an inequality since $\tilde{V}\in T_{X}\mathscr{K}$. From these
relations, the evolution equations for $(\tilde{X},\tilde{V})$ and (7.6), we
then, omitting the subscript $L^{2}(0,1)$ from norms and inner products for
brevity, obtain

$\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\left(\frac{c_{1}}{2}\lVert\tilde{X}\rVert^{2}+c_{2}\langle\tilde{X},\tilde{U}\rangle+\frac{1}{2}\lVert\tilde{U}\rVert^{2}\right)\leq-
c_{2}\lambda^{2}\lVert\tilde{X}\rVert^{2}+(c_{1}-2\kappa
c_{2}-\lambda^{2})\langle\tilde{X},\tilde{U}\rangle+(c_{2}-2\kappa)\lVert\tilde{U}\rVert^{2}$

for $c_{1},c_{2}\geq 0$. Choosing $c_{1}=\lambda^{2}$ and $c_{2}=\kappa$, we
can apply Grönwall’s inequality to find

$\lambda^{2}\lVert\tilde{X}(t)\rVert^{2}+2\kappa\langle\tilde{X}(t),\tilde{U}(t)\rangle+\lVert\tilde{U}(t)\rVert^{2}\leq\left(\lambda^{2}\lVert\tilde{X}(0)\rVert^{2}+2\kappa\langle\tilde{X}(0),\tilde{U}(0)\rangle+\lVert\tilde{U}(0)\rVert^{2}\right)e^{-2\kappa
t},$

or, recalling $\lambda^{2}=\kappa^{2}+\sigma^{2}$,

$\sigma^{2}\lVert\tilde{X}(t)\rVert^{2}+\lVert\tilde{U}(t)+\kappa\tilde{X}(t)\rVert^{2}\leq\left(\sigma^{2}\lVert\tilde{X}(0)\rVert^{2}+\lVert\tilde{U}(0)+\kappa\tilde{X}(0)\rVert^{2}\right)\mathrm{e}^{-2\kappa
t}.$

That is, both $\tilde{X}$ and $\tilde{V}$ tend exponentially fast to zero in
the $L^{2}$-norm with rates proportional to $\gamma=2\kappa$. Then, as a
Lagrangian solution satisfies $V(t,m)=v(t,X(t,m))$ for a.e. $t$, we use (3.26)
to conclude. ∎

####  7.2.5. Lyapunov stability of undamped system

In the undamped cases of our examples, we have seen that the solutions did not
tend to a steady state. One could then instead ask, if two solutions are close
in some sense, will they remain close?

To this end, take two solutions $(X_{1},V_{1})$ and $(X_{2},V_{2})$ with
corresponding prescribed velocities $U_{1}$ and $U_{2}$. We further introduce
the averaged quantities $\bar{x}_{i}$, $\bar{v}_{i}$, as well as
$\hat{X}_{i}\coloneqq X_{i}-\bar{x}_{i}$, $\hat{V}_{i}\coloneqq
V_{i}-\bar{v}_{i}$, $\hat{U}_{i}\coloneqq U_{i}-\bar{v}_{i}$ for
$i\in\\{1,2\\}$. From the equations we then get

$\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\left(\hat{X}_{1}-\hat{X}_{2}\right)=\hat{V}_{1}-\hat{V}_{2},\qquad\frac{\mathop{}\\!\mathrm{d}}{\mathop{}\\!\mathrm{d}t}\left(\hat{U}_{1}-\hat{U}_{2}\right)=-\lambda^{2}\left(\hat{X}_{1}-\hat{X}_{2}\right),$

leading to

$\displaystyle\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\left(\lambda^{2}\lVert\hat{X}_{1}-\hat{X}_{2}\rVert^{2}+\lVert\hat{U}_{1}-\hat{U}_{2}\rVert^{2}\right)$
$\displaystyle=2\lambda^{2}\langle\hat{X}_{1}-\hat{X}_{2},\hat{V}_{1}-\hat{V}_{2}\rangle-2\lambda^{2}\langle\hat{U}_{1}-\hat{U}_{2},\hat{X}_{1}-\hat{X}_{2}\rangle$
$\displaystyle=-2\lambda^{2}\langle\hat{X}_{1}-\hat{X}_{2},U_{1}-V_{1}-(U_{2}-V_{2})\rangle$
$\displaystyle=2\lambda^{2}\left(\langle\hat{X}_{1},U_{2}-V_{2}\rangle+\langle\hat{X}_{2},U_{1}-V_{1}\rangle-\sum_{i=1}^{2}\langle\hat{X}_{i},U_{i}-V_{i}\rangle\right).$

Now, since $\hat{X}_{i}\in\mathscr{H}_{X_{i}}$ and $U_{i}-V_{i}\in
N_{X_{i}}\mathscr{K}\subset\mathscr{H}_{X_{i}}^{\perp}$, we have
$\langle\hat{X}_{i},U_{i}-V_{i}\rangle=0$. For the first cross-term, we note
that $\hat{X}_{1}(t,\cdot)$ is nondecreasing, so it belongs to the tangent
cone $T_{X_{2}}\mathscr{K}$ by Lemma 3.9, and then it follows from (3.6) that
$\langle\hat{X}_{1},U_{2}-V_{2}\rangle\leq 0$. The same argument applies to
the second cross-term, and we have shown

$\frac{\mathop{}\\!\mathrm{d}^{+}}{\mathop{}\\!\mathrm{d}t}\left(\lambda^{2}\lVert\hat{X}_{1}-\hat{X}_{2}\rVert^{2}+\lVert\hat{U}_{1}-\hat{U}_{2}\rVert^{2}\right)\leq
0,$

that is,

$\lambda^{2}\lVert\hat{X}_{1}(t)-\hat{X}_{2}(t)\rVert^{2}+\lVert\hat{U}_{1}(t)-\hat{U}_{2}(t)\rVert^{2}\leq\lambda^{2}\lVert\hat{X}_{1}(0)-\hat{X}_{2}(0)\rVert^{2}+\lVert\hat{V}_{1}(0)-\hat{V}_{2}(0)\rVert^{2}.$

In particular, comparing $(X_{1}(t),V_{1}(t))=(X(t),V(t))$ with the
‘equilibrium state’
$(\tilde{X}_{2}(t),\tilde{V}_{2}(t))=(2m-1+\bar{x}(t),\bar{v}(t))$, writing
$\tilde{X}=\hat{X}-(2m-1)$, we have the stability result

$\lambda^{2}\lVert\tilde{X}(t)\rVert^{2}+\lVert\hat{V}(t)\rVert^{2}\leq\lambda^{2}\lVert\tilde{X}(t)\rVert^{2}+\lVert\hat{U}(t)\rVert^{2}\leq\lambda^{2}\lVert\tilde{X}(0)\rVert^{2}+\lVert\hat{V}(0)\rVert^{2}.$

In particular, for this undamped case we have $\bar{v}(t)\equiv\bar{v}(0)$
such that $\bar{x}(t)=\bar{x}(0)+t\bar{v}(0)$, and so the above result can be
restated as

$\lambda^{2}\lVert\tilde{X}(t)\rVert^{2}+\lVert
V(t)\rVert^{2}\leq\lambda^{2}\lVert\tilde{X}(t)\rVert^{2}+\lVert
U(t)\rVert^{2}\leq\lambda^{2}\lVert\tilde{X}(0)\rVert^{2}+\lVert
V(0)\rVert^{2}.$

#### Acknowledgments

JAC was supported by the Advanced Grant Nonlocal-CPD (Nonlocal PDEs for
Complex Particle Dynamics: Phase Transitions, Patterns and Synchronization) of
the European Research Council Executive Agency (ERC) under the European
Union’s Horizon 2020 research and innovation programme (grant agreement No.
883363) and the EPSRC grant number EP/V051121/1. STG was supported by the
Research Council of Norway through grant no. 286822, Wave Phenomena and
Stability — a Shocking Combination (WaPheS), and the Swedish Research Council
under grant no. 2021-06594 while in residence at Institut Mittag-Leffler in
Djursholm, Sweden during the fall semester of 2023.

## References

  * [1] L. Ambrosio, N. Gigli, and G. Savaré.  Gradient flows in metric spaces and in the space of probability measures.  Lectures in Mathematics ETH Zürich. Birkhäuser Verlag, Basel, second edition, 2008.
  * [2] M. Bhatnagar, H. Liu, and C. Tan.  Critical thresholds in the Euler-Poisson-alignment system.  J. Differential Equations, 375:82–119, 2023.
  * [3] F. Bolley, Y. Brenier, and G. Loeper.  Contractive metrics for scalar conservation laws.  J. Hyperbolic Differ. Equ., 2(1):91–107, 2005.
  * [4] G. A. Bonaschi, J. A. Carrillo, M. Di Francesco, and M. A. Peletier.  Equivalence of gradient flows and entropy solutions for singular nonlocal interaction equations in 1D.  ESAIM Control Optim. Calc. Var., 21(2):414–441, 2015.
  * [5] F. Bouchut and F. James.  Duality solutions for pressureless gases, monotone scalar conservation laws, and uniqueness.  Comm. Partial Differential Equations, 24(11-12):2173–2189, 1999\.
  * [6] Y. Brenier.  $L^{2}$ formulation of multidimensional scalar conservation laws.  Arch. Ration. Mech. Anal., 193(1):1–19, 2009.
  * [7] Y. Brenier, W. Gangbo, G. Savaré, and M. Westdickenberg.  Sticky particle dynamics with interactions.  J. Math. Pures Appl. (9), 99(5):577–617, 2013.
  * [8] Y. Brenier and E. Grenier.  Sticky particles and scalar conservation laws.  SIAM J. Numer. Anal., 35(6):2317–2328, 1998.
  * [9] A. Bressan.  Hyperbolic systems of conservation laws, volume 20 of Oxford Lecture Series in Mathematics and its Applications.  Oxford University Press, Oxford, 2000.  The one-dimensional Cauchy problem.
  * [10] H. Brézis.  Opérateurs maximaux monotones et semi-groupes de contractions dans les espaces de Hilbert.  North-Holland Mathematics Studies, No. 5. North-Holland Publishing Co., Amsterdam-London; American Elsevier Publishing Co., Inc., New York, 1973\.
  * [11] J. Carrillo.  Entropy solutions for nonlinear degenerate problems.  Arch. Ration. Mech. Anal., 147(4):269–361, 1999.
  * [12] J. A. Carrillo, Y.-P. Choi, and S. P. Perez.  A review on attractive-repulsive hydrodynamics for consensus in collective behavior.  In Active particles. Vol. 1. Advances in theory, models, and applications, Model. Simul. Sci. Eng. Technol., pages 259–298. Birkhäuser/Springer, Cham, 2017.
  * [13] J. A. Carrillo, Y.-P. Choi, E. Tadmor, and C. Tan.  Critical thresholds in 1D Euler equations with non-local forces.  Math. Models Methods Appl. Sci., 26(1):185–206, 2016.
  * [14] J. A. Carrillo, Y.-P. Choi, and E. Zatorska.  On the pressureless damped Euler-Poisson equations with quadratic confinement: critical thresholds and large-time behavior.  Math. Models Methods Appl. Sci., 26(12):2311–2340, 2016.
  * [15] J. A. Carrillo, M. Di Francesco, and C. Lattanzio.  Contractivity of Wasserstein metrics and asymptotic profiles for scalar conservation laws.  J. Differential Equations, 231(2):425–458, 2006.
  * [16] J. A. Carrillo, M. Di Francesco, and C. Lattanzio.  Contractivity and asymptotics in Wasserstein metrics for viscous nonlinear scalar conservation laws.  Boll. Unione Mat. Ital. Sez. B Artic. Ric. Mat. (8), 10(2):277–292, 2007.
  * [17] F. Cavalletti, M. Sedjro, and M. Westdickenberg.  A simple proof of global existence for the 1D pressureless gas dynamics equations.  SIAM J. Math. Anal., 47(1):66–79, 2015.
  * [18] C. M. Dafermos.  Polygonal approximations of solutions of the initial value problem for a conservation law.  J. Math. Anal. Appl., 38:33–41, 1972.
  * [19] C. M. Dafermos.  Hyperbolic conservation laws in continuum physics, volume 325 of Grundlehren der mathematischen Wissenschaften [Fundamental Principles of Mathematical Sciences].  Springer-Verlag, Berlin, fourth edition, 2016.
  * [20] A. de La Fortelle.  A study on generalized inverses and increasing functions Part I: generalized inverses.  hal-01255512, 2016\.
  * [21] M. Di Francesco.  Scalar conservation laws seen as gradient flows: known results and new perspectives.  In Gradient flows: from theory to application, volume 54 of ESAIM Proc. Surveys, pages 18–44. EDP Sci., Les Ulis, 2016.
  * [22] W. E, Y. G. Rykov, and Y. G. Sinai.  Generalized variational principles, global weak solutions and behavior with random initial data for systems of conservation laws arising in adhesion particle dynamics.  Comm. Math. Phys., 177(2):349–380, 1996.
  * [23] S. Engelberg, H. Liu, and E. Tadmor.  Critical thresholds in Euler-Poisson equations.  volume 50, pages 109–157. 2001.  Dedicated to Professors Ciprian Foias and Roger Temam (Bloomington, IN, 2000).
  * [24] D. Golovaty and T. Nguyen.  On the existence, uniqueness and stability of entropy solutions to scalar conservation laws.  J. Differential Equations, 253(5):1341–1375, 2012.
  * [25] H. Holden and N. H. Risebro.  Front Tracking for Hyperbolic Conservation Laws, volume 152 of Applied Mathematical Sciences.  Springer, Heidelberg, second edition, 2015.
  * [26] F. Huang and G. Tang.  Global well-posedness of pressureless Euler-Poisson system in one dimension.  SSRN preprint ssrn:4539247, 2023.
  * [27] F. Huang and Z. Wang.  Well posedness for pressureless flow.  Comm. Math. Phys., 222(1):117–146, 2001.
  * [28] R. Hynd.  Lagrangian coordinates for the sticky particle system.  SIAM J. Math. Anal., 51(5):3769–3795, 2019.
  * [29] M. Josephy.  Composing functions of bounded variation.  Proc. Amer. Math. Soc., 83(2):354–356, 1981.
  * [30] K. H. Karlsen and M. Ohlberger.  A note on the uniqueness of entropy solutions of nonlinear degenerate parabolic equations.  J. Math. Anal. Appl., 275(1):439–458, 2002.
  * [31] S. N. Kružkov.  First order quasilinear equations with several independent variables.  Mat. Sb. (N.S.), 81 (123):228–255, 1970.
  * [32] P. G. LeFloch and S. Xiang.  Existence and uniqueness results for the pressureless Euler-Poisson system in one spatial variable.  Port. Math., 72(2-3):229–246, 2015.
  * [33] T. M. Leslie and C. Tan.  Sticky particle Cucker–Smale dynamics and the entropic selection principle for the 1D Euler-alignment system.  Comm. Partial Differential Equations, 48(5):753–791, 2023.
  * [34] K.-A. Lie.  Front tracking for one-dimensional quasilinear hyperbolic equations with variable coefficients.  Numer. Algorithms, 24(3):275–298, 2000.
  * [35] L. Natile and G. Savaré.  A Wasserstein approach to the one-dimensional sticky particle system.  SIAM J. Math. Anal., 41(4):1340–1365, 2009.
  * [36] T. Nguyen and A. Tudorascu.  Pressureless Euler/Euler-Poisson systems via adhesion dynamics and scalar conservation laws.  SIAM J. Math. Anal., 40(2):754–775, 2008.
  * [37] T. Nguyen and A. Tudorascu.  One-dimensional pressureless gas systems with/without viscosity.  Comm. Partial Differential Equations, 40(9):1619–1665, 2015.
  * [38] O. A. Oleĭnik.  Discontinuous solutions of non-linear differential equations.  Uspehi Mat. Nauk (N.S.), 12(3(75)):3–73, 1957.  Translated in Amer. Math. Soc. Transl. (2), 26:95–172, 1963\.
  * [39] O. A. Oleĭnik.  Uniqueness and stability of the generalized solution of the Cauchy problem for a quasi-linear equation.  Uspehi Mat. Nauk, 14(2(86)):165–170, 1959.  Translated in Amer. Math. Soc. Transl. (2), 33:285–290, 1963\.
  * [40] E. Y. Panov.  On generalized entropy solutions of the Cauchy problem for a first-order quasilinear equation in the class of locally integrable functions.  Izv. Ross. Akad. Nauk Ser. Mat., 66(6):91–136, 2002.
  * [41] F. Santambrogio.  Optimal transport for applied mathematicians, volume 87 of Progress in Nonlinear Differential Equations and their Applications.  Birkhäuser/Springer, Cham, 2015.  Calculus of variations, PDEs, and modeling.
  * [42] S. F. Shandarin and Y. B. Zel’dovich.  The large-scale structure of the universe: turbulence, intermittency, structures in a self-gravitating medium.  Rev. Modern Phys., 61(2):185–220, 1989.
  * [43] M. Suder and A. Tudorascu.  On the Lagrangian description and uniqueness for the one-dimensional pressureless Euler system.  Indiana Univ. Math. J., 72(2):785–819, 2023.
  * [44] E. Tadmor and C. Tan.  Critical thresholds in flocking hydrodynamics with non-local alignment.  Philos. Trans. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci., 372(2028):20130401, 22, 2014.
  * [45] E. Tadmor and D. Wei.  On the global regularity of subcritical Euler-Poisson equations with pressure.  J. Eur. Math. Soc. (JEMS), 10(3):757–769, 2008.
  * [46] E. Tadmor and D. Wei.  A variational representation of weak solutions for the pressureless Euler-Poisson equations.  arXiv preprint: arXiv:1102.5579, 2011.
  * [47] C. Villani.  Topics in optimal transportation, volume 58 of Graduate Studies in Mathematics.  American Mathematical Society, Providence, RI, 2003.
  * [48] A. I. Vol’pert.  Spaces ${\rm BV}$ and quasilinear equations.  Mat. Sb. (N.S.), 73 (115):255–302, 1967.
  * [49] Y. B. Zel’dovich.  Gravitational instability: An approximate theory for large density perturbations.  Astronom. and Astrophys., 5:84–89, 1970.