They have not provided any derivation of the statement for hydrodynamic or MHD flow. It would be great if one could provide a derivation of the above statement for hydrodynamic or MHD flow.
Perhaps a better reference is a review by Holzer [1989].
Start by looking at a stead-state pressure balance of a fluid flow given by:
$$
\rho u^{2} + P + \frac{ B^{2} }{ 2 \mu_{o} } + P_{f} = constant \tag{0}
$$
where $\rho$ is the mass density of the flow, $u$ is the bulk flow velocity normal to some boundary of interest, $B$ is the magnetic field magnitude, $\mu_{o}$ is the permeability of free space, $P$ is the thermal pressure, and $P_{f}$ is a term used to encapsulate the pressure contributions from things like neutral gases and cosmic rays.
In general the first term (i.e., the dynamic pressure) in Equation 0 will dominate all of these terms except for very close to the Sun. We can assume an asymptotic, steady-state flow where the speed, $u$, is constant but the mass density, $\rho$, drops as $r^{-2}$ (e.g., like adiabatic expansion). So we let $\rho = \rho_{E} r^{-2}$ and $u = u_{E}$, where subscript $E$ corresponds to values measured near Earth. Then we can equate the insterstellar pressure, $P_{I}$, with the ram pressure to get:
$$
R_{t} = \sqrt{ \frac{ \rho_{E} u_{E}^{2} }{ P_{I} } } \tag{1}
$$
We can take this further and assume the interaction produces a strong shock wave, which from the Rankine–Hugoniot relations we know that:
$$
\begin{align}
\frac{ u_{up} }{ u_{dn} } & = \frac{ \rho_{dn} }{ \rho_{up} } = 4 \tag{2a} \\
P_{dn} & = \frac{ 3 \rho_{up} u_{up}^{2} }{ 4 } \tag{2b}
\end{align}
$$
where the subscript $up$($dn$) refers to values upstream(downstream) of the shock. Note to derive Equation 2b we use Equation 0 but ignore contributions from the magnetic field and $P_{f}$ in both regions, which is generally okay (not the greatest, but okay for now). We also ignore the thermal pressure in the upstream as being small, such that one has:
$$
\rho_{up} u_{up}^{2} \approx \rho_{dn} u_{dn}^{2} + P_{dn} \tag{3}
$$
We then replace $\rho_{dn}$ and $u_{dn}$ using Equation 2a to get Equation 2b.
Further downstream of the termination shock there will be a point of stagnation flow, which is called the heliopause. After this boundary we can assume that thermal dominates, so if Equation 0 is to hold we must have:
$$
P_{I} \approx \rho_{dn} u_{dn}^{2} + P_{dn} \tag{4}
$$
Note that the more accurate representation includes a factor of $\tfrac{ 1 }{ 2 }$ with the dynamic pressure. So taking Equation 4 and using Equations 2a and 2b, we find that:
$$
P_{I} = \frac{ 7 }{ 8 } \rho_{up} u_{up}^{2} \tag{5}
$$
If we follow a similar line of reasoning that led us to Equation 1 for adiabatic, stationary flow, we can show the location of a shock, $R_{s}$, occurs at:
$$
R_{s} = \sqrt{ \frac{ 7 \rho_{E} u_{E}^{2} }{ 8 P_{I} } } \tag{6}
$$
which only differs from Equation 1 by ~12.5%.
...I do not understand why reverse shock will appear in the system.
A reverse shock forms when an obstacle stands against a supersonic flow, thus the outward normal of the shock is anti-parallel to the incident flow. In contrast, an interplanetary shock moves with the solar wind so the outward normal is parallel to the flow, thus is considered a forward shock.
Also why a contact discontinuity separates those two shocked medium. I think the contact discontinuity is same as the stagnation point in the paper though I have to read it further in order to understand it.
Yes, a contact discontinuity is a stagnation point where no mass flux or flow go across some boundary. In plasmas, the stagnation point is actually formed by something called a tangential discontinuity. The jump conditions for these are:
$$
\begin{align}
\left\{ \mathbf{u}_{t} \right\} & \neq 0 \tag{7a} \\
\left\{ \mathbf{B}_{t} \right\} & \neq 0 \tag{7b} \\
\left\{ \rho \right\} & \neq 0 \tag{7c} \\
U_{n} & = 0 \tag{7d} \\
B_{n} & = 0 \tag{7e} \\
\left\{ P + \frac{ B^{2} }{ 2 \mu_{o} } \right\} & = 0 \tag{7f}
\end{align}
$$
where subscript $t$($n$) denotes the transverse(normal) direction relative to the shock normal unit vector and the $\{ Q \} = Q_{dn} - Q_{up}$ denote the change in the contained parameters across the boundary. In fluid dynamics, yes the stagnation point would be like a contact discontinuity but in MHD the stagnation points behind shocks are tangential discontinuities.
As to why a stagnation point is even necessary, the flow must move around the obstacle causing the shock in the first place, not through. Therefore, the normal flow must stop at some point and redirect along the surface of the stagnation to get around the obstacle otherwise there'd be a diverging mass pile-up at the obstacle.
Note that under some conditions, the boundary between two regions can actually have flow across it (e.g., open magnetosphere), i.e., $U_{n} \neq 0$, then the boundary is called a rotational discontinuity. Unlike the tangential discontinuity, there is no density jump across a rotational discontinuity. Basically, the magnetic field and bulk flow velocity transverse components just rotate. In fact, a rotational discontinuity is sometimes considered a special type of Alfven wave.
So the heliopause can behave like both of these, depending on the conditions at any given location.
References
