\section{\textbf{Introduction}}
\label{Sec:1}

In most sports involving the race against other participants for time,
drag plays a crucial role.
It is defined as the aerodynamic force acting in an axial direction on
the athlete as he or she moves against the surrounding fluid, such as water or air.
As a result, it has become common practice both in sports and motorsports to
benefit from the presence of other athletes or vehicles to alleviate the resistance of the fluid
by drafting in the associated wake region (i.e. \emph{slipstreaming}).
In the wake there is not just drastically lower fluid velocity, but also a negative
pressure coefficient leading to suction, which additionally benefits the trailing competitor.

While in motorsports with its high speeds and considerable research efforts,
drafting has long been investigated (e.g. \citet{katz1995race, romberg1971aerodynamics}), it is less researched in other sports.
However, \citet{rundell1996effects} investigated drafting during speed skating and
showed large positive effects on metabolic activity, heart rate, and lactate response.
Also in swimming, drafting was associated with considerable energy saving and
the optimal drafting distance was established at $0$-$\SI{0.5}{\meter}$ behind the leading swimmer,
causing power savings of up to $31\%$ \citep{bassett1991metabolic, chatard2003drafting}.
Furthermore, cycling is a sport of large benefits through aerodynamic drafting and riding in
formations, such as tightly packed pelotons \citep{edwards2007aerodynamic, blocken2018aerodynamic, malizia2020bicycle}.
Apart from testing in wind tunnels, typically applying similarity principles through scaled
models at similar Reynolds numbers, computational fluid dynamics (CFD) has become
a widely used method to accurately assess drag distributions across athletes.
\citet{blocken2013cfd} examined the effect of two drafting cyclists at various upper-body
positions through CFD and found drag reductions of up
to $27.1\%$ for the trailing cyclist, while also the leading cyclist's drag was decreased by up to $2.6\%$.
Moreover, drag reductions of $42\%$ to $48\%$ have been measured for drafting during cycling in the velodrome \citep{fitton2018impact}.
\citet{fuss2018slipstreaming} recently suggested an analytical model for slipstreaming in gravity-powered sports
such as ski cross and found significant advantages in gained velocity and glide distance with respect to a leading skier,
especially in a tucked body position.

Although running happens typically at fairly low speeds,
drafting shows also here a measurable impact on aerodynamic forces.
In a major study on the subject, \citet{pugh1970oxygen} demonstrated drag savings of up
to $80\%$  when running at a middle-distance speed of $\SI{4.46}{\meter\per\second}$ against a slight head wind.
More recently, particularly in light of attempts to break the $2$-hour barrier of the marathon distance (e.g. \emph{INEOS 1:59 Challenge}, \emph{Nike Breaking2}),
increasing research efforts were directed to minimize drag for a main runner in formations.
\citet{hoogkamer2018modeling} used a reduced-order model to establish a sustainable velocity
of $\SI{5.93}{\meter\per\second}$ using cooperative drafting in a four-runner team.
\citet{polidori2020numerical} applied CFD to compute the drag and power savings of Kenenisa Bekele,
while running the Berlin marathon and using cooperative drafting.
It was shown that up to $57.3\%$ of aerodynamic power were gained in the ideal case of running directly
behind a front pacer at a distance of $\SI{1.3}{\meter}$ when compared to the case of running alone.

The current study aims at shedding more light on the aerodynamic effects of running alone and in formations
by using a realistic model of a female runner and a validated CFD methodology.
Four different formations are chosen to reflect various running scenarios and their implications on the
drag on the main runner as well as her pacers are discussed.
The effect of the resulting drag reduction on the metabolic rate and running economy
is computed by using a mechanical power model of running.
As a result, possible performance predictions are made in terms of improvements of
running speed and time.
Finally, by applying the conclusions from this study, direct implications for an optimized drafting strategy are achieved.

\section{\textbf{Methods}}
\label{Sec:2}

\subsection{\textbf{Numerical method}}
\label{Sec:2A}

For the computation of the flow field variables and the resulting aerodynamic
forces on the runner, computational fluid dynamics (CFD) using the finite-volume method is applied.
The Navier-Stokes equations describing the conservation of momentum are
discretized in its Reynolds-averaged (RANS) form.

Turbulence models have an important effect on aerodynamic force computations, mainly because of their modelling of momentum transfer near surfaces and in boundary layers. This in turn influences the locations of flow separation, which is crucial for correct drag predictions.
Here, turbulence is treated with the $k$-$\omega$ model in its shear-stress-transport (SST) formulation according to \citet{menter1994two}.
The model has been repeatedly shown to be suitable for external aerodynamics
and to give superior performance to the Spalart-Allmaras and the $k$-$\epsilon$ models
in flow scenarios with adverse pressure gradients and free shear layers,
which do occur in our setup \citep{menter1992improved,hellsten1998some}.
\citet{defraeye2010computational} tested the performance of various turbulence
treatments such as the most common RANS models and large eddy simulation
with respect to drag predictions for cycling aerodynamics.
They found the $k$-$\omega$ SST model to perform best in comparison
with experimental data (i.e. force and moment areas) from wind tunnel measurements,
with average discrepancies of $6\%$ and consistently below $11\%$.

\begin{figure*}[htbp]
\begin{centering}
\begin{overpic}[width=0.90\textwidth]{Figure1.pdf}
\end{overpic}
\par
\end{centering}
\caption{Geometry of a formation of runners surrounded by the computational domain
and the imposed boundary conditions (a). The axial spacing of $1.2\,m$ and the lateral spacing
of $0.7\,m$ between the runners are indicated (b).}
\label{Fig:2A-1}
\end{figure*}

Moreover, \citet{crouch2016dynamic} has shown for cycling aerodynamics that the effect of
the dynamical motion of the limbs on instantaneous drag is only minor.
Additionally, wind tunnel experiments by \citet{inoue2016wind} revealed an approximately
$10\%$ increase in drag values of solo running when using a moving-belt system as compared with a stationary setup.
In their measurements the formation of unrealistic ground boundary layers was thus inhibited, as it is the case in this study.
We use a stationary model of an athlete together with a moving ground boundary,
which was found to be sufficient for the prediction of the global aerodynamic forces.

Furthermore, discretization in this study is performed with a second-order scheme
for the spatial derivatives. Pressure and velocity are calculated in a coupled approach.
The fluid is defined as air at ambient temperature of $T=\SI{293.15}{K}=\SI{20}{\celsius}$
and a kinematic viscosity of $\nu=\SI{1.516e-5}{\meter\squared\per\second}$.
Additionally, the velocity inlet of the domain pictured in Fig. \ref{Fig:2A-1}(a)
imposes an axial flow velocity at the chosen running speeds, while the pressure outlet
applies an ambient pressure of $p= \SI{101325}{\pascal}=\SI{1}{\atm}$.
The ground is defined as wall moving at the running speed, which ensures
that no boundary layers form at the bottom that could invalidate the results.
A Dirichlet boundary condition, $u_{i}=0$, is applied for
all components of the flow velocity directly at the ground.
The surrounding walls are symmetry boundaries with frictionless properties.
This is achieved through a Neumann boundary condition $\partial u_{i} / \partial x_{i}=0$.

Tab. \ref{Tab:2A-1} gives the various speeds of this study with the respective Reynolds numbers
of the flow, using the runner's height for the definition of the characteristic length.
While all four speeds are considered for the assessment of single-runner aerodynamics (cf. Sec. \ref{Sec:3A}),
the drafting formations are studied at a running speed of $\SI{21}{\kilo\meter\per\hour}$ (cf. Sec. \ref{Sec:3B}-\ref{Sec:3D}).
This is a high speed occurring during elite marathon running and allows for good comparison between
the considered formations.
In fact, the current official marathon world record for men is 2:01:39 by Eliud Kipchoge
(September 16, 2018 during the Berlin Marathon; at an average speed of $\sim \SI{21}{\kilo\meter\per\hour}$)
and for women is 2:14:04 by Brigid Kosgei
(October 13, 2019 during the Chicago Marathon; at an average speed of $\sim \SI{19}{\kilo\meter\per\hour}$).
Here, we chose a running velocity at the upper range of these numbers to assess the ideal potential
of drafting with respect to metabolic quantities, drafted velocity, and time savings.

\begin{table}[htbp]
\caption{Running speeds and associated Reynolds numbers $Re=\left( U_{\infty} L \right)/\nu$
at a characteristic length of $L=\SI{1.65}{\meter}$.}
\centering
\begin{tabular}{@{}llll@{}}
\toprule
Running speed ($\SI{}{\kilo\meter\per\hour}$) & Reynolds number \\
\midrule
$15$ ($\SI{4.16}{\meter\per\second}$) & $\SI{452770}{}$ \\
$18$ ($\SI{5}{\meter\per\second}$) & $\SI{544195}{}$ \\
$21$ ($\SI{5.83}{\meter\per\second}$) & $\SI{634531}{}$ \\
$36$ ($\SI{10}{\meter\per\second}$) & $\SI{1088390}{}$ \\
\bottomrule
\end{tabular}
\label{Tab:2A-1}
\end{table}

The inherent challenge of simulations lies in ensuring that physical effects
are computed as accurately as possible.
Here, that means that flow separation lines and points, as well as boundary layer thickness and characteristics
are captured for the investigated Reynolds numbers, which lie in the turbulent regime (cf. Tab. \ref{Tab:2A-1}).
In order to avoid inaccurate prediction of aerodynamic forces on athletes at high Reynolds numbers,
as described by \citet{meile2006aerodynamics}, it is crucial to optimize near-surface grid resolutions,
such that no numerically induced jumps of separation points occur.
This is important, as an artificially delayed flow separation would lead to a sharp drop in the
pressure drag, which is the dominant component of total drag of bluff bodies (such as the human body).
A thorough verification and validation of the applied numerical approach is presented in \ref{App:A}.
The solution of the discretized continuity and momentum equations is performed
with the solver \emph{ANSYS Fluent 19.2}, post-processing of flow variables is achieved using
the post-processing tool \emph{ParaView} and drag distribution and
resulting power values are calculated with the free scientific programming
language \emph{GNU Octave}.

\subsection{\textbf{Geometry and mesh}}
\label{Sec:2B}

For the numerical setup we apply the realistic, three-dimensional  geometry of a female
runner with a height of $L=\SI{1.65}{\meter}$, a projected frontal area of $A\approx\SI{0.44}{\meter\squared}$,
and a body mass of $m\approx\SI{55}{\kilogram}$.
The geometry represents an average of 59 individual track athletes, scanned in-house and
positioned into a running pose. The obtained pose-invariant statistical model for the female human body
was reconstructed from the high-quality scan data using the method outlined by \citet{colaianni2014pose}.
The axial and lateral spacing between the runners is chosen according to estimates
of realistic running conditions.
\citet{polidori2020numerical} for instance uses an axial distance of $\SI{1.3}{\meter}$ and
a lateral shoulder-to-shoulder distance of $\SI{0.3}{\meter}$, which are taken from
an actual drafting situation during a marathon race.

The computational domain shown in Fig. \ref{Fig:2A-1}(a) has a distance of $\SI{7}{\meter}$
from the main runner to the upstream boundary and $\SI{13}{\meter}$ to the downstream boundary.
It contains far-field refinement as well as three refinement stages close to the runner's surfaces.
Furthermore, 10 prism layers (i.e. inflation layers) are inserted directly at all solid boundaries
with a steady growth factor of the prism layer thickness of $1.2$.
These measures ensure that the dimensionless wall distance of $y^{+}\approx1$ holds in the
entire domain and all boundary layers are well resolved.
Due to the local mesh refinement the cell count varies sharply between formations depending
on the number of runners included. It rises from $40$ million cells for the single-runner case
up to $107$ million cells for the final formation including a total of four runners.
The near-surface mean cell size of $\SI{3.125e-3}{\meter}$ and the outer-boundary mean cell size of $\SI{0.196}{\meter}$
for the final grid settings have been chosen following a grid convergence study (cf. \ref{Sec:A1}).
Homogeneous, gradual cell growth is applied between the near and far field.
Locally, especially in the boundary layer regions, perpendicular mesh dimensions are significantly lower (more than by a factor 10)
than the near-surface mean cell size, which is the cell size of the finest grid level just outside the prism layers.
A depiction of the considered formations of this study is shown in Fig. \ref{Fig:3B-1}-\ref{Fig:3B-2}.

\subsection{\textbf{Calculation of power values}}
\label{Sec:2C}

Various empirical models exist for the calculation of the mechanical power
output while running.
\citet{fukunaga1980effect} found through experiments on athletic runners
on force platforms a relation between running velocity $u$ and the
exerted power for forward motion as
\begin{equation}
P_{R,0} = 0.436u^{2.01},
\label{Eq:2C-1}
\end{equation}
with $u$ in $\SI{}{\meter\per\second}$ and $P_{R}$ in $\SI{}{\watt\per\kilogram}$.
In a study using motion tracking, \citet{cavagna1977mechanical} also
took into account additional energy expenditure created by the
movement of the limbs. They arrive at a function for the specific running power per unit mass as
\begin{equation}
P_{R,0} = 9.42 + 4.73u + 0.266u^{1.993}.
\label{Eq:2C-2}
\end{equation}
Here, $u$ is given in $\SI{}{\kilo\meter\per\hour}$ and $P_{R,0}$ in $\SI{}{\calorie\per\kilogram\per\minute}$.
For our considered cases, both approaches give similar results with
deviations in the range of only $2$-$5\%$, however,
Eq. (\ref{Eq:2C-2}) is used in the following power calculations due to
being the more complete model of running motion.
By `more complete' we mean that it takes into account both internal and
external work performed for the modelling of the total mechanical work
\begin{equation}
W_{tot} = W_{ext}+W_{int}.
\label{Eq:2C-2-1}
\end{equation}
Here, external work $W_{ext}$ is due to acceleration and lift of the centre of mass,
while internal work refers to both the translational and rotational acceleration of the
limbs relative to the trunk.
This approach has also been suggested more recently by the work of
\citet{saibene2003biomechanical}, \citet{pavei2019comprehensive}, and \citet{gray2020model}.

Using Eq. (\ref{Eq:2C-2}), we arrive at the power
\begin{equation}
P_{R} = P_{R,0} \cdot \frac{4.1868}{60} \cdot m
\label{Eq:2C-3}
\end{equation}
in $W$ for a runner of mass $m$ in $\SI{}{\kilogram}$, in the following simply called \emph{running power}.
Furthermore, there is a certain aerodynamic power $P_{A}$, used
to overcome drag $F_{D}$, which is defined as
\begin{equation}
P_{A} = F_{D} \cdot u = \frac{1}{2} \rho u^{3} C_{D} A.
\label{Eq:2C-4}
\end{equation}
Here, the expression for $F_{D}$ assumes stagnant air and the absence of crosswinds.
As a result, the total mechanical power generated while running
is the sum of the running power $P_{R}$ and the aerodynamic power $P_{A}$:
\begin{equation}
P_{tot} = P_{R}+P_{A}.
\label{Eq:2C-5}
\end{equation}
Using the total mechanical power $P_{tot}$ and the associated efficiency $\epsilon$
at a certain run speed by \citet{cavagna1977mechanical},
the metabolic power
\begin{equation}
P_{meta} = \frac{P_{tot}}{\epsilon},
\label{Eq:2C-6}
\end{equation}
which gives the metabolic energy expenditure during running, can be computed.

It must be emphasized here that the science behind efficiency of locomotion
and in particular of walking and running is far from concluded.
In both running and walking the muscles utilize energy stored during a previous phase of stretching
in the following phase of contraction.
As stated by \citet{cavagna1977mechanical}, the efficiency of walking reaches a maximum
of $0.35$-$0.40$ at intermediate speeds due to the expected properties of the contractile component of muscle.
This maximum of efficiency follows the force-velocity relation and the trends due to initial efficiency of muscle \citep{hill1964efficiency}.
Also in cycling the efficiency of locomotion was found to reach maximum values at intermediate speeds ($\sim 0.22$ according to the early study of \citet{dickinson1929efficiency},
$\sim 0.25$ according to \citet{ettema2009efficiency}).
In running, however, the efficiency increases steadily with speed (from $0.45$ to $0.70$), suggesting that
the positive work during that activity derives mainly from the passive recoil of elastic muscle tissue and
not from the active shortening of contractile muscle components \citep{cavagna1977mechanical, saunders2004factors}.
In support of this hypothesis, \citet{cavagna1968positive} found that the useful effect of pre-streched
muscle on the performed work increases with the speed of stretching and shortening.
Furthermore, \citet{komi1978utilization} investigated this utilization of elastic energy during
eccentric-concentric contraction of muscles in human motion and found
that up to $90\%$ of energy produced in the pre-stretching phase is recovered in case of
counter-movement jumps.
More recently, \citet{komi2000stretch} used data based on in-vivo force measurements, buckle-transducer
technique, and optic-fiber technique to show that the stretch-shortening cycle in human skeletal muscle
leads to considerable performance enhancement compared to simple concentric action.
With respect to running economy, \citet{hunter2015muscle} showed a clear positive correlation
of stretch-shortening cycle potentiation and increases of running economy
and recommended eccentric force development (e.g. by resistance training) for performance improvements.
Nonetheless, we point out that efficiency and total mechanical power models
are inherently difficult to define due to the sometimes conflicting understanding of internal and external
work performed during human locomotion and the question, whether internal work should appear in
total mechanical power and efficiency expressions, as has been suggested by
\citet{winter1979new} and more recently by \citet{minetti2001bipedalism}
(see e.g. \citet{ettema2009efficiency} for a discussion of these issues with regards to cycling,
and \citet{williams2000dynamics} with regards to running).

\section{\textbf{Results}}
\label{Sec:3}

\subsection{\textbf{Single-runner aerodynamics}}
\label{Sec:3A}

Using the geometry of the single runner,
the drag forces acting on various body parts are computed over the range of relevant velocities (cf. Fig. \ref{Fig:3A-1}(b)).
It can be seen that the drag acting on torso and legs is highest, with both parts roughly affected
by the same aerodynamic loads. Each make for a fraction of $39$-$40\%$ of the total drag on the runner.
Furthermore, arms are affected considerably less, with a share of approximately $13\%$,
and the head with $8\%$ of the total drag.
These estimates give an overview of the sensitivity of the various parts of the
body towards drag loads and are directly related to the frontal area of each part (cf. Eq. (\ref{Eq:2C-4})).
Fig. \ref{Fig:3A-1}(a) gives the total drag on the runner for the four simulated running speeds.
While the drag for lower speeds of $15$-$\SI{21}{\kilo\meter\per\hour}$, which are relevant for long distance
disciplines like marathon, lies in the range of $3.39$-$\SI{6.52}{\newton}$, it sharply rises to $\SI{18.45}{\newton}$
at the higher considered speed of $\SI{36}{\kilo\meter\per\hour}$. This demonstrates the (quadratically) increasing relevance of
drag reduction at higher speeds.
Moreover, the associated power values are calculated using the definitions of Sec. \ref{Sec:2C}.
Fig. \ref{Fig:3A-1}(c) shows the share of aerodynamic and running power in the total mechanical power
over the considered range of velocity.
However, the actual energy per unit time, which needs to be made available during motion, is the
metabolic power (cf. Eq. (\ref{Eq:2C-6})).

\begin{figure}[htbp]
\begin{centering}
\begin{overpic}[width=0.50\textwidth]{Figure2.pdf}
\end{overpic}
\par
\end{centering}
\caption{Total drag over the range of considered running speeds (a),
as well as drag acting on the body parts of the runner (b).
Total mechanical power, consisting of running power computed according to the empirical
model by \citet{cavagna1977mechanical} and aerodynamic power, as a function of speed (c).
Polynomial functions of second order (aerodynamic forces and running power) and of third order (aerodynamic and total mechanical power)
are used for fitting the data.}
\label{Fig:3A-1}
\end{figure}

Here, the aerodynamic component provides a fraction of the metabolic power of
$2.6\%$, $3.4\%$, $4.3\%$, and $8.5\%$ for
$\SI{4.16}{\meter\per\second}$ ($\SI{15}{\kilo\meter\per\hour}$), $\SI{5}{\meter\per\second}$ ($\SI{18}{\kilo\meter\per\hour}$), $\SI{5.83}{\meter\per\second}$
($\SI{21}{\kilo\meter\per\hour}$), and $\SI{10}{\meter\per\second}$ ($\SI{36}{\kilo\meter\per\hour}$).
These can be considered as threshold values for the maximum savings in metabolic energy consumption
during drafted running.

The resulting wake regions are depicted in Fig. \ref{Fig:3A-2}. They are visualized by
isosurfaces of total pressure, giving the boundary between the negative pressure region
in the wake and the positive pressure region due to recovery of the flow.
There is a notable decrease of the wake length from $\SI{1.14}{\meter}$ at $\SI{15}{\kilo\meter\per\hour}$ to
$\SI{0.98}{\meter}$ at $\SI{36}{\kilo\meter\per\hour}$ due to increased turbulent mixing.
The wake length is a good indication of the most beneficial range for drafting,
as the suction provided by the negative total pressure in the wake slightly pulls a following runner.
However, even at a length larger than that there will be a considerable drag-reducing effect due to flow
deceleration of the head wind.

\begin{figure*}[htbp]
\begin{centering}
\vspace{0.2cm}
\begin{overpic}[width=1.00\textwidth]{Figure3.pdf}
\end{overpic}
\par
\end{centering}
\caption{Wake region behind the runner visualized by isosurfaces of total pressure at $p_{tot}=0$.}
\label{Fig:3A-2}
\end{figure*}

\subsection{\textbf{Drag comparison of formations}}
\label{Sec:3B}

At first we consider the performance of the investigated formations with respect to
their effect on the main runner's drag at a running speed of $\SI{21}{\kilo\meter\per\hour}$.
It can be concluded from Fig. \ref{Fig:3B-1} that formation $2$ leads to the
largest reduction of drag of $75.6\%$, followed by formation $1$ with $70.1\%$.
Formations $3$ and $4$ cause significantly less drag reduction of $41.3\%$ and $33.4\%$, respectively.

\begin{figure}[htbp]
\begin{centering}
\begin{overpic}[width=0.50\textwidth]{Figure4.pdf}
\end{overpic}
\par
\end{centering}
\caption{Sketch of the investigated formations (a) together with their effect on the
main runner's drag (b). The formations are compared to the case of the single runner.
Axial and lateral spacing between the runners follows the values introduced in Fig. \ref{Fig:2A-1}(b).}
\label{Fig:3B-1}
\end{figure}

Furthermore, considering the drag distribution across the whole body of the main runner,
it can be recognized that the dominant part of the drag acts on legs and torso (cf. Sec. \ref{Sec:3A}).
Thus, it is crucial for a drafting formation to reduce the drag in these areas.
Looking at Fig. \ref{Fig:3B-2}, it can be seen that formations $1$-$2$ achieve this, while
formations $3$-$4$ cannot adequately shield the body of the main runner from oncoming air.
Additionally, it is shown that while the drag reduction on head and arms is fairly equal
for formations $1$-$3$, formation $4$ shows almost no improvement for the drag on the arms.
The reasons for the varying aerodynamic performances of formations $1$-$4$ are elaborated further in
Sec. \ref{Sec:3C}.

\begin{figure*}[htbp]
\begin{centering}
\begin{overpic}[width=1.00\textwidth]{Figure5.pdf}
\end{overpic}
\par
\end{centering}
\caption{Drag acting on the various body parts of the main runner alongside the total drag for
all involved runners, including pacers, for the formations 1 (a), 2 (b), 3 (c), and 4 (d).}
\label{Fig:3B-2}
\end{figure*}

Using the results for formation $1$, the distance between main runner and pacer was further varied
to give the drag values shown in Fig. \ref{Fig:3B-3}.
An exponential function of the type $y = a + b \cdot \exp \left( -c / x \right)$ with $a = 0$, $b = 5.63$, and $c = 1.28$
was used to fit the data and has the desired property of asymptotically approaching zero drag at small distances
and the single runner's maximum drag at large distances.
The values of the parameters were estimated by least-squares optimization using the Levenberg-Marquardt algorithm \citep{levenberg1944method,marquardt1963algorithm}.
The absolute axial limit of separation between the two runners, given the dimension and stride length of the applied geometry,
was found to be approximately $\SI{0.6}{\meter}$.
In reality, however, this limit will barely be reached and the runners will be unable to approach
already considerably earlier.
\citet{pugh1970oxygen} estimates a natural limit of axial spacing of $\SI{1}{\meter}$ for middle-distance running
and states a possible drag reduction at this distance of up to $80\%$.
In this study we find a drag of $\SI{1.45}{\newton}$ acting on the main runner at a distance of $\SI{1}{\meter}$, which corresponds to $22.2\%$
of the drag on the single runner and thereby lies in good agreement with the results by \citet{pugh1970oxygen}.

\begin{figure}[htbp]
\begin{centering}
\begin{overpic}[width=0.50\textwidth]{Figure6.pdf}
\end{overpic}
\par
\end{centering}
\caption{Drag on the main runner as a function of the distance to the leading pacer in formation 1.
The numerical data is fitted using an exponential function with a determination coefficient of $R^{2}=0.99$.
In reality, there is an axial limit which prohibits the further approach of the main runner to the pacer.
This limit is estimated from geometric considerations as $\SI{0.6}{\meter}$. \citet{pugh1970oxygen} estimates
this minimum distance during drafting at $\SI{1.0}{\meter}$.}
\label{Fig:3B-3}
\end{figure}

\subsection{\textbf{Aerodynamic effects on runners within formations}}
\label{Sec:3C}

Air resistance while running manifests itself in static pressure increase
as the air is being decelerated and diverted at the front of the runner.
Fig. \ref{Fig:3C-1} provides insight to the magnitude of the pressure acting
on the main runner for the various formations. The pressure coefficient
\begin{equation}
c_{p} = \frac{p-p_{\infty}}{\frac{1}{2} \rho U_{\infty}^{2}} = \frac{p-p_{\infty}}{p_{0}-p_{\infty}},
\label{Eq:3B-1}
\end{equation}
with the free-stream pressure of $p_{\infty}= \SI{101325}{\pascal} = \SI{1}{\atm}$
and a relative free-stream velocity of $U_{\infty} = \SI{5.83}{\meter\per\second} = \SI{21}{\kilo\meter\per\hour}$,
indicates regions of high or low relative pressure.

\begin{figure*}[htbp]
\begin{centering}
\begin{overpic}[width=1.00\textwidth]{Figure7.pdf}
\end{overpic}
\par
\end{centering}
\caption{Pressure coefficient on the main runner for the single-runner case (a) compared to
the formations 1-4 (b)-(e).}
\label{Fig:3C-1}
\end{figure*}

Fig. \ref{Fig:3C-1} confirms the beneficial pressure distribution for formations $1$ and $2$.
This means that there are no large areas of high pressure coefficient at the front of the runner,
nor any large regions of low negative pressure coefficient at the back, which
would be an indication for suction force slightly pulling the runner back.
Formation $3$ on the other hand shows a high positive pressure coefficient at the
front, which already hints at the leakage of oncoming flow between the front pacers
that can be seen in Fig. \ref{Fig:3C-2}(c) and Fig. \ref{Fig:3C-3}(c).
While formation $4$ does not display any significantly high pressure peaks at the
front, it is the elevated negative pressure coefficient at the back that makes it less
advantageous than formations $1$-$2$ (cf. Sec. \ref{Sec:3B}).

\begin{figure*}[htbp]
\begin{centering}
\vspace{0.4cm}
\begin{overpic}[width=1.00\textwidth]{Figure8.pdf}
\end{overpic}
\par
\end{centering}
\caption{Streamlines seeded at the front of the leading pacers showing the diversion
of air around the main runner for all four formations.
Streamlines are coloured by velocity magnitude and the runners' surfaces by pressure coefficient.}
\label{Fig:3C-2}
\end{figure*}

\begin{figure*}[htbp]
\begin{centering}
\begin{overpic}[width=1.00\textwidth]{Figure9.pdf}
\end{overpic}
\par
\end{centering}
\caption{Velocity magnitude contours in the horizontal plane section through the runners' centre.
Formations 1-4 are shown in images (a)-(d).}
\label{Fig:3C-3}
\end{figure*}

Fig. \ref{Fig:3C-2}-\ref{Fig:3C-3} both demonstrate the reasons for the large differences
in aerodynamic performance between the formations.
In case of formations $1$ and $2$ the main runner is well embedded inside the wake region of
the front pacer and therefore only affected by low-speed, almost stagnant air relative to herself.
Formation $2$ has the additional advantage of a rear pacer, whose high-pressure stagnation
region at the front positively affects the low-pressure wake region of the main runner by
raising the pressure coefficient to a value closer to zero, thus lowering suction.
This can be seen when comparing Fig. \ref{Fig:3C-4}(a) and Fig. \ref{Fig:3C-4}(b).
Furthermore, Fig. \ref{Fig:3C-2}(c) and Fig. \ref{Fig:3C-3}(c) clearly demonstrate a major weakness
of formation $3$, which is the passing air between the front pacers.
This leads to a large region of high stagnation pressure directly at the area of the impacting air
on the main runner, which is also visible in Fig. \ref{Fig:3C-1}(d) and  Fig. \ref{Fig:3C-4}(c).
Also formation $4$ is suboptimal due to air being accelerated between the main runner and the
pacers on the side. This leads to the aforementioned effect of higher drag at the arms
compared to the other formations (cf. Sec. \ref{Sec:3B}).
Additionally, this effect of locally increased flow velocity between the runners leads to a shortening
of the wake region of slow air for the main runner and thus negatively impacts her drafting conditions (cf. Fig. \ref{Fig:3C-3}(d)).

\begin{figure*}[htbp]
\begin{centering}
\begin{overpic}[width=0.75\textwidth]{Figure10.pdf}
\end{overpic}
\par
\end{centering}
\caption{Pressure coefficient in the vertical plane section through the main runner's centre.
Formations 1-4 are depicted in images (a)-(d).}
\label{Fig:3C-4}
\end{figure*}

\subsection{\textbf{Metabolic power savings and performance predictions}}
\label{Sec:3D}

Using the drag values from Sec. \ref{Sec:3B}, we can compute the associated aerodynamic
power and the total mechanical power using Eq. (\ref{Eq:2C-2})-(\ref{Eq:2C-5}).
Tab. \ref{Tab:3D-1} lists the results for the single runner and the four formations.
Additionally, relative changes in aerodynamic power and total power are given by
comparison with the results of the single-runner case.
The resulting metabolic power needed for running in the considered formations is computed by
using the efficiency of running at the given velocity as detailed in \citet{cavagna1977mechanical}.
Furthermore, running economy as the specific metabolic rate per unit mass is given together with its relative changes in Tab. \ref{Tab:3D-2}.
As described by \citet{kipp2019extrapolating}, the relationship between improvements in running economy
and improvements in velocity is not a linear, but rather a curvilinear one.
This means that for a given improvement in running economy, velocity gains will be smaller at higher speeds and larger at lower speeds.
At elite marathon speeds of around $~\SI{5.5}{\meter\per\second}$ there is roughly a $2/3^{\text{rds}}$ increase in velocity for each percent increase in running economy.
By applying the conversion factors of \citet{kipp2019extrapolating} we arrive at estimates for the velocity gains for each formation in Tab. \ref{Tab:3D-3}. The relative changes in speed are $2.1\%$ (formation 1), $2.3\%$ (formation 2), $1.3\%$ (formation 3), and $1.1\%$ (formation 4).
\begin{table*}[htbp]
\caption{Aerodynamic power $P_{A}$, and total mechanical power $P_{tot}$
computed with Eq. (\ref{Eq:2C-4}) and (\ref{Eq:2C-5}) and by computing the
running power $P_{R}=\SI{858.1}{\watt}$ at $\SI{21}{\kilo\meter\per\hour}$ ($\SI{5.83}{\meter\per\second}$) using the model of \citet{cavagna1977mechanical}.
Savings $\Delta P_{A}$ and $\Delta P_{tot}$ are given with respect to the case of the runner alone.}
\centering
\begin{tabular}{@{}llllll@{}}
\toprule
Formation & $F_{D}$ ($\SI{}{\newton}$) & $P_{A}$ ($\SI{}{\watt}$) & $\Delta P_{A}$ & $P_{tot}$ ($\SI{}{\watt}$) & $\Delta P_{tot}$ \\
\midrule
Single runner & $6.52$ & $38.0$ & - & $896.1$ & - \\
$1$ & $1.96$ & $11.4$ & $70.0\%$ & $869.5$ & $3.0\%$ \\
$2$ & $1.58$ & $9.2$ & $75.8\%$ & $867.3$ & $3.2\%$ \\
$3$ & $3.84$ & $22.4$ & $41.1\%$ & $880.5$ & $1.7\%$ \\
$4$ & $4.33$ & $25.3$ & $33.4\%$ & $883.4$ & $1.4\%$\\
\bottomrule
\end{tabular}
\label{Tab:3D-1}
\end{table*}
\begin{table*}[htbp]
\caption{Metabolic power $P_{meta}$ according to Eq. (\ref{Eq:2C-6}) and running economy $RE$, as well as percent improvements of the
running economy $\Delta RE$ for each of the formations with respect to the single-runner case.
Efficiency of locomotion for a running velocity of $\SI{21}{\kilo\meter\per\hour}$ ($\SI{5.83}{\meter\per\second}$) is taken as $\epsilon=0.64$ \citep{cavagna1977mechanical}.}
\centering
\begin{tabular}{@{}llll@{}}
\toprule
Formation & $P_{meta}$ ($\SI{}{\watt}$) & $RE$ ($\SI{}{\watt\per\kilogram}$) & $\Delta RE$ \\
\midrule
Single runner & $1400.2$ & $25.5$ & - \\
1 & $1358.6$ & $24.7$ & $3.1\%$ \\
2 & $1355.2$ & $24.6$ & $3.5\%$ \\
3 & $1375.8$ & $25.0$ & $2.0\%$ \\
4 & $1380.3$ & $25.1$ & $1.6\%$ \\
\bottomrule
\end{tabular}
\label{Tab:3D-2}
\end{table*}
\begin{table*}[htbp]
\caption{Drafted running speeds $u_{draft}$ for all formations alongside potential time $t_{draft}$
and time savings through drafting $\Delta t_{draft}$ over the marathon distance of $\SI{42195}{\meter}$ by comparison with the undrafted single-runner case.}
\centering
\begin{tabular}{@{}lllll@{}}
\toprule
Formation & $u_{draft}$ ($\SI{}{\meter\per\second}$) & $t_{draft}$ ($\SI{}{\second}$) & $\Delta t_{draft}$ ($\SI{}{\second}$) \\
\midrule
1 & $5.95$ & $7092$ ($\approx\SI{1.970}{\hour}$) & $142$ ($\approx\SI{2.4}{\minute}$) \\
2 & $5.96$ & $7080$ ($\approx\SI{1.967}{\hour}$) & $154$ ($\approx\SI{2.6}{\minute}$) \\
3 & $5.91$ & $7140$ ($\approx\SI{1.983}{\hour}$) & $98$ ($\approx\SI{1.6}{\minute}$) \\
4 & $5.89$ & $7164$ ($\approx\SI{1.990}{\hour}$) & $70$ ($\approx\SI{1.2}{\minute}$) \\
\bottomrule
\end{tabular}
\label{Tab:3D-3}
\end{table*}

\section{Discussion and conclusions}

Drafting and the resulting aerodynamic drag reduction cause a notable effect on running
and the energy expenditure used for it.
This study shows the overall aerodynamics and associated energy expenditure during running
for a range of relevant speeds. Considerable shares of up to $8.5\%$ at a velocity of $\SI{36}{\kilo\meter\per\hour}$ are computed.

Additionally, differences in drag, pressure distribution across the runners,
and potential power savings are demonstrated for four formations that could arise during running competitions.
As shown in Sec. \ref{Sec:3B}, simple formations of axial spacing of the pacers
(i.e. one pacer in front, as in formation 1, or one pacer in front and one in the back, as in formation 2) already lead to considerable
reduction of over $70\%$ in drag at a distance of $\SI{1.2}{\meter}$.
Furthermore, formations 1 and 2 also show a clear positive effect on leading and trailing
pacers, since their overpressure and underpressure regions interact with the main runner.
This effect has also been noted by \citet{beaumont2019does} together with a measurable
advantage for the pacers with respect to oxygen consumption.

We further deduce an empirical relationship of drag force acting on the main runner versus
her distance to a leading pacer, which can be used for future prediction of a pacer's aerodynamic effects (cf. Fig. \ref{Fig:3B-3}).
The results are congruent with experimental data by \citet{pugh1970oxygen},
who predicts drag decrease of up to $80\%$ when running behind a pacer at $\SI{1}{\meter}$.
\citet{davies1980effects}, who claimed possible drag savings in the range of $80$-$85\%$,
also confirms these results.

Moreover, the drag distribution at various body parts is examined in Sec. \ref{Sec:3B}.
From this and the explanations in Sec. \ref{Sec:3C}, we can see a large sensitivity of the
drag acting on the main runner if pacers are positioned laterally, such as in formations $3$ and $4$.
In formation $3$, flow leaks between the front pacers, causing large pressure peaks at the runner's
centre (cf. Fig. \ref{Fig:3C-1}).
Formation $4$ on the other hand has the negative effect of pacers at the runner's
sides and the resulting acceleration of air between the runners leads to a relatively low
drag reduction of only $33.4\%$ compared to formation $1$ with $70.1\%$.
This suggests that runners at the side should be avoided if optimal drag savings are the goal.

In essence, there are three aerodynamic effects of drafting that lead to a reduction of
axial force (i.e. drag) while running:
\begin{enumerate}[(i)]
\item \emph{Shielding against fast, oncoming air.} \\
Therefore, there is no sharp deceleration of air and resulting formation of areas of high stagnation
pressure on the runner's surface.
\item \emph{Suction effects due to the region of negative total pressure in the leading runner's wake.} \\
Any object moving through air forms a wake region of relatively low pressure.
This underpressure can be harnessed to the benefit of a trailing runner.
The wake region of negative total pressure has been found in this study to extend to about $\SI{1.14}{\meter}$ at $\SI{15}{\kilo\meter\per\hour}$ down to $\SI{0.98}{\meter}$ at
$\SI{36}{\kilo\meter\per\hour}$ from the trailing runner's back.
\item \emph{Delay of separation points due to turbulent free shear layers.} \\
Turbulence levels are naturally increased in the high-shear-stress regions of the wake,
where air is displaced against the stagnant surroundings.
This raises momentum transfer in the boundary layers of the runner and delays separation due
to formation of a turbulent boundary layer from a laminar one. As a result, pressure drag (i.e. form drag)
drops, which is the dominant source of drag of bluff bodies.
\end{enumerate}

Eventually, the aerodynamic power savings are computed in this study and listed together with the
total mechanical power and the metabolic power output in Sec. \ref{Sec:3D}.
The derived metabolic power and running economy of up to $\SI{1400.2}{\watt}$ and $\SI{25.5}{\watt\per\kilogram}$
for the single runner at a speed of $\SI{21}{\kilo\meter\per\hour}$ agrees well with values from literature for trained athletes \citep{kipp2017does, batliner2018does}.
The largest improvements of running economy occur for formation $1$ with $3.1\%$
and formation $2$ with $3.5\%$.
These metabolic gains are slightly higher than the $2.84\%$ gain computed
by \citet{polidori2020numerical} for the best position of Kenenisa Bekele's drafting strategy
during the 2019 Berlin marathon.
The reason for this is most likely the larger axial spacing of $\SI{1.3}{\meter}$ to the front pacer
as opposed to $\SI{1.2}{\meter}$ in the current study.
Furthermore, a particular benefit of this study is the additional performance prediction in
terms of improvements of velocity and time for the investigated formations.
By applying the curvilinear relationship between gains in velocity and running economy suggested by \citet{kipp2019extrapolating},
we deduce possible speed gains in the range of $1.1$-$2.3\%$.
Time savings of $\SI{1.2}{\minute}$ (formation 4) to $\SI{2.6}{\minute}$ (formation 2) are achieved over the marathon distance.

In this study we focus on the impact of aerodynamic drafting on
running economy.
However, there are a multitude of other factors that influence running economy and performance.
With respect to biomechanics, some of these factors are (lower) limb mass distribution, Achilles tendon moment arm
and lower body musculotendinous structures for reutilization of elastic energy,
as well as running style and gait patterns (e.g. so-called \emph{pose running}
to benefit from elastic recoil forces in the lower limbs), among others \citep{williams1987relationship, saunders2004factors, barnes2015running}.

Concluding from this work, it is possible to give targeted suggestions for formation
patterns that are best suited for cooperative drafting.
It is shown that axially positioned pacers give the best performance and that laterally positioned pacers
should be running at a distance further away from the main runner to avoid local flow acceleration, while
still contributing positively to drag reduction.
Such an approach has been applied with the recent successful world record attempt by Eliud Kipchoge
during the \emph{INEOS 1:59 Challenge}, where a staggered formation in the shape of an inverted V has been chosen.
Future studies could explore a larger number of formations by CFD and subsequently
derive reduced-order models for parameters such as number of pacers, or axial and lateral spacing.
\bigskip

\noindent
\textbf{Acknowledgements} \\
This work was realised through research budget from ADIDAS and KTM E-TECHNOLOGIES.
Computational resources were provided by the parallel computing cluster at KTM in Mattighofen, Austria.
\bigskip

\noindent
\textbf{Conflict of interest statement} \\
The authors have nothing to disclose that would have biased this work.


\section{Introduction}
In this paper, we continue studies of critical phenomena, namely shock waves and phase transitions, appearing in solutions to Euler equations describing flows of gases \cite{LR-ljm-shock,LR-sym}. Such effects have always been paid a great deal of attention both because of their mathematical beauty \cite{Arn1,Arn2,Arn3} and practical applications \cite{Zeld}. The tendency of studying such phenomena is remaining nowadays as well, see, for example \cite{Huang}, where the case of Chaplygin gases is considered, \cite{Ros-Tab,Chat}, where the authors discuss weak shock waves, which is actually the case considered in the present paper, and also it is worth mentioning \cite{Polud1,Polud2}, where the influence of turbulence on detonations is investigated. The properties of global solvability for Euler equations and singularities of their solutions were also studied in a series of works \cite{Tun1,Tun2,Tun3,Tun-Bog}.

Our approach to studying and finding singular properties of solutions to nonlinear PDEs is essentially based on a geometrical theory of PDEs \cite{KrVin,KVL,Ovs,Olver}. Namely, it is known that Euler equations considered here are a particular case of Jacobi type systems (see, for example,\cite{KLR}), which have a natural representation in terms of differential 2-forms on 0-jet space. This observation goes back to a seminal paper \cite{Lych}. One of advantages of this approach is that we need to deal with geometrical structures on 0-jet space instead of 1-jet space, where the equations in question have a natural representation. This idea has also found applications in incompressible hydrodynamics \cite{RoulRub1,RoulRub2}.

Also this approach allows to extend the notion of a smooth solution to a generalized solution understood as an integral manifold of the mentioned forms, which makes it possible to find solutions that are not globally given by functions. Then, singularities of projections of multivalued solutions to the space of independent variables is exactly what corresponds to a formation of shock waves \cite{LychSing}. This concept has been used to describe shock waves in non-stationary filtration problems \cite{AKL-dan,AKL-ifac,AKL-gsa}.

 We combine this observation with another approach to finding multivalued solutions to nonlinear PDEs, which is adding a differential constraint compatible with the original system \cite{KrugLych} that appeared to be fruitful, in particular, in applications to 2-dimensional gas dynamics \cite{Lych-Yum-1}, the Khokhlov-Zabolotskaya equation \cite{Lych-Yum-2}, and also the Hunter-Saxton equation \cite{Eiv}. In the present paper, compared with \cite{LR-ljm-shock,LR-sym}, we give a more accurate description of finding such constraints and discuss two methods of finding them. One of them is based on the mentioned specific geometry that Euler equations have, and another one is based on the theory of differential invariants and quotient PDEs \cite{KLych,Eiv,DLT-sym}. It is worth saying that the last one is more general, however, we decided to discuss both of them to give a more complete picture of the geometry underlying non-stationary Euler equations in one spatial dimension.

The paper is organized as follows. First, we briefly discuss thermodynamics in terms of contact and symplectic geometries \cite{Gibbs,Mrug,Rup,Lych-therm,DLT1,LR-ljm-1}, since it is significant in description of flows of gases. Then, we turn to Euler equations and describe two methods of finding multivalued solutions, using their specific geometric description and using the concept of a quotient PDE. Finally, we get a new class of exact multivalued solutions for an arbitrary thermodynamic model and illustrate them on flows of ideal gases. We find caustics and shock wave fronts. It is worth saying that one can elaborate these solutions for various thermodynamic models, for instance, the van der Waals model taking into account phase transitions, which will be addressed in future papers.

The system of PDEs describing the motion of gases consists of three conservation laws \cite{Batch,Lan-Lif}.

\begin{itemize}
\item Conservation of momentum

\begin{equation}
\label{Euler}
\rho(u_{t}+uu_{x})=-p_{x},
\end{equation}

\item Conservation of mass

\begin{equation}
\label{mass}
\rho_{t}+(\rho u)_{x}=0,
\end{equation}

\item Conservation of entropy along the flow

\begin{equation}
\label{energy}
s_{t}+us_{x}=0,
\end{equation}

\end{itemize}
where $u(t,x)$ is the velocity, $\rho(t,x)$ is the density, $p(t,x)$ is the pressure, $s(t,x)$ is the specific entropy.

One can see that system (\ref{Euler})-(\ref{energy}) is incomplete, which is natural, since the very medium we describe has not yet been specified. This is usually done either by incompressibility condition, when the density is put constant, or by \textit{state equations}, representing the relations between various thermodynamic quantities. Here, we concentrate on the last approach.

\section{Thermodynamics}
We start with the discussion of equilibrium thermodynamics of gases. Any thermodynamical system in equilibrium is described by the following thermodynamical quantities: $e$ --- specific inner energy, $v=\rho^{-1}$ --- specific volume (or inverse of density), $T$ --- temperature, $p$ --- pressure, and $s$ --- specific entropy. The main law of thermodynamics, which is the energy conservation law, states that the differential 1-form
\begin{equation}
\label{theta}
\theta=-ds+T^{-1}de-pT^{-1}\rho^{-2}d\rho
\end{equation}
must vanish. This drives us to the notion of \textit{thermodynamic states} understood as maximal integral manifolds of \eqref{theta}.

More precisely, consider the contact manifold $\widehat{\Phi}=W\times W^{\ast}\times\mathbb{R}\simeq\mathbb{R}^{5}$, where $W\simeq\mathbb{R}^{2}(e,\rho)$, $W^{\ast}\simeq\mathbb{R}^{2}(p,T)$, and $s\in\mathbb{R}$. Then, a thermodynamic state is a Legendrian manifold $\widehat{L}\subset(\widehat{\Phi},\theta)$, such that
\begin{equation*}
\theta|_{\widehat{L}}=0.
\end{equation*}
This exactly means that the energy conservation law holds on $\widehat{L}$.

Let us choose $(e,\rho)$ as local coordinates on $\widehat{L}$. Then, the two-dimensional manifold $\widehat{L}\subset(\widehat{\Phi},\theta)$ is given by three relations
\begin{equation}
\label{legmani}
\widehat{L}=\left\{s=S(e,\rho),\, T=\frac{1}{S_{e}},\, p=-\rho^{2}\frac{S_{\rho}}{S_{e}}\right\}
\end{equation}
for some function $S(e,\rho)$.

However, in practice there are no ways of determining $S(e,\rho)$, which motivates us to switch to the Lagrangian viewpoint. Namely, consider projection
\begin{equation*}
\pi\colon\widehat{\Phi}\to\Phi,\quad\pi(p,T,e,\rho,s)=(p,T,e,\rho).
\end{equation*}
A pair $(\Phi,\Omega)$ is a symplectic manifold with the structure form
\begin{equation*}
\Omega=d\theta=d(T^{-1})\wedge de-d(pT^{-1}\rho^{-2})\wedge d\rho.
\end{equation*}

Then, a thermodynamic state is a Lagrangian manifold $\pi(\widehat{L})=L\subset(\Phi,\Omega)$, such that
\begin{equation*}
\Omega|_{L}=0.
\end{equation*}
In the symplectic space $(\Phi,\Omega)$ the Lagrangian manifold $L$ is given by \textit{state equations}:
\begin{equation*}L=\left\{f(p,T,e,\rho)=0,\,g(p,T,e,\rho)=0\right\},\end{equation*}
such that the Poisson bracket $[f,g]$ with respect to the structure form $\Omega$
\begin{equation*}[f,g]\,\Omega\wedge\Omega=df\wedge dg\wedge\Omega\end{equation*}
vanishes on $L$:
\begin{equation}
\label{compat}
[f,g]=0\text{ on }L.
\end{equation}
Condition $[f,g]=0$ on $L$ is called the \textit{compatibility condition} for state equations.

If one chooses $(T,\rho)$ as local coordinates on the Lagrangian manifold $L$, that is
\begin{equation*}
L=\left\{p=P(T,\rho),\,e=E(T,\rho)\right\}.
\end{equation*}
The condition $[f,g]=0$ on $L$ leads to the equation
\begin{equation*}(-\rho^{-2}T^{-1}P)_{T}=(T^{-2}E)_{\rho},
\end{equation*}
and therefore the following theorem is valid:
\begin{theorem}
The Lagrangian manifold $L$ is given by means of \textit{the Massieu-Planck potential} $\phi(\rho,T)$
\begin{equation}
\label{lag-mani}
p=-\rho^{2}T\phi_{\rho},\quad e=T^{2}\phi_{T}.
\end{equation}
\end{theorem}
\begin{remark}
One can build up a Legendrian manifold \eqref{legmani} from a given Lagrangian one by resolving the overdetermined system of equations on function $S(e,\rho)$:
\begin{equation}
\label{sys-ent}
T=\frac{1}{S_{e}},\quad p=-\rho^{2}\frac{S_{\rho}}{S_{e}},
\end{equation}
where $p(e,\rho)$ and $T(e,\rho)$ are specified once \eqref{lag-mani} is given, and the compatibility condition for \eqref{sys-ent} is \eqref{compat}. Indeed, the specific entropy is expressed in terms of the Massieu-Planck potential as follows \cite{LR-ljm-1}:
\begin{equation*}
s=\phi+T\phi_{T}.
\end{equation*}
\end{remark}

The manifold $\Phi=W\times W^{\ast}$ is equipped with the differential quadratic form \cite{Lych-therm}
\begin{equation*}
\kappa=d(T^{-1})\cdot de-\rho^{-2}d(pT^{-1})\cdot d\rho.
\end{equation*}
Not all the points on the Lagrangian manifold correspond to real physical states, but only those where (applicability condition)
\begin{equation*}
\kappa|_{L}<0.
\end{equation*}
Writing down $\kappa|_{L}$ in terms of the Massieu-Planck potential, we get:
\begin{equation*}
\kappa|_{L}=-(2T^{-1}\phi_{T}+\phi_{TT})dT\cdot dT+(2\rho^{-1}\phi_{\rho}+\phi_{\rho\rho})d\rho\cdot d\rho.
\end{equation*}
Taking into account \eqref{lag-mani} we observe that the applicability condition $\kappa|_{L}<0$ is
\begin{equation*}
e_{T}>0,\quad p_{\rho}>0,
\end{equation*}
which is known to be the conditions of thermodynamic stability.

In the context of the geometrical approach to thermodynamics, \textit{thermodynamic processes} are understood as contact transformations of $\widehat{\Phi}$, preserving the Legendrian manifold $\widehat{L}$. They are generated by contact vector fields, and their integral curves $l\subset\widehat{L}$ we will also call \textit{thermodynamic processes}.
\section{Euler equations}
Let us note that system (\ref{Euler})-(\ref{energy}) becomes complete once we fix a thermodynamic state \eqref{legmani}.

We will apply a \textit{homentropicity ansatz}, which means that we put $s(t,x)=s_{0}=\mathrm{const}$. On the one hand, we get rid of equation \eqref{energy}, on the other hand, we are able to express all the thermodynamic quantities in terms of the density $\rho$. Indeed, considering $s_{0}=\phi+T\phi_{T}$ as an equation for $T(\rho)$, we observe that the derivative of its right-hand side with respect to $T$ is positive in an applicable domain, and therefore due to the implicit function theorem, this equation determines $T(\rho)$ uniquely. By means of state equations we also get $p=p(\rho)$. Summarizing above discussion, we end up with the following two-component system of PDEs:

\begin{equation}
\label{EU}
\mathcal{E}=\begin{cases}
 \rho_{t}+(\rho u)_{x}=0,
   \\
\displaystyle u_{t}+uu_{x}+\frac{p^{\prime}(\rho)}{\rho}\rho_{x}=0.
 \end{cases}
\end{equation}
\subsection{Geometrical structures associated to Euler equations}
Let $E=J^{0}(t,x,u,\rho)$ be the space of 0-jets, and $M=\mathbb{R}^{2}(t,x)$. Then, following \cite{KLR,Lych}, one can associate two differential 2-forms
\begin{equation*}
\begin{cases}
 \omega_{1}=\rho dt\wedge du+udt\wedge d\rho-dx\wedge d\rho,
   \\
\displaystyle \omega_{2}=udt\wedge du+\frac{p^{\prime}(\rho)}{\rho}dt\wedge d\rho-dx\wedge du
 \end{cases}
\end{equation*}
with system \eqref{EU}. Indeed, let $\Omega^{2}(E)$ be a module of differential 2-forms on $E$, then any differential 2-form $\omega\in\Omega^{2}(E)$ generates an operator
\begin{equation*}
\Delta_{\omega}\colon C^{\infty}(M)\to\Omega^{2}(M),\quad \Delta_{\omega}(f)=\omega|_{\Gamma^{0}(f)},
\end{equation*}
where $\Gamma^{0}(f)\subset E$ is a graph of the vector-function $f$. The system $\mathcal{E}$ can now be written as
\begin{equation*}
\Delta_{\omega_{1}}(f)=0,\quad\Delta_{\omega_{2}}(f)=0,
\end{equation*}
where $f=(u(t,x),\rho(t,x))$.

A 2-dimensional manifold $N\subset E$ is said to be a \textit{(multivalued) solution} of $\mathcal{E}$, if $\omega_{1}|_{N}=\omega_{2}|_{N}=0$.

Note that the form $\omega_{2}$ is closed and non-degenerate, and therefore may serve as a symplectic structure on $E$, which means that any multivalued solution to $\mathcal{E}$ is a Lagrangian manifold.

Let us fix a volume form on $E$ as $q=dt\wedge dx\wedge du\wedge d\rho$ and introduce a bilinear operator
\begin{equation*}
P\colon\Omega^{2}(E)\times\Omega^{2}(E)\to C^{\infty}(E),\quad \alpha_{1}\wedge\alpha_{2}=P(\alpha_{1},\alpha_{2})q,\quad\alpha_{i}\in\Omega^{2}(E).
\end{equation*}
Let $P_{\omega}=\|P(\omega_{i},\omega_{j})\|$ be its matrix in the basis $\langle\omega_{1},\omega_{2}\rangle$. Then, system $\mathcal{E}$ is said to be \textit{hyperbolic}, if $\det(P_{\omega})<0$, \textit{elliptic}, if $\det(P_{\omega})>0$ and \textit{parabolic}, if $\det(P_{\omega})=0$. In our case the matrix $P_{\omega}$ has the following form:
\begin{equation*}
P_{\omega}=
\begin{pmatrix}
2\rho & 0\\
0 & -2\rho^{-1}p^{\prime}(\rho)
\end{pmatrix},
\end{equation*}
and therefore the condition for the system $\mathcal{E}$ to be of hyperbolic type is $p^{\prime}(\rho)>0$.

It is worth saying that the applicability condition $e_{T}>0$ is satisfied for a considerable number of real gas models on the entire Lagrangian manifold, and we will consider only such models, while another applicability condition $p_{\rho}>0$ violates in some domains of Lagrangian manifolds of real gases, for instance, van der Waals model \cite{LR-sym}. Thus the applicability condition $p_{\rho}>0$ is essential for us both in the context of thermodynamics and in the context of type of the Euler system.
\begin{theorem}
If a thermodynamic process curve $l\subset\widehat{L}$ lies in the domain of negativity of the form $\kappa|_{L}$, then the system $\mathcal{E}$ is hyperbolic.
\end{theorem}

Obviously, the forms $\omega_{1}$ and $\omega_{2}$ are defined up to a non-degenerate linear transformation. Indeed, the forms
\begin{eqnarray*}
\label{forms}
\widehat{\omega_{1}}&=&a_{11}\omega_{1}+a_{12}\omega_{2},\\
\widehat{\omega_{2}}&=&a_{21}\omega_{1}+a_{22}\omega_{2},
\end{eqnarray*}
where $a_{11}a_{22}-a_{12}a_{21}\ne 0$, define the Euler system as well. In a hyperbolic case one can choose $a_{ij}$ in such a way that
\begin{equation}
\label{eff}
\omega_{1}\wedge\omega_{2}=0,\quad \omega_{1}\wedge\omega_{1}=-\omega_{2}\wedge\omega_{2}.
\end{equation}
Differential 2-forms satisfying relations \eqref{eff} are also called \textit{effective} \cite{KLR}. For the case of Euler system, we get the following theorem \cite{LR-ljm-shock}:
\begin{theorem}
Let $\mathcal{E}$ be a system of hyperbolic type. Then, it can be given by 2-forms
\begin{eqnarray*}
\omega_{1}&=&A(\rho)(\rho dt\wedge du+udt\wedge d\rho-dx\wedge d\rho),\\
\omega_{2}&=&udt\wedge du+\rho A^{2}(\rho)dt\wedge d\rho-dx\wedge du,
\end{eqnarray*}
where $A(\rho)=\rho^{-1}\sqrt{p^{\prime}(\rho)}$, and 2-forms $\omega_{1}$, $\omega_{2}$ satisfy relations (\ref{eff}).
\end{theorem}

Let us introduce another operator
\begin{equation*}
A_{\omega}\colon D(E)\to D(E),\quad X\rfloor\omega_{2}=A_{\omega}(X)\rfloor\omega_{1},
\end{equation*}
where $D(E)$ is the module of vector fields on $E$, and $X\in D(E)$.

If one chooses $\langle\partial_{t},\partial_{x},\partial_{\rho},\partial_{u}\rangle$ as a basis in $D(E)$, one gets the matrix of this operator
\begin{equation*}
W=\frac{1}{\rho A(\rho)}\begin{pmatrix}
u & -1 & 0 & 0\\
u^{2}-\rho^{2}A^{2}(\rho) & -u & 0 & 0\\
0 &0 & 0 & \rho A^{2}(\rho)\\
0 &0 & \rho & 0
\end{pmatrix},
\end{equation*}
and moreover, in a hyperbolic case $A_{\omega}^{2}=\mathrm{id}$.

Let $a\in E$, $\mathcal{C}_{\pm}(a)$ be eigenspaces of the operator $A_{\omega}(a)$, then
\begin{equation*}
T_{a}E=\mathcal{C}_{-}(a)\oplus\mathcal{C}_{+}(a),
\end{equation*}
and \textit{characteristic distributions} $\mathcal{C}_{+}=\langle X_{+},Y_{+}\rangle$ and $\mathcal{C}_{-}=\langle X_{-},Y_{-}\rangle$ are generated by vector fields \cite{LR-ljm-shock}
\begin{eqnarray*}
X_{\pm}&=&\pm A(\rho)\partial_{u}+\partial_{\rho},\\
Y_{\pm}&=&(\mp \rho A(\rho)+u)^{-1}\partial_{t}+\partial_{x}.
\end{eqnarray*}
In the case when both distributions $\mathcal{C}_{+}$ and $\mathcal{C}_{-}$ are integrable, one can explicitly solve the Cauchy problem (see, for example, \cite{KLR}). The integrability conditions for $\mathcal{C}_{+}$ and $\mathcal{C}_{-}$ are given by the following theorem \cite{LR-ljm-shock}:
\begin{theorem}
Distributions $\mathcal{C}_{+}$ and $\mathcal{C}_{-}$ are integrable if
\begin{equation}
\label{intr}
p(\rho)=c_{0}\rho^{3}+c_{1},
\end{equation}
where $c_{0}$ and $c_{1}$ are constants.
\end{theorem}

Nevertheless, we will not specify the dependence $p(\rho)$. To construct solutions in general case, we will need the following theorem \cite{KLR}:
\begin{theorem}
\label{thm-1}
A two-dimensional manifold $N\subset E$ is a multivalued solution to $\mathcal{E}$, if and only if the tangent spaces $T_{a}N$ for all $a\in N$ have one-dimensional intersections $h_{+}(a)$ and $h_{-}(a)$ with planes $\mathcal{C}_{+}(a)$ and $\mathcal{C}_{-}(a)$:
\begin{equation*}
T_{a}N=h_{+}(a)\oplus h_{-}(a).
\end{equation*}
\end{theorem}
The lines $h_{\pm}(a)=\mathcal{C}_{\pm}(a)\cap T_{a}N$ are called \textit{characteristic directions}. The statement of the above theorem is graphically shown in Fig.~\ref{char-dir}.

\begin{figure}[h!]
\centering
\includegraphics[scale=.5]{char1}
\caption{Characteristic directions.}
\label{char-dir}
\end{figure}

We will look for a one-parametric family of solutions defined by some 3-dimensional manifold
\begin{equation*}
N_{0}=\left\{F(t,x,u,\rho)=0\right\}
\end{equation*}
with integrable distribution $h=\langle h_{+},h_{-}\rangle$.

Let us take two vector fields $V_{\pm}\in\mathcal{C}_{\pm}$, tangent to $N_{0}$, i.e. $V_{\pm}\rfloor dF=0$ and choose $F(t,x,u,\rho)$ in such a way that $V_{+}\wedge V_{+}\wedge[V_{+},V_{-}]=0$. Restricting $V_{\pm}$ the $N_{0}$, we get an integrable distribution $h=\langle h_{+},h_{-}\rangle$, where $h_{+}=V_{+}|_{N_{0}}$ and $h_{-}=V_{-}|_{N_{0}}$.

Consider a particular case $F(t,x,\rho,u)=t-f(\rho,u)/\rho$. Then the integrability condition $V_{-}\wedge V_{+}\wedge[V_{-},V_{+}]=0$ leads to the following equation for $f(\rho,u)$:
\begin{equation}
\label{wave}
f_{uu}-A^{-2}(\rho)f_{\rho\rho}=0,
\end{equation}
which in some cases can be reduced to the wave equation with constant coefficients.
\begin{theorem}
Equation (\ref{wave}) is equivalent to the wave equation with constant coefficients if
\begin{equation*}
A(\rho)=(\beta_{1}\rho+\beta_{2})^{-2},
\end{equation*}
where $\beta_{1}$, $\beta_{2}$ are constants.
\end{theorem}
In this case we are able to find all 3-dimensional manifolds $N_{0}$ of the form $N_{0}=\left\{t-f(\rho,u)/\rho=0\right\}$.

Taking any of solutions to \eqref{wave} we get the corresponding integrable distribution $h=\langle h_{+},h_{-}\rangle$, and according to Theorem~\ref{thm-1}, its integrals give us multivalued solutions.
\subsection{Integrability via quotients}
There is another interpretation of \eqref{wave} as a quotient PDE for the system $\mathcal{E}$. To show this, let us first recall some concepts from the geometrical theory of PDEs following \cite{KLych,KrVin,KVL}.

Let $(x_{1},\ldots,x_{n})$ be a collection of independent variables and let $(u^{1},\ldots,u^{m})$ be a collection of dependent variables. Let $\mathcal{E}\subset J^{k}(\mathbb{R}^{n})$ be a system of PDEs with the symmetry Lie algebra $\mathfrak{g}$. By $\mathcal{E}_{l}\subset J^{l}(\mathbb{R}^{n})$ for some $l\ge k$ we will denote the prolongation of $\mathcal{E}$. Then, under some conditions, the field of rational differential $\mathfrak{g}$-invariants is finitely generated. More precisely, the global Lie-Tresse theorem \cite{KLych} is valid.
\begin{theorem}[Kruglikov, Lychagin]
Let $\mathcal{E}_{k}\subset J^{k}(\mathbb{R}^{n})$ be an algebraic formally integrable differential equation and let $\mathfrak{g}$ be its algebraic symmetry Lie algebra. Then, there exist rational differential $\mathfrak{g}$-invariants $a_{1},\ldots,a_{n},b^{1},\ldots,b^{N}$ of order $\le l$, such that the field of rational $\mathfrak{g}$-invariants is generated by rational functions of these functions and Tresse derivatives $\frac{d^{|\alpha|}b^{j}}{da^{\alpha}}$.
\end{theorem}
Remark that conditions of the global Lie-Tresse theorem are not very restrictive, and its statement is true for a considerable number of PDE systems and their symmetry Lie algebras, in particular, for Euler equations.

The Tresse derivatives mentioned in the Lie-Tresse theorem are constructed as follows (we omit here some technical details while emphasizing on general concepts, and we refer to \cite{KLych} for details). Let us take $n$ horizontally independent differential invariants $a_{1},\ldots,a_{n}$, which means that
\begin{equation}
\label{indep-cond}
\widehat{d}a_{1}\wedge\ldots\wedge\widehat{d}a_{n}\ne0
\end{equation}
in some Zariski-open set in $J^{l+1}(\mathbb{R}^{n})$. In \eqref{indep-cond}, $\widehat{d}$ is the total differential:
\begin{equation*}
\widehat{d}f=\sum\limits_{i=1}^{n}\frac{df}{dx_{i}}dx_{i},
\end{equation*}
where
\begin{equation*}
\frac{d}{dx_{i}}=\frac{\partial}{\partial x_{i}}+\sum\limits_{j,\sigma}u^{j}_{\sigma i}\frac{\partial}{\partial u_{\sigma}^{j}}.
\end{equation*}
Tresse derivatives are constructed as partial derivatives with respect to invariants $a_{1},\ldots,a_{n}$:
\begin{equation*}
\frac{d}{da_{i}}=\sum\limits_{j}A_{ij}\frac{d}{dx_{j}},\quad \frac{da_{j}}{da_{i}}=\delta_{ij},
\end{equation*}
where $\delta_{ij}$ is the Kronecker delta. Condition \eqref{indep-cond} guarantees the existence of solution to the system of linear equations on $A_{ij}$. By applying the Tresse derivatives to differential invariants we get new differential invariants.

In general, the algebra of differential invariants is not freely generated, there are relations between invariants, called \textit{syzygies}.

Let $s$ be a solution to $\mathcal{E}$, and consider restrictions:
\begin{equation*}
a_{1}|_{s},\ldots,a_{n}|_{s},b^{1}|_{s},\ldots,b^{N}|_{s},
\end{equation*}
which locally can be viewed as $(n+N)$ functions on an $n$-dimensional manifold, therefore (locally)
\begin{equation*}
b^{j}|_{s}=B^{j}(a_{1}|_{s},\ldots,a_{n}|_{s}),\quad j=1,\ldots,N.
\end{equation*}
In fact, $B^{j}$ depend on the equivalence class of $s$ (where the equivalence relation is defined by the symmetry Lie algebra), rather than on $s$ itself.

Removing restrictions, we get:
\begin{equation}
\label{rel}
b^{j}=B^{j}(a_{1},\ldots,a_{n}).
\end{equation}
Functions $B^{j}$ can be found from \textit{quotient PDEs}. Let us apply the Tresse derivatives to \eqref{rel}. We get
\begin{equation*}
J_{ji}=\frac{db^{j}}{da_{i}}=B^{j}_{a_{i}}(a_{1},\ldots,a_{n}),\ldots,J_{j\alpha\beta}=\frac{d^{2}b^{j}}{da_{\alpha}da_{\beta}}=B^{j}_{a_{\alpha}a_{\beta}}(a_{1},\ldots,a_{n}),\ldots
\end{equation*}
Since $b^{j}$ are invariants, $\frac{d}{da_{i}}$ are invariant derivations, functions $J_{ji},J_{j\alpha\beta},\ldots$ are invariants too. Finding syzygies between invariants $a_{i}$, $b^{j}$, $J_{ji}$, $J_{j\alpha\beta},\ldots$ we get a relation (perhaps, a number of them)
\begin{equation*}
F(a_{1},\ldots,a_{n},B^{1},\ldots,B^{N},B^{j}_{a_{i}},B^{j}_{a_{\alpha}a_{\beta}},\ldots)=0,
\end{equation*}
called a \textit{quotient PDE}, which is a PDE on functions $B^{j}(a_{1},\ldots,a_{n})$.

Let us collect some of the most important properties of quotient PDEs:
\begin{enumerate}
\item Solution to a quotient PDE is a $\mathfrak{g}$-orbit of $s$;
\item Solutions to quotient PDEs provide us with differential constraints $b^{j}-B^{j}(a_{1},\ldots,a_{n})=0$, compatible with the original system $\mathcal{E}$.
\end{enumerate}
Finding such constraints, we reduce the integration of $\mathcal{E}$ to the integration of a completely integrable Cartan distribution $\mathcal{C}_{\mathcal{E}}$ with the same symmetry algebra.

Let us now apply these ideas to the Euler system $\mathcal{E}$.
\begin{theorem}
The symmetry Lie algebra of the system of Euler equations \eqref{EU} is generated by vector fields
\begin{equation*}
\begin{split}
X_{1}=\frac{\partial}{\partial x},\quad &X_{2}=\frac{\partial}{\partial t},\quad X_{3}=t\frac{\partial}{\partial x}+\frac{\partial}{\partial u},\quad X_{4}=t\frac{\partial}{\partial t}+x\frac{\partial}{\partial x},\quad {}\\&
X_{5}=C_{1}e^{-C_{1}u}Y\frac{\partial}{\partial t}+(C_{1}Yu+Y^{\prime}\rho+Y)e^{-C_{1}u}\frac{\partial}{\partial x},
\end{split}
\end{equation*}
where $Y=Y(\rho)$ is a solution to the ODE
\begin{equation*}
-\rho A^{2}(\rho)C_{1}^{2}Y+\rho Y^{\prime\prime}+2Y^{\prime}=0,
\end{equation*}
and $C_{1}$ is a constant.
\end{theorem}
Let us choose $\mathfrak{g}=\langle\partial_{x}\rangle$. Then, we have the following 0-order invariants:
\begin{equation*}
a_{1}=u,\quad a_{2}=\rho,\quad b^{1}=t.
\end{equation*}
Put $b^{1}=t=B(u,\rho)$. Then, Tresse derivatives will be:
\begin{eqnarray*}
\frac{d}{du}&=&\frac{\rho_{x}}{\rho u_{x}^{2}-\rho A^{2}\rho_{x}^{2}}\frac{d}{dt}+\frac{\rho u_{x}+u\rho_{x}}{\rho u_{x}^{2}-\rho A^{2}\rho_{x}^{2}}\frac{d}{dx},\\
\frac{d}{d\rho}&=&\frac{u_{x}}{\rho A^{2}\rho_{x}^{2}-\rho u_{x}^{2}}\frac{d}{dt}+\frac{A^{2}\rho\rho_{x}+uu_{x}}{\rho A^{2}\rho_{x}^{2}-\rho u_{x}^{2}}\frac{d}{dx}.
\end{eqnarray*}
Finding syzygies between $u$, $\rho$, $\frac{dB}{du}$, $\frac{dB}{d\rho}$, $\frac{d^{2}B}{du^{2}}$, $\frac{d^{2}B}{d\rho^{2}}$, $\frac{d^{2}B}{du d\rho}$, we get a quotient PDE:
\begin{equation*}
A^{2}\rho\frac{d^{2}B}{du^{2}}-\rho\frac{d^{2}B}{d\rho^{2}}-2\frac{dB}{d\rho}=0.
\end{equation*}
Putting $B(u,\rho)=f(u,\rho)/\rho$, we get
\begin{equation}
\label{QUO-PDE}
f_{uu}-A^{-2}f_{\rho\rho}=0,
\end{equation}
which coincides with \eqref{wave}.

Let us take any solution $f(u,\rho)$ to \eqref{QUO-PDE}. Then, differentiating $t-f(u,\rho)/\rho=0$ with respect to $t$ and $x$, we get two more PDEs, and together with Euler equations we get 5 relations on 1-jet space, determining the submanifold $\widehat{\mathcal{E}}\subset J^{1}(\mathbb{R}^{2})$.

So, $\dim J^{1}(\mathbb{R}^{2})=8$, $\dim\widehat{\mathcal{E}}=3$, $\dim\mathcal{C}_{\widehat{\mathcal{E}}}=2\Longrightarrow\mathrm{codim}\,\mathcal{C}_{\widehat{\mathcal{E}}}=1$, and the Cartan distribution $\mathcal{C}_{\widehat{\mathcal{E}}}$ is given by a differential 1-form
\begin{equation}
\label{Cart-form}
\begin{split}
\omega=&-\frac{\rho  (f_{\rho} \rho -f)}{\rho^2 A^2 f_{u}^2-f_{\rho}^2 \rho^2+2 f_{\rho} f\rho-f^2}dx+\frac{(f_{\rho} \rho -f) (f_{u} u -f_{\rho} \rho +f)}{\rho^2 A^2 f_{u}^2-f_{\rho}^2 \rho ^2+2 f_{\rho} f \rho -f^2}du+{}\\&+\frac{(-\rho ^2 A^2 f_{u}+f_{\rho} \rho  u -f u ) (f_{\rho} \rho -f)}{\rho^3 A^2 f_{u}^2-\rho^3 f_{\rho}^2+2 \rho^2 f f_{\rho}-f^2 \rho}d\rho,
\end{split}
\end{equation}
such that $\omega\wedge d\omega=0$ for any solution $f(u,\rho)$ to \eqref{wave}. Integrals of the form $\omega$ give us solutions to $\widehat{\mathcal{E}}$.
\begin{remark}
The distribution generated by $\omega$ coincides with that generated by $\langle h_{-},h_{+}\rangle$.
\end{remark}
\subsection{Solutions}
We can see that finding differential constraint to integrate the Euler system can be performed in one of two equivalent ways described in this paper. We will proceed with the integration of $\mathcal{E}$ by finding integrals of the form \eqref{Cart-form}.

Note that the vector field $\partial_{x}$ is an infinitesimal transversal symmetry of \eqref{Cart-form}, and therefore the differential form $\varkappa=\omega/\omega(\partial_{x})$ is closed and therefore locally exact, i.e. $\varkappa=dJ$ for some function $J\in C^{\infty}(E)$, and $J=\mathrm{const}$ together with $t-f(u,\rho)/\rho=0$ give us solutions to $\mathcal{E}$.

Let us apply the following separation ansatz to \eqref{wave}: $f(u,\rho)=\rho\mu(u)+\nu(\rho)$. Then, we get two ODEs on $\mu(u)$ and $\nu(\rho)$:
\begin{equation*}
\mu^{\prime\prime}=\lambda,\quad \nu^{\prime\prime}=\lambda\rho A^{2},
\end{equation*}
for some constant $\lambda$, and we get the first quadrature for $\mathcal{E}$:
\begin{equation}
\label{sol-t}
t=\frac{\lambda u^{2}}{2}+\alpha_{0}u+\frac{\alpha_{2}}{\rho}+\frac{\lambda}{\rho}\int Q(\rho)d\rho+t_{0},
\end{equation}
where $\alpha_{j}, t_{0}$ are constants,  $Q(\rho)=\int\rho A^{2}d\rho$.

The differential 1-form $\varkappa$ equals
\begin{equation}
\label{varkappa}
\begin{split}
\varkappa=&dx-\left(\lambda u^{2}+\alpha_{0}u-\lambda Q(\rho)+\frac{\lambda}{\rho}\int Q(\rho)d\rho+\frac{\alpha_{2}}{\rho}\right)du+{}\\&+\left(\frac{\lambda u}{\rho^{2}}\int Q(\rho)d\rho-\frac{\lambda u}{\rho}Q(\rho)+\rho A^{2}(\rho)(\lambda u+\alpha_{0})+\frac{\alpha_{2}u}{\rho^{2}}\right)d\rho.
\end{split}
\end{equation}
Integrating \eqref{varkappa}, we get the second quadrature for $\mathcal{E}$:
\begin{equation}
\label{sol-x}
x=\frac{\lambda u^{3}}{3}-\lambda uQ(\rho)+\frac{\alpha_{0}u^{2}}{2}+\frac{\lambda u}{\rho}\int Q(\rho)d\rho+\frac{\alpha_{2}u}{\rho}-\alpha_{0}Q(\rho)+x_{0},
\end{equation}
where $x_{0}$ is a constant.

Formulae \eqref{sol-t},\eqref{sol-x} give us a multivalued solution to $\mathcal{E}$. Note that it can be applied to any thermodynamic state model, because $p(\rho)$ is not specified.

Singularities of projections of the multivalued solution $N$ given by \eqref{sol-t},\eqref{sol-x} to the space of independent variables appear where the differential 2-form $dx\wedge dt$ degenerates. Equation $(dx\wedge dt)|_{N}=0$ gives us a curve, called \textit{caustic}. Parametric equations for the caustic are
\begin{eqnarray*}
x(\rho)&=&x_{0}-\frac{Z_{\pm}^{3}}{3\lambda^{2}A^{3}\rho^{6}}+\frac{\alpha_{0}Z_{\pm}^{2}}{2\lambda^{2}A^{2}\rho^{4}}\pm\frac{\alpha_{0}Z_{\pm}}{\lambda\rho}\mp\frac{Z_{\pm}^{2}}{\lambda A\rho^{3}}-\alpha_{0}Q,\\
t(\rho)&=&\frac{Z_{+}^{2}}{2\lambda A^{2}\rho^{4}}+\frac{Z_{+}}{\rho}\left(1-\frac{\alpha_{0}}{\lambda\rho A}\right)-\alpha_{0}\rho A+\lambda Q+t_{0},
\end{eqnarray*}
where
\begin{equation*}
Z_{\pm}=\alpha_{0}\rho^{2}A\mp\lambda\rho Q\pm\lambda\int Qd\rho\pm\alpha_{2}.
\end{equation*}

To cut a multivalued solution given by \eqref{sol-t},\eqref{sol-x} into single valued branches representing the corresponding discontinuous solution, one needs a conservation law. Following \cite{Lan-Lif}, we will use the mass conservation law. Let us choose $(t,\rho)$ as coordinates on the solution $N$. Resolving \eqref{sol-t} with respect to $u$, we get two roots
\begin{equation}
\label{u-rho-t-sol}
u=U_{\pm}(\rho,t)=-\frac{\alpha_{0}}{\lambda}\pm\frac{1}{\lambda\rho}\sqrt{2\rho\left(\rho\lambda(t-t_{0})+\frac{\rho\alpha_{0}^{2}}{2}-\lambda\alpha_{2}-\lambda^{2}\int Qd\rho\right)}.
\end{equation}

Substituting this relation to \eqref{mass}, we get
\begin{equation*}
\rho_{t}+(\rho U_{\pm})_{x}=0,
\end{equation*}
and therefore the conservation law is
\begin{equation}
\label{con-law}
\Theta_{\pm}=\rho dx-\rho U_{\pm}dt.
\end{equation}
The restriction $\Theta_{\pm}|_{N}$ of \eqref{con-law} to solution \eqref{sol-t},\eqref{sol-x} is a closed form, and therefore locally $\Theta_{\pm}|_{N}=dH_{\pm}$, and the potential $H_{\pm}(\rho,t)$ equals
\begin{equation*}
H_{\pm}(\rho,t)=\mp\frac{1}{\lambda\sqrt{\rho}}\sqrt{2\left(\rho\lambda(t-t_{0})+\frac{\rho\alpha_{0}^{2}}{2}-\lambda\alpha_{2}-\lambda^{2}\int Qd\rho\right)}\left(\rho\lambda Q-\lambda\int Q d\rho-\alpha_{2}\right).
\end{equation*}

The points where solution has discontinuity (shock wave front) is found from the system of equations:
\begin{equation*}
H_{\pm}(\rho_{1},t)=H_{\pm}(\rho_{2},t), \quad x(\rho_{1},t)=x(\rho_{2},t),
\end{equation*}
which is to be resolved at each time moment $t$ with respect to $\rho_{1,2}$ --- the values of densities on each branch of a discontinuous solution, and where $x(\rho,t)$ is found as a result of substitution of \eqref{u-rho-t-sol} to \eqref{sol-x}.

\section{Solutions for an ideal gas, caustics, shock waves}
Here, we illustrate solutions obtained in the previous section on an ideal gas model. The Legendrian manifold for an ideal gas is given by
\begin{equation*}
\widehat{L}=\left\{p=R\rho T,\, e=\frac{n}{2}RT,\, s=R\ln\left(\frac{T^{n/2}}{\rho}\right)\right\},
\end{equation*}
where $R$ is the universal gas constant, and $n$ is the degree of freedom.

The differential quadratic form $\kappa|_{L}$ for ideal gases is
\begin{equation*}
\kappa|_{L}=-\frac{Rn}{2}\frac{dT^{2}}{T^{2}}-R\rho^{-2}d\rho^{2},
\end{equation*}
and it is negative on the entire Lagrangian manifold. Therefore, the system $\mathcal{E}$ in case of ideal gases is hyperbolic for any process $l\subset\widehat{L}$.

By fixing the entropy constant $s=s_{0}$, we get the following expressions for $T(\rho)$ and $p(\rho)$:
\begin{equation*}
T(\rho)=\exp\left(\frac{2s_{0}}{Rn}\right)\rho^{2/n},\quad p(\rho)=R\exp\left(\frac{2s_{0}}{Rn}\right)\rho^{2/n+1},
\end{equation*}
and the function $A(\rho)=\rho^{-1}\sqrt{p^{\prime}(\rho)}=A_{0}\rho^{m}$ (see also \cite{LR-ljm-shock}), where
\begin{equation*}
A_{0}=\sqrt{R\left(1+\frac{2}{n}\right)\exp\left(\frac{2s_{0}}{Rn}\right)},\quad m=\frac{1}{n}-1.
\end{equation*}

In the case of an ideal gas, we get the following quadratures:
\begin{equation*}
t=\frac{\lambda u^{2}}{2}+\alpha_{0}u+\frac{\alpha_{2}}{\rho}+\frac{\lambda A_{0}^{2}\rho^{2m+2}}{2(2m^{2}+5m+3)}+t_{0},
\end{equation*}
\begin{equation*}
\begin{split}
x=&\frac{1}{\rho(2m^{2}+5m+3)}\left((m+1)(2m+3)\left(\alpha_{2}u+\rho\left(\frac{\lambda u^{3}}{3}+\frac{\alpha_{0}u^{2}}{2}+x_{0}\right)\right)-\right.{}\\&-\left.A_{0}^{2}\rho^{2m+3}\left(m(\lambda u+\alpha_{0})+\lambda u+\frac{3\alpha_{0}}{2}\right)\right).
\end{split}
\end{equation*}

Equations for the caustics are
\begin{eqnarray*}
x(\rho)&=&x_{0}-\frac{Z_{\pm}^{3}}{3\lambda^{2}A_{0}^{3}\rho^{3m+6}}+\frac{\alpha_{0}Z_{\pm}^{2}}{2\lambda^{2}A_{0}^{2}\rho^{2m+4}}\pm\frac{\alpha_{0}Z_{\pm}}{\lambda\rho}\mp\frac{Z_{\pm}^{2}}{\lambda A_{0}\rho^{m+3}}-\frac{\alpha_{0}A_{0}^{2}\rho^{2m+2}}{2(m+1)},\\
t(\rho)&=&\frac{Z_{+}^{2}}{2\lambda A_{0}^{2}\rho^{2m+4}}+\frac{Z_{+}}{\rho}\left(1-\frac{\alpha_{0}}{\lambda A_{0}\rho^{m+1}}\right)-\alpha_{0}A_{0}\rho^{m+1}+\frac{\lambda A_{0}^{2}\rho^{2m+2}}{2(m+1)}+t_{0},
\end{eqnarray*}
where
\begin{equation*}
Z_{\pm}=\alpha_{0}A_{0}\rho^{m+2}\mp\frac{\lambda A_{0}^{2}\rho^{2m+3}}{2m+3}\pm\alpha_{2}.
\end{equation*}

The graphs of the multivalued solutions are presented in Fig.~\ref{mult}, where we used substitution $A_{0}=1,\, \alpha_{0}=1,\,t_{0}=1,\,\alpha_{2}=-2,\, x_{0}=0$.

\begin{figure}[h]

\begin{minipage}[h]{0.32\linewidth}
\center{\includegraphics[scale=0.18]{init}}
\end{minipage}
\hfill
\begin{minipage}[h]{0.32\linewidth}
\center{\includegraphics[scale=0.18]{mom2}}
\end{minipage}
\hfill
\begin{minipage}[h]{0.32\linewidth}
\center{\includegraphics[scale=0.18]{mom3}}
\end{minipage}
\caption{Density graph for $n=3$ at moments $t=0$, $t=2.7$, $t=3.75$}
\label{mult}
\end{figure}

The discontinuity lines together with caustics are shown in Fig.~\ref{shock}. One can observe a very interesting phenomenon. Smooth initial datum evolutes into formation of two cusps, moving towards each other. However, as numerical computations show, two discontinuities meet only asymptotically when $t\to+\infty$.

\begin{figure}[h!]
\centering
\includegraphics[scale=.5]{wave}
\caption{Caustics (black) and shock wave fronts (red).}
\label{shock}
\end{figure}

\section*{Acknowledgements}
The author is grateful to Valentin Lychagin for helpful suggestions during the preparation of the paper. This work was partially supported by the Foundation for the Advancement of Theoretical Physics and Mathematics ``BASIS'' (project 19-7-1-13-3) and by the Russian Science Foundation (project 21-71-20034).