---
abstract: 'It is proposed to evaluate the effects of thermodynamic metastability on fluid dynamics by comparing two different ideal-hydrodynamics solutions — one obtained with the fully equilibrium equation of state using the Maxwell construction, and the other in what we call the phase-flip approximation.  The latter is based on the assumption of instantaneous decay of metastable states upon reaching the spinodal.  The proposed method is applied to the classical problem of the centered rarefaction wave by expansion into vacuum, for which exact analytical solutions exist in both approximations.  It is shown that the rapid decay of metastable states leads to the formation of a rarefaction shock in the expanding flow.  Implications for the laser-heating experiments are discussed.' 
author:
- 'Mikhail M. Basko'
bibliography:
- 'Basko\_pfhyd18-biblio.bib'
title: 'Centered rarefaction wave with a liquid-gas phase transition in the approximation of “phase-flip” hydrodynamics'
---

Introduction
============

A typical situation in the majority of laser-matter experiments is where a laser beam heats up a surface layer of a solid (liquid) sample, which then begins to expand into vacuum or a low-pressure environment.  In many cases, important for both practical applications and fundamental research, the expanding material passes through the liquid-gas phase coexistence region in the thermodynamic parameter space.  As characteristic examples, two different groups of experiments with ultra-short laser pulses [@Linde_Sokow-Tinten.1997; @Vinokhodov_Koshelev.2016] can be pointed out.  In both these examples, an adequate description of fluid dynamics in the two-phase region proved to be crucial for the correct interpretation of the observed results [@Sokol-Tint_Bialkowski.1998; @Basko_Krivokor.2017]. 

Any attempt to employ the one-fluid hydrodynamics for modeling such experiments stumbles over a dilemma of what equation of state (EOS) ought to be used in the two-phase region below the critical point of the liquid-gas phase transition.  One approach, which is fully self-consistent and poses neither conceptual nor technical difficulties, is to employ the *fully equilibrium* (EQ) EOS, obtained by applying the well-known Maxwell rule.  The EQ EOS was successfully used [@Sokol-Tint_Bialkowski.1998; @Inogamov_Anisimov.1999] to explain the Newton rings observed in Ref.  , but only partially so by modeling [@Basko_Krivokor.2017] the fragmentation dynamics of liquid tin microspheres, hit by picosecond laser pulses [@Vinokhodov_Koshelev.2016]. 

At the same time, an obvious flaw of the EQ hydrodynamics is the tacit assumption that the boiling of liquid (condensation of vapor) begins immediately upon crossing the outer boundary of the phase coexistence region — the *binodal*.  In reality, especially in the laser-driven experiments where the hydrodynamic time scale can be as short as ${\:\raisebox{.25ex}{$<$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}0.1$ ns, the fluid elements can penetrate rather deeply into the metastability region and boil up (condense) quite close to the boundary of absolute thermodynamic instability — the *spinodal*.  To remedy this deficiency of the EQ hydrodynamics, we propose an approach, based on comparing the EQ solution for a given problem with an alternative hydrodynamic solution for the same problem, obtained in what we call the *phase-flip* (PF) approximation. 

In the PF approximation we assume that each fluid element obeys the metastable (MS) EOS all the way down to the spinodal, where an instantaneous and irreversible relaxation to the EQ EOS takes place.  The justification for this assumption is given in Section \[s:pf\].  On the one hand, the PF approximation is simpler than other models invoking the kinetics of phase transition (see, for example, Ref.   and references therein) because one is allowed to stay within the framework of the one-fluid ideal hydrodynamics.  On the other hand, allowing for maximum possible penetration into the region of metastable states, it can be expected to provide a good estimate for the maximum possible influence of metastability on the flow pattern under specific experimental conditions. 

This paper is devoted to the application of the PF approach to one of the classical flows in the ideal hydrodynamics — the centered rarefaction wave by expansion into vacuum [@Courant_Friedrichs1977; @LL-H87; @Zeldovich_Raizer2012].  Although the principal ideas behind the PF approximation have been discussed earlier [@Saurel_Petitpas.2008], it is, to the best of our knowledge, for the first time that an exact analytical PF solution is presented in full detail (Section \[s:crw\]) and applied to draw specific qualitative and quantitative conclusions in regard to certain types of laser-driven experiments.  This is done by making use of a specific form of the generalized van der Waals EOS introduced in Section \[s:GW\], which is reasonably well suited for the description of real materials. 

Liquid-gas phase diagram for the generalized van der Waals equation of state \[s:GW\]
=====================================================================================

All the numerical results in this paper are obtained with a generalized version of the van der Waals equation of state (GWEOS), where the power exponent $n$ in the attractive term is treated as a free parameter [@Martynyuk1991; @Martynyuk1993; @Basko2018].  This variant of a two-phase EOS preserves the convenient property of the original van der Waals EOS that it can be cast in the universal reduced form, where the dimensionless quantities $\rho$, $\theta$ and $p$ are, respectively, the density, temperature and pressure normalized by their values $\rho_{cr}$, $T_{cr}$ and $P_{cr}$ at the critical point; the latter can be considered as three independent free parameters of GWEOS taken, for example, from experiment. 

In its reduced form, GWEOS can be presented as $$\begin{aligned}
\label{GW:p(v,tet)=}
  p(v,\theta) &=& \frac{\alpha \theta}{v-\kappa^{-1}} -\frac{\kappa}{v^n},
  \\ \label{GW:e(v,tet)=}
  e(v,\theta) &=& c_V\alpha \theta -\frac{1}{2} \kappa(\kappa-1) v^{1-n},
  \\ \label{GW:s(v,tet)=}
  s(v,\theta) &=& \alpha \left[ c_V\left(1+ \ln\theta\right)
  +\ln\left(v-\kappa^{-1}\right) \right],\end{aligned}$$ where $v = \rho^{-1}$ is the reduced specific volume, $e$ is the specific internal energy normalized by $P_{cr}/\rho_{cr}$, $s$ is the specific entropy normalized by $P_{cr}/\left(\rho_{cr}T_{cr}\right)$, and $$\label{GW:alpha=}
  c_V >0, \quad \alpha=\kappa-\kappa^{-1}, \quad \kappa =\frac{n+1}{n-1} > 1$$ are dimensionless constants.  The square of the reduced sound velocity, normalized by $P_{cr}/\rho_{cr}$, is given by $$\label{GW:cs2=}
  c_s^2 \equiv \left(\frac{\partial p}{\partial \rho}\right)_s =
  \frac{\left(1+c_V^{-1}\right) \alpha\, \theta v^2}{\left(v-\kappa^{-1} \right)^2} -\frac{n\kappa}{v^{n-1}}.$$ Physically meaningful are the values $n>1$, implying $\kappa >1$. The original van der Waals EOS is recovered for $n=2$. Note that GWEOS is defined over a finite density interval $0 < \rho < \kappa$, i.e. for $v> \kappa^{-1}$.

Beside $\rho_{cr}$, $T_{cr}$ and $P_{cr}$, GWEOS has two more free parameters, namely, the heat capacity at constant volume $c_V$ and the power exponent $n$. Since $c_V$ is assumed to be constant, we use its ideal-gas value in the limit of $v\to \infty$; for the monoatomic substances $c_V=3/2$. A physically sensible recipe for evaluation of $n$ could be by fitting the experimental values of the dimensionless ratio $$\label{GW:Lam=}
  \Lambda = \frac{E_{coh}}{T_{cr}} = \frac{(n+1)^{n+1}}{4n(n-1)^n},$$ where $E_{coh}$ is the cohesive energy per atom (molecule).  For a wide variety of monoatomic and diatomic substances the measured $\Lambda$ values [@Grigoriev1997] fall in the range $\Lambda \approx 4.0$–5.3, which implies $n = 1.4$–1.65.  Somewhat different physical considerations [@Martynyuk_B1995] have led to practically the same values of $n$. Accordingly, in this work we adopt a universal value of $n=1.5$. Note that this is one aspect where GWEOS is superior to the original van der Waals EOS, which yields a not quite realistic value of $\Lambda =3.375$.

![( Color online) GWEOS with $n=c_V=1.5$ on the thermodynamic $(v,p)$ plane.  Metastable states (shaded cyan) lie between the binodal $bi$ (blue) and the spinodal $sp$ (magenta) curves.  The area between the spinodal $sp$ and the cold curve $\theta=0$ (shaded grey and grey-orange) is the region of absolute thermodynamic instability.  The region of non-hyperbolicity with $c_s^2 <0$ is shaded grey-orange.  Two isotherms $\theta=0.75$ and $\theta=0.92$ are shown as thin black curves in both the MS (solid) and EQ (dash-dotted) versions of GWEOS.  Note that for $p<0$ the vertical axis has a different scale \[linear in $\ln(1-p)$\], which explains the kinks on isotherms and the spinodal by crossing the $p=0$ line.  \[f:1\]](fig_1_pv.eps){width="\myfigwidth"}

Main features of GWEOS for $n=1.5$ ($\kappa=5$) are illustrated in Fig.  \[f:1\].  Below the critical point $CP$ at $p=v=\theta=1$ any isotherm (\[GW:p(v,tet)=\]) with $\theta < 1$ has a segment $v_{sp,l}< v< v_{sp,g}$ of absolute thermodynamic instability, where $\partial p(v,\theta)/\partial v >0$. On the $(v,p)$ and $(v,\theta)$ planes the region of absolute instability lies below the spinodal curve defined by $$\begin{aligned}
\label{GW:p_sp(v)=} &&
  p=p_{sp}(v) = \frac{(n+1)v-n}{v^{n+1}}, \\ \label{GW:tet_sp(v)=} &&
  \theta = \theta_{sp}(v) = \frac{1}{v^{n+1}}
  \left(\frac{v-\kappa^{-1}}{1-\kappa^{-1}} \right)^2.\end{aligned}$$ For any $p<1$, the liquid, $v= v_{sp,l}(p)< 1$, and the vapor, $v= v_{sp,g}(p)> 1$, branches of the spinodal are found as the two roots of Eq.  (\[GW:p\_sp(v)=\]).  According to the basic thermodynamic principles, thermodynamic states, calculated from Eqs.  (\[GW:p(v,tet)=\])–(\[GW:cs2=\]) inside the interval $v_{sp,l}< v< v_{sp,g}$, cannot be realized in nature. 

Equilibrium thermodynamics tells us [@LL-SP96 §81,83] that, for any given $\theta< 1$, the liquid-vapor phase transition occurs at a fixed pressure $p<1$ across a finite interval $v_{bi,l}< v< v_{bi,g}$ of the phase coexistence.  The functions $v= v_{bi,l}(p)< 1$ and $v= v_{bi,g}(p)> 1$, obtained by applying the Maxwell rule [@LL-SP96 §84,85] to isotherms (\[GW:p(v,tet)=\]), make up, respectively, the liquid and the vapor branches of the binodal curve $bi$ in Fig.  \[f:1\].  Accordingly, the Maxwell construction applied to Eqs.  (\[GW:p(v,tet)=\])–(\[GW:cs2=\]) yields the EQ version of GWEOS for the entire phase coexistence region below the binodal.  The full set of formulae for calculating the EQ GWEOS is given in Ref.  .  Note that one always has $v_{bi,l}(p)< v_{sp,l}(p)$ and $v_{bi,g}(p)> v_{sp,g}(p)$.

In the region between the binodal and spinodal curves (shaded cyan in Fig.  \[f:1\]) the original Eqs.  (\[GW:p(v,tet)=\])–(\[GW:cs2=\]) are still physically meaningful: they represent the metastable thermodynamic states and provide the alternative MS option of GWEOS in this region.  The region of non-hyperbolicity, where $c_s^2< 0$ and the hydrodynamic equations cease to be hyperbolic, lies always below the spinodal inside the region of absolute thermodynamic instability. 

![( Color online) Four expansion isentropes on the $(v,p)$ plane for GWEOS with $n=c_V=1.5$. Upon crossing the binodal at point $B$, each isentrope splits into the metastable MS (black thin solid) and the equilibrium EQ (black dash-dotted) branches.  Every MS branch finally reaches the spinodal at point $D$; its further continuation (black dots) into the instability region is assumed to be meaningless from the viewpoint of hydrodynamics.  \[f:2\]](fig_2_s-on-pv_v2.eps){width="\myfigwidth"}

In this work we analyze a centered rarefaction wave into vacuum from a certain initial state $(v_0,p_0)$, denoted as point $O$ in Figs.  \[f:2\]–\[f:4\].  If the flow remains continuous, it evolves along the expansion isentrope starting at $O$. Figure \[f:2\] shows four examples of such isentropes for the GWEOS with $n= c_V=1.5$. When the origin point $O$ lies above the binodal, the isentrope crosses the binodal at a certain point $B$ with $v=v_B> v_0$. In ideal hydrodynamics, further expansion to $v> v_B$ can proceed either along the EQ branch of this isentrope (dash-dotted curves in Fig.  \[f:2\]), or along its MS branch (thin solid curves) until the latter crosses the spinodal at a certain point $D$; we assume that below $D$ the MS branch of the expanding isentrope becomes physically meaningless.  An important property of the EQ isentropes inside the phase coexistence region is that they never come out from under the binodal. 

In GWEOS, the equation of an MS isentrope, passing through a certain point $(v_0,p_0)$ on the $(v,p)$ plane, is $$\label{GW:p_s(v)=}
  p = p_s(v) =\left(p_0 +\frac{\kappa}{v_0^n} \right) \left(\frac{v_0-\kappa^{-1}}{v-\kappa^{-1}} \right)^{\gamma} - \frac{\kappa}{v^n},$$ where $\gamma = 1+c_V^{-1}$. If the condition $n< \gamma$ is fulfilled (as is the case for $n= c_V= 1.5$), any expanding MS isentrope earlier or later crosses the spinodal because the specific entropy $s$, calculated from Eq.  (\[GW:s(v,tet)=\]) along the spinodal (\[GW:tet\_sp(v)=\]), monotonically increases with $v$. The latter provides one more argument in favor of GWEOS because the original van der Waals EOS with $c_V=1.5$ does not have this property. 

Two of the four isentropes in Fig.  \[f:2\] represent special cases.  The critical isentrope $s=s_{cr}\equiv s(1,1)$ crosses both the binodal and the spinodal at the critical point $CP$; the metastable segment $BD$ is absent, and further expansion proceeds along the EQ-GWEOS branch.  The special property of the isentrope $s= s_{\ast} \equiv s(v_{\ast},\theta_{\ast})$ is that it crosses the spinodal at the zero pressure $p_D =0$. The parameters of the crossing point (for $n=1.5$) are $$\label{GW:v_star=}
  v_{\ast} =\frac{n}{n+1} =0.6, \quad
  \theta_{\ast}= \frac{1}{4} \left(\frac{n+1}{n}\right)^{n+1} \approx 0.89652.$$ This state corresponds to the highest possible temperature of a sharp liquid surface in vacuum.  All isentropes with $s< s_{\ast}$ have an equilibrium metastable boundary with vacuum at point $E$, where $p_s(v_E) =0$ and $\kappa^{-1}< v_E <n/(n+1)$. For $s> s_{\ast}$ no stable or metastable boundary with vacuum at a non-zero density is possible. 

Phase-flip approximation \[s:pf\]
=================================

In a hydrodynamic flow, the thermodynamic state of each fluid element traces a certain trajectory in the $(v,p)$ parametric plane.  If such a trajectory crosses the binodal and enters the liquid-gas coexistence region, a non-trivial problem of the adequacy of the implemented EOS arises.  The fully equilibrium EQ EOS, which is based on the Maxwell construction and would be the obvious choice, is poorly justified for fast dynamic processes because it assumes an infinitely fast relaxation to the full equilibrium immediately below the binodal, and, as a consequence, does not allow negative pressures.  In practice, negative pressures do often play an important role, like in cavitation, or by providing a stabilizing mechanism for preservation of a sharp liquid-vacuum boundary up to relatively high surface temperatures — an important issue for certain laser-driven experiments. 

A better founded approach would be, once below the binodal, to proceed along the MS segment $BD$ of the thermodynamic trajectory and, simultaneously, keep track of the relaxation timescale $\tau_r$ from the MS to the EQ state.  The MS-decay time $\tau_r$ must be compared with the dynamic timescale $\tau_h$ of the specific experiment.  Clearly, so long as $\tau_r \gg \tau_h$, the MS branch should be a more adequate choice of EOS than the EQ one. 

The dynamic timescale $\tau_h$ is determined either by the rate of hydrodynamic expansion or by such external factors as the boundary conditions, the rate of external heating, duration of observation, etc.  By hydrodynamic expansion of laser-heated samples into vacuum the physically relevant timescale can be as short as $\tau_h \simeq \lambda/c_s \simeq 10^{-10}$ s, where $\lambda \simeq 1$ $\mu$m is the typical laser wavelength, and $c_s \simeq 10^6$ cm/s is the typical sound speed.  For nanosecond laser pulses it can amount to $\tau_h \sim 10^{-9}$–$10^{-7}$ s. Under quasi-static conditions much larger values of $\tau_h$ are possible. 

In superheated liquids the relaxation time $\tau_r$ can be evaluated on the basis of the theory of homogeneous bubble nucleation [@Skripov1974; @Blander_Katz1975; @Skripov_Skripov1979].  According to this theory, a superheated liquid state decays due to thermodynamic fluctuations, resulting in spontaneous creation of vapor bubbles with supercritical radii $r \geq r_c$ that subsequently grow in size.  The critical radius $r_c$ in a given state $(p,T)$ is defined by the condition $$\label{pf:r_c=}
  p_{sat}(T) = p+ \frac{2\sigma}{r_c},$$ where $p_{sat}(T)$ is the saturated vapor pressure, and $\sigma =\sigma(\rho,T)$ is the surface tension; here and throughout the text before Eq.  (\[pf:Hug=\]) the thermodynamic variables are in conventional units.  The energy barrier $W_c$ for the creation of a critical bubble is evaluated as $$\label{pf:W_c=}
  W_c = \frac{4\pi}{3} \sigma r_c^2 = \frac{16\pi\sigma^3}{3\left[p_{sat}(T)-p \right]^2}.$$ Then the MS${}\rightarrow{}$EQ transition time can be evaluated as [@Faik_Basko.2012] $$\label{pf:tau_r=}
  \tau_r = \nu_V \left(V_c J \right)^{-1}, \quad J = B\, \exp\left( -\frac{W_c}{T}\right),$$ where $0 \leq \nu_V= \nu_V(\rho,T) \leq 1$ is the volume fraction of vapor in the corresponding EQ state, $V_c =4\pi r_c^3/3$ is the volume of the critical bubble, $J$ \[cm${}^{-3}$ s${}^{-1}$\] is the spontaneous creation rate of bubbles with radius $r_c$, and $B = B(\rho,T)$ is a slowly varying preexponential factor. 

The key factor, determining $\tau_r$, is the Gibbs exponent $\exp\left(W_c/T\right)$, which, in particular, yields $\tau_r =\infty$ at the binodal.  As a thermodynamic trajectory penetrates the MS region and approaches the spinodal, $\tau_r$ decreases and falls below $\tau_h$ somewhere close to the spinodal.  Note that even in quasi-static laboratory experiments at $p \ll P_{cr}$ superheats up to $T {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}0.9 T_{cr}$ have been observed [@Skripov1974; @Blander_Katz1975]. 

Theoretical estimates for specific substances [@Martynyuk1977; @Skripov_Skripov1979; @Faik_Basko.2012] indicate that the transition from the condition $\tau_r \gg \tau_h$ to its opposite $\tau_r \ll \tau_h$ occurs near the spinodal in the course of a steep drop of $\tau_r$ by several orders of magnitude over a very narrow interval of temperatures $\Delta T/T_{cr} \simeq 1$–3% (for quasi-isobaric trajectories with $p \ll P_{cr}$).  In the immediate vicinity of the spinodal the model of homogeneous bubble nucleation breaks down [@Skripov1974]; at the spinodal itself the timescale $\tau_r$ can roughly be estimated as not larger than a few tens of intermolecular collision times, i.e. $\tau_r {\:\raisebox{.25ex}{$<$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}10^{-12}$–$10^{-11}$ s in superheated liquids — which is practically always shorter than the relevant hydrodynamic timescales $\tau_h$ in laser-heating experiments. 

A very sharp transition from the situation with $\tau_r \gg \tau_h$ to that of $\tau_r \ll \tau_h$ suggests that one can make a simplifying assumption of $\tau_r = \infty$ along the initial stretch of the MS-EOS trajectory below the binodal until the transition point is reached, beyond which $\tau_r=0$ can be assumed.  Earlier, a mathematical criterion for the moment of such transition was derived in Ref.   in the framework of the homogeneous nucleation theory.  Because it is assumed to occur instantaneously, this MS${}\rightarrow{}$EQ transition may be called a *phase flip (PF)*; in literature also the term *phase explosion* is used [@Martynyuk1977]. 

In this work we introduce a further simplification by assuming that the instantaneous and irreversible MS${}\rightarrow{}$EQ transition — the phase flip — occurs just at the intersection with the spinodal.  More precisely, if the thermodynamic trajectory of a fluid element enters the MS region through the liquid binodal branch, the phase flip occurs when crossing the liquid branch of the spinodal; respectively, if the MS region is entered from the vapor side, it occurs upon reaching the vapor spinodal branch.  We refer to this variant of hydrodynamic treatment of fluids with a liquid-gas phase transition as the *PF hydrodynamics*. 

![( Color online) Two examples of the phase-flip Hugoniot $hFGHJh_{PF}$ (red), calculated from Eq.  (\[pf:Hug=\]) for the initial state $D$ lying on the liquid (a) and the vapor (b) branches of the spinodal in GWEOS with $n=c_V=1.5$. In both cases the pre-flip state $D$ is presumed to be reached along an expansion isentrope $OBD$. On the upper graph (a) the isentrope $OBD$ is artificially tilted counterclockwise to avoid visual coalescence with the Hugoniot segment $h$ \[as on graph (b)\].  Dash-dotted curve $h_{sp}$ (orange) is the locus of isochoric post-flip states $G$ as the initial state $D$ slides along the spinodal $sp$ (magenta).  \[f:3\]](fig_3a_Hug-liq_v2.eps "fig:"){width="\myfigwidth"} ![( Color online) Two examples of the phase-flip Hugoniot $hFGHJh_{PF}$ (red), calculated from Eq.  (\[pf:Hug=\]) for the initial state $D$ lying on the liquid (a) and the vapor (b) branches of the spinodal in GWEOS with $n=c_V=1.5$. In both cases the pre-flip state $D$ is presumed to be reached along an expansion isentrope $OBD$. On the upper graph (a) the isentrope $OBD$ is artificially tilted counterclockwise to avoid visual coalescence with the Hugoniot segment $h$ \[as on graph (b)\].  Dash-dotted curve $h_{sp}$ (orange) is the locus of isochoric post-flip states $G$ as the initial state $D$ slides along the spinodal $sp$ (magenta).  \[f:3\]](fig_3b_Hug-gas_v2.eps "fig:"){width="\myfigwidth"}

Generally, the dynamic response of a fluid to a jump in its thermodynamic properties would lead to a flow discontinuity, i.e. generally a phase flip should occur inside a PF shock front.  In principle, it may be either a compression or a rarefaction PF shock.  Conservation of mass, momentum and energy implies that, as in the theory of conventional shocks, the fluid states in front (initial) and behind (final) a PF shock must obey the Hugoniot relation [@Courant_Friedrichs1977; @LL-H87] $$\label{pf:Hug=}
  e_2 -e_1 = \frac{1}{2} (p_1+p_2) (v_1-v_2).$$ The initial state (point $D$ in Fig.  \[f:3\]) is fully defined by its specific volume $v_1$: because it lies on the spinodal, we have $\theta_1= \theta_{sp}(v_1)$, while $p_1 = p(v_1,\theta_1)$, $e_1= e(v_1,\theta_1)$ and $s_1= s(v_1,\theta_1)$ are given by Eqs.  (\[GW:p(v,tet)=\])–(\[GW:s(v,tet)=\]).  The final state $p_2 = p_{EQ}(v_2,\theta_2)$, $e_2 =e_{EQ}(v_2,\theta_2)$ and $s_2 =s_{EQ}(v_2,\theta_2)$ must be calculated from the EQ EOS by invoking the Maxwell construction.  Once the initial state $v_1$ is fixed and $v_2$ is varied, Eq.  (\[pf:Hug=\]) defines a unique PF Hugoniot curve $p_2(v_2)$. Two characteristic examples of the PF Hugoniot with initial states on the liquid and the vapor spinodal branches are shown in Figs.  \[f:3\]a and \[f:3\]b. 

Because the final state is calculated from a different EOS branch, the PF Hugoniot bears a close resemblance to the well-known detonation-deflagration Hugoniot curve [@Courant_Friedrichs1977; @LL-H87].  Similar to the latter, the segment $GH$ between the vertical and horizontal lines, passing through the initial state $D$, has no physical meaning because the square of the mass flux $$\label{pf:j^2=}
  j^2 =(\rho_1\eta_1)^2 =(\rho_2\eta_2)^2 =\frac{p_2-p_1}{v_1-v_2}$$ across a PF front would be negative there; here $\eta_1$ and $\eta_2$ are, respectively, the flow velocities before and after the PF front in the reference frame comoving with this front.  For the problem of adiabatic unloading into vacuum, studied in this work, only the segment below the $H$ point is relevant. 

The phase flip at a constant density from $D$ to $G$ can be singled out as a special case: it occurs, for example, by a quasistatic volumetric heating (cooling) of an extended uniform liquid (vapor) volume at constant pressure [@Faik_Basko.2012].  In this case, when $v_1=v_2$, Eq.  (\[pf:Hug=\]) yields a “PF image” $p_2(v_1)$ of the spinodal — the curve $h_{sp}$ in Fig.  \[f:3\], which is the locus of final states $G$ when the initial state $D$ slides along the spinodal.  Note that the $h_{sp}$ curve lies everywhere between the spinodal and the binodal \[in both the $(v,p)$ and $(v,\theta)$ parametric planes\], which implies that the isochoric explosive boiling of superheated liquid, as well as the rapid condensation of supersaturated vapor, are both accompanied by an upward jump in pressure and temperature. 

In summary, the EQ and PF approximations represent two opposite extremes in the hydrodynamic description of fluid dynamics in the region of liquid-gas phase transition.  Once both the EQ and the PF solutions are found for one and the same problem, a comparison between the two can provide valuable information about a possible impact of metastable states.  Because the EQ option ignores metastable states altogether, whereas the PF approximation allows maximum penetration into the MS region, one can expect that such a comparison should provide a good measure for maximum possible dynamical effects of metastability. 

Centered rarefaction wave by expansion into vacuum \[s:crw\]
============================================================

Self-similarity ansatz
----------------------

The unloading of an initially uniform half-space bordering on vacuum is described by the one-dimensional self-similar solution for a centered rarefaction wave (CRW), where one family of characteristics in the $(x,t)$ plane is a set of straight lines originating from a single center [@Courant_Friedrichs1977; @LL-H87; @Zeldovich_Raizer2012].  A remarkable property of CRW is that it remains self-similar for any equation of state.  If the rarefaction flow is continuous, it is also isentropic.  If the initial condition is $\rho(t,x) =\rho_0$, $p(t,x) =p_0$, $\theta(t,x) =\theta_0$ for $x<0$ at $t=0$, than all the flow parameters are functions of a single independent variable $\xi=x/t$, and the solution takes the form $$\begin{aligned}
\label{crw:xi=u-c}
  \xi &\equiv&  \frac{x}{t} = u(\rho)-c_s(\rho), \\ \label{crw:u=}
  u(\rho) &=& \int\limits_{\rho}^{\rho_0} \frac{dp(\rho')}{\rho' c_s(\rho')} = \int\limits_{\rho}^{\rho_0} c_s(\rho')\,d\ln\rho',\end{aligned}$$ where $u$ is the flow velocity, and $c_s$ is the isentropic sound velocity.  In Eqs.  (\[crw:xi=u-c\]) and (\[crw:u=\]) both $p(\rho)$ and $c_s(\rho)$ are considered to be known functions of $\rho$ along the expansion isentrope $s=s_0$. Once the integral in Eq.  (\[crw:u=\]) is calculated for the given EOS, Eq.  (\[crw:xi=u-c\]) defines in implicit form the functions $\rho(t,x)$ and $u(t,x)$.

For subsequent analysis it will be convenient to introduce the notion of the *phase* of CRW as an arbitrary point on its profile with some fixed values of the material velocity $u$, density $\rho \equiv v^{-1}$, and other thermodynamic variables.  Equation (\[crw:xi=u-c\]) tells us that the phase is uniquely defined by the value of $\xi$, which, in its turn, can be called the *phase velocity* because points with fixed $\rho$ and $p$ propagate in space with the velocity $\xi$. The head of the rarefaction wave, where $\rho= \rho_0$ and $u=0$, propagates into the unperturbed matter with the phase velocity $$\label{crw:xi_0=}
  \xi_0 =-c_s(\rho_0) \equiv -c_{s0}.$$ The difference $\eta \equiv u-\xi$ is the material velocity with respect to a fixed wave phase, i.e. the *phase-relative* flow velocity.  One can prove that the absolute magnitude of the phase-relative velocity in any continuous adiabatic flow is equal to the local speed of sound $c_s$, as is expressed by Eq.  (\[crw:xi=u-c\]) for the CRW case. 

The EQ solution
---------------

The CRW solution for the EQ EOS has been described in detail in Ref.  .  It is continuous in space, and the thermodynamic trajectories of all fluid elements follow one and the same expansion EQ isentrope $s=s_0 \equiv s_{EQ}(v_0,\theta_0)$ defined by the initial state $O$, as is shown in Fig.  \[f:2\].  Qualitatively, the spatial CRW profiles look similar for all unloading EQ isentropes.  As an illustration, Fig.  \[f:4\] shows the EQ density profiles $\rho(\xi)$ (black dash-dotted) for two different $s_0$ values.  More data on the EQ-CRW solutions for five different values of $s_0$ are given in Table \[t:1\]. 

  ----------------------------- ----------------- ----------------- ----------------------------------- --------------- ---------------
  Case                                  1                 2                          3                         4               5
                                 $s_0< s_{\ast}$   $s_0= s_{\ast}$   $s_{\ast}\!  <\!  s_0\!  <\!  s_{cr}$   $s_0= s_{cr}$   $s_0> s_{cr}$
  \[4pt\]                                                                                                               
  \[-4pt\] $v_B$                     0.27468           0.56574                    0.65070                      1            2.06392
  $p_B$                              0.01055           0.73952                    0.85803                      1            0.79175
  $\theta_B$                         0.54055           0.95166                    0.97499                      1            0.96226
  $c_{sB_+}$                         6.64741           2.87141                    2.63834                   2.23607         2.05383
  $c_{sB_-}$                         0.00925           0.49076                    0.60481                   0.91158         1.27836
  $\Gamma_{B_-}$                     0.03677           0.6091                     0.6809                    0.6909          0.9617
  \[3pt\] $v_D$                       0.321              0.6                      0.67482                      1              2.5
  $p_D$                              -11.948              0                         0.5                        1            0.48067
  $\theta_D$                         0.39186           0.89652                    0.94169                      1            0.83642
  $c_{sD_+}$                         2.97070           2.54066                    2.46711                   2.23607         1.77828
  \[3pt\] $\Delta\xi_{sh,EQ}$        6.63816           2.38065                    2.03353                   1.32448         0.77547
  $\Delta\xi_{lvt,EQ}$               0.00034           0.2989                     0.4118                    0.6298          1.2294
  ----------------------------- ----------------- ----------------- ----------------------------------- --------------- ---------------

  : \[t:1\] Parameters of the EQ-CRW solutions for five representative cases with different values of the initial entropy $s_0$, calculated for GWEOS with $n=c_V=1.5$. For notation see the text. 

A salient feature of the EQ-CRW profiles is the *binodal shelf* — a uniform layer $B_+ B_-$ with constant values of $u$, $\rho$ and other thermodynamic variables.  The binodal shelf is formed at point $B$ (see Fig.  \[f:2\]), where the unloading isentrope crosses the binodal at $\rho= \rho_B \equiv v_B^{-1}$ and the isentropic speed of sound $c_s$ experiences a jump from a higher upstream value $c_{sB_+}$, calculated along the MS-EOS branch, to a lower downstream value $c_{sB_-}$, calculated along the EQ-EOS branch.  The flow velocity $u_B=u(\rho_B)$ on the shelf is given by Eq.  (\[crw:u=\]). 

Qualitatively, the formation of a binodal shelf can be explained as follows: while material enters the state $\rho=\rho_B$ from the upstream side with a phase-relative velocity $\eta_{B_+}= c_{sB_+}$, it leaves this state with a lower phase-relative velocity $\eta_{B_-}= c_{sB_-}$; hence, material accumulates in state $B$ at a rate $\rho_B(c_{sB_+}- c_{sB_-})> 0$, and the width of the resulting uniform layer grows in time at a speed $$\label{crw:Dxi_B=}
  \Delta \xi_{sh,EQ} \equiv \xi_{B_-}-\xi_{B_+} = c_{sB_+}- c_{sB_-} >0.$$ That is to say, a single point $B$ on the $(v,p)$ plane is represented by a finite interval $\xi_{B_-}\leq \xi\leq \xi_{B_+}$ of the phase velocities. 

At the boundary with vacuum, the phase velocity $\xi_e =u_e$ coincides with the flow velocity $u_e=u(0)$ calculated from Eq.  (\[crw:u=\]).  Normally, the integral in (\[crw:u=\]) converges in the limit of $\rho \to 0$, though relatively slowly.  The layer $\xi_{B_-} <\xi< \xi_e$ downstream from the binodal shelf is composed of a disperse mixture of liquid and vapor.  At low entropies $s_0 {\:\raisebox{.25ex}{$<$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}s_{\ast}$, this *liquid-vapor tail* exhibits a steep density drop, clearly manifested by Fig.  \[f:4\]a.  From Eqs.  (\[crw:xi=u-c\]) and (\[crw:u=\]) one calculates the effective width of this tail to be $$\label{crw:Dxi_lvt=}
  \Delta \xi_{lvt,EQ} \equiv\left| \frac{d\ln\rho}{d\xi}
  \right|^{-1}_{\xi =\xi_{B_-}+0} = c_{sB_-}\Gamma_{B_-},$$ where $\Gamma_{B_-}$ is the value of the fundamental gasdynamic derivative [@Thomson1971] $$\label{crw:Gam=}
  \Gamma =1+\left(\frac{\partial \ln c_s}{\partial\ln\rho}\right)_s,$$ calculated along the EQ branch of the expansion isentrope immediately below point $B$. In the present EOS, the tail width $\Delta \xi_{lvt,EQ}$ is a monotonically increasing function of the specific volume $v_B$ (see Table \[t:1\]) over a broad interval $0.2< v_B {\:\raisebox{.25ex}{$<$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}40$. As $p_B\to 0$ ($v_B \to 0.2$) on the liquid branch, $\Delta \xi_{lvt,EQ}$ very rapidly approaches zero, and the tail becomes physically insignificant — which is expected to be true for any realistic EOS. 

![( Color online) EQ (dash-dotted black) and PF (solid red) density profiles in the self-similar centered rarefaction wave, calculated for cases 2 (a) and 3 (b) in Table \[t:1\].  \[f:4\]](fig_4a_sss.eps "fig:"){width="0.5\myfigwidth"} ![( Color online) EQ (dash-dotted black) and PF (solid red) density profiles in the self-similar centered rarefaction wave, calculated for cases 2 (a) and 3 (b) in Table \[t:1\].  \[f:4\]](fig_4b_sss.eps "fig:"){width="0.5\myfigwidth"}

The PF solution
---------------

Depending on how the initial entropy $s_0$ compares with $s_{\ast}$, the PF-CRW solution can be of one of the two types, illustrated by Figs.  \[f:4\]a and \[f:4\]b respectively.  When $s_0 \leq s_{\ast}$, the boundary condition $p(\rho_E) =0$ is satisfied at point $E$ (see Fig.  \[f:2\]) with a non-zero density $\rho_E \equiv v_E^{-1} \geq 1+n^{-1}$ before the MS isentrope reaches the spinodal.  In this case, displayed by Fig.  \[f:4\]a, no phase flip occurs, the liquid-vapor tail is absent, and the boundary with vacuum remains a sharp surface of superheated liquid.  Because material enters the boundary state $E$ from the upstream side with the phase-relative velocity $\eta_{E_+} =c_{sE}$, it accumulates there with the rate $\rho_E c_{sE}$ and forms a shelf $E_+E_-$ of width $$\label{crw:Dxi_E=}
  \Delta\xi_{sh,PF} \equiv \xi_{E_-}-\xi_{E_+} =c_{sE}.$$ That is, for $-c_{s0}\leq \xi \leq \xi_{E_+}$ the PF solution is given by Eqs.  (\[crw:xi=u-c\]) and (\[crw:u=\]), while at $\xi_{E_+}< \xi< \xi_{E_-}$ the flow is uniform with $\rho=\rho_E$ and $u=u_e=u(\rho_E)$; the phase velocities at the two ends of the $E_+E_-$ shelf are $\xi_{E_+}= u_e-c_{sE}$ and $\xi_{E_-}\equiv \xi_e= u_e$. Parameters of state $E$ for the illustrative cases 1 and 2 (in case 2 points $E$ and $D$ coincide) are listed in Table \[t:2\]. 

  ----------------------------- ----------------- ----------------- --------------------------- ---------------
  Case                                  1                 2                      3                     5
                                 $s_0< s_{\ast}$   $s_0= s_{\ast}$   $s_{\ast} < s_0 < s_{cr}$   $s_0> s_{cr}$
  \[4pt\]                                                                                       
  \[-4pt\] $v_E$                     0.27470             0.6                     –                     –
  $\theta_E$                         0.54046           0.89652                   –                     –
  $c_{sE}$                           6.64519           2.54066                   –                     –
  \[3pt\] $v_J$                         –                 –                   2.57976               5.20258
  $p_J$                                 –                 –                   0.30343               0.30995
  $\theta_J$                            –                 –                   0.82902               0.83164
  $c_{sJ}$                              –                 –                   0.82871               1.30758
  $\Gamma_{J}$                          –                 –                   0.8711                1.0074
  $s_J-s_0$                             –                 –                   0.17477               0.18982
  $\nu_{m,J}$                           –                 –                   0.31055               0.70077
  \[3pt\] $\Delta\xi_{sh,PF}$        6.64519           2.54066                2.25034               1.14995
  $\Delta\xi_{lvt,PF}$                  0                 0                   0.7219                1.3173
  $\delta u_{sh}$                   0.000436           0.15895                0.09289               0.36660
  ----------------------------- ----------------- ----------------- --------------------------- ---------------

  : \[t:2\] Parameters of the PF-CRW solutions for four representative cases from Table \[t:1\].  For notation see the text. 

At higher initial entropies $s_0> s_{\ast}$, the unloading isentrope crosses the spinodal at point $D$ with $p_D>0$, and the vacuum boundary condition can only be satisfied after relaxation to the EQ EOS.  Because the assumed instantaneous MS${}\rightarrow{}$EQ transition incurs a jump in pressure and/or density, it must take place inside a flow discontinuity, i.e. inside a PF shock front with the pre-shock state $(v_1,p_1) = (v_D,p_D)$ at point $D$, and the post-shock state $(v_2,p_2)$ on the corresponding PF Hugoniot displayed in Fig.  \[f:3\].  In our case it will be a rarefaction shock with $p_2< p_D$ because the CRW evolves from a step-like initial pressure drop where all fluid elements are accelerated in one direction towards vacuum, i.e. the post-shock state $(v_2,p_2)$ must lie on the “deflagration” branch $HJh_{PF}$ of the PF Hugoniot. 

Further on, we argue that the post-shock state must be exactly at the Chapman-Jouguet point $J$, where the ray $DJ$ is tangent to the PF Hugoniot, the outflow velocity $\eta_2$ from the shock front equals the local speed of sound $c_{sJ}$, and the mass flux across the front $$\label{crw:j_J=}
  j\equiv \rho_2\eta_2 =j_{DJ} = \left( \frac{v_J-v_D}{p_D-p_J} \right)^{1/2} =\rho_J c_{sJ}$$ is maximum for the given initial state $D$ [@Courant_Friedrichs1977; @LL-H87].  Indeed, the conservation of mass and momentum dictates that, for any given $j$, the transition from the MS state $D$ to the EQ state 2 inside the shock front must proceed along the Rayleigh line $$\label{crw:p_tran(v)=}
  p(v) =p_D-(v-v_D)/j^2.$$ If $j< j_{DJ}$, the line (\[crw:p\_tran(v)=\]) crosses the PF Hugoniot at two points, say at $J'$ with $v_{J'}<v _J$ and at $J''$ with $v_{J''}>v _J$. Physically meaningful would only be the first solution $J'$ because the kinetics of the irreversible MS${}\rightarrow{}$EQ transition would ensure that the final EQ state be already reached there.  On the other hand, because the shock-outflow velocity $\eta_2$ would in this case be subsonic, $\eta_2< c_{sJ'}$, the zero pressure at the outer boundary would inevitably “suck out” the fluid from the shock front at $J'$ until the post-shock state $J$ with exactly the sonic outflow were reached; more detailed arguments to this point can be found in §20 of chapter XI in Ref.  .  Parameters of the Chapman-Jouguet point $J$ for our illustrative cases 3 and 5 are listed in Table \[t:2\].  Case 4 with $s_0=s_{cr}$ is a degenerate one in the sense that points $B$, $D$ and $J$ coalesce with the critical point $CP$, and the PF and EQ solutions coincide.  Note that the post-shock entropy $s_2(v_2)$ on the deflagration branch of the PF Hugoniot has its local maximum $s_J= s_2(v_J)$ at point $J$ [@Courant_Friedrichs1977]. 

Having uniquely determined the post-shock state $J$, we notice that matter enters state $D$ at the rate $\rho_D c_{sD_+}$, and leaves it through the shock front at the rate $j_{DJ}= \rho_J c_{sJ}$. Normally (see Fig.  \[f:3\]) the Rayleigh line $DJ$ slopes down less steeply than the MS isentrope at point $D$, which means that $\rho_D c_{sD_+}> j_{DJ}$ and material accumulates in state $D$ upstream from the PF shock.  In result, the *spinodal shelf* $D_+D_-$ of a uniform flow with $\rho=\rho_D$ is formed, whose width grows in time at a speed $$\label{crw:Delxi_sh,PF=}
  \Delta\xi_{sh,PF} \equiv \xi_{D_-}- \xi_{D_+} = c_{sD_+} -\frac{v_D}{v_J} c_{sJ}.$$

The overall density profile of the PF solution for $s_0> s_{\ast}$ is displayed in Fig.  \[f:4\]b.  Similar to the EQ solution, the head of the rarefaction wave at $-c_{s0} \leq \xi \leq \xi_{D_+}$ is described by Eqs.  (\[crw:xi=u-c\]) and (\[crw:u=\]), where the sound speed $c_s(\rho)$ is calculated along the initial MS isentrope $s=s_0$. The flow velocity on the spinodal shelf at $\xi_{D_+}< \xi < \xi_{D_-}$ is $u_D=u(\rho_D)$. At $\xi =\xi_{D_-} =\xi_J$ the fluid passes through a rarefaction jump from state $D$ to state $J$, in which the entropy increases from $s_0$ to $s_J$, the density drops from $\rho_D=v_D^{-1}$ to $\rho_J=v_J^{-1}$, and the flow velocity jumps from $u_D$ to $u_D +c_{sJ}(1-v_D/v_J)$. When $s_{\ast}< s_0< s_{cr}$, the pure liquid, entering the PF shock, undergoes instantaneous partial vaporization (explosive boiling) inside the shock front; for $s_0> s_{cr}$ the pure (dry) vapor, entering the PF shock, undergoes instantaneous partial condensation.  The values $\nu_{m,J}$ of the mass fraction of vapor in state $J$ behind the shock are listed in Table \[t:2\]. 

Downstream from the PF shock we have a liquid-vapor tail at $\xi_J< \xi <\xi_e=u_e$, analogous to that in the EQ solution.  There the flow is again continuous and isentropic, and the density and velocity profiles are found from $$\begin{aligned}
\label{crw:xi_dtPQ=u-c}
  \xi &\equiv&  \frac{x}{t} = u(\rho)-c_s(\rho), \\ \label{crw:u_dtPQ=}
  u(\rho) &=& u_D +c_{sJ}\left(1-\frac{v_D}{v_J}\right) + \int\limits_{\rho}^{\rho_J} c_s(\rho')\,d\ln\rho',\end{aligned}$$ where $c_s(\rho)$ at $\rho\leq \rho_J$ must be calculated along the EQ expansion isentrope $s=s_J$; the outer edge velocity $u_e=u(0)$. In Fig.  \[f:3\] the isentrope $s=s_J$ (dash-dotted black curve) visually merges with the corresponding segment $h_{PF}$ of the PF Hugoniot; note that the entropy increment $s_J-s_0$ in the PF shock (see Table \[t:2\]) is relatively small.  By analogy with Eq.  (\[crw:Dxi\_lvt=\]), the effective width of the liquid-vapor tail is given by $$\label{crw:Dxi_lvtPF=}
  \Delta\xi_{lvt,PF} = c_{sJ}\Gamma_J,$$ where $\Gamma_J$ is the fundamental gasdynamic derivative (\[crw:Gam=\]) at point $J$. The values of $\Gamma_J$ and $\Delta\xi_{lvt,PF}$ for the illustrative cases 3 and 5 are listed in Table \[t:2\].  It should be emphasized that for all $s_0> s_{\ast}$ the PF-CRW profiles look qualitatively similar irrespective of how $s_0$ compares with $s_{cr}$, i.e. irrespective of whether the PF shock is of the vaporization or condensation type. 

Comparison between the EQ and PF solutions
------------------------------------------

When comparing the EQ and PF profiles in CRW, we focus our attention on two main aspects that might be important for the interpretation of the experimental data.  The feature that both solutions have in common is a horizontal shelf on spatial profiles, representing a layer of uniform flow.  Only quantitative characteristics of the EQ and PF shelves are different.  The main qualitative difference is that, except for the degenerate case of $s_0=s_{cr}$, the PF flow contains a discontinuity: it is either a density jump at the outer boundary with vacuum or a rarefaction shock downstream from the uniform layer. 

When $s_0<s_{cr}$ and the expansion isentrope crosses the liquid binodal branch, the quantitative difference between the PF and EQ shelves is not large: it disappears in both limits of $p_B \to 0$ and $p_B \to 1$, and becomes maximum somewhere at $s_{\ast} {\:\raisebox{.25ex}{$<$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}s_0 < s_{cr}$ (at $0.7 {\:\raisebox{.25ex}{$<$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}p_B <1$ for our EOS).  Generally, the PF shelf is broader, $\Delta\xi_{sh,PF} > \Delta\xi_{sh,EQ}$, but has a somewhat lower density ($\rho_E$ or $\rho_D$) than that ($\rho_B$) in the EQ case.  Over the whole liquid branch, the values of both these parameters (cf.  Tables \[t:1\] and \[t:2\]) differ by no more than a few percent.  Also, the difference $$\label{crw:du_sh=}
  \delta u_{sh} = \int\limits^{\rho_B}_{\max\{\rho_E,\rho_D\}} c_s(\rho')\,d\ln\rho'$$ between the flow velocities on the PF and EQ shelves, listed in the last row of Table \[t:2\], remains small for $s_0< s_{cr}$. On the one hand, these results confirm the statement by Inogamov *et al.* [ @Inogamov_Petrov.1999] that penetration into the metastable region of a superheated liquid should have only a minor effect on the unloading profiles.  On the other hand, the same comparison for supercritical isentropes $s_0>s_{cr}$ (case 5 in Tables \[t:1\] and \[t:2\]) reveals that rapid condensation of the supersaturated vapor leads to a considerable distortion of the CRW profiles, which becomes more pronounced with the increasing entropy $s_0$.

The new quality, brought in by the PF approximation, is the presence of a discontinuity, which formally disappears only for the critical isentrope $s_0=s_{cr}$. Important practical implications of this discontinuity may be not so in the mechanical as in the optical properties of the rarefaction flows because in laser-driven experiments the absorption of the incoming laser light is sensitive to the density gradient in the expanding material.  But even in this respect, the PF solution may still be practically indistinguishable from the EQ one if the effective width $\Delta\xi_{lvt,EQ}$ of the EQ-CRW tail, defined in Eq.  (\[crw:Dxi\_lvt=\]), is small enough — which is definitely the case in the limit of $p_B \to 0$ along the liquid branch of the binodal. 

Firstly, we note that, from the viewpoint of laser absorption, the PF density profiles for practically all subcritical isentropes $s_0<s_{cr}$ can be considered as ending with a sharp liquid surface of superheated liquid bordering on vacuum.  Indeed, for $s_0< s_{\ast}$ ($p_B <0.74$, $\theta_B< 0.95$) this is literally the case, while for $s_0>s_{\ast}$, even as close to the critical point as $p_B = 0.86$, $\theta_B= 0.975$ (case 3 in Table \[t:1\]), we still observe a density jump by a large factor of $v_J/v_D =3.82$. The PF density jump disappears only in the immediate vicinity of the critical point, where the PF approximation becomes inaccurate anyway.  Hence, the impact of metastability on the laser absorption can be claimed significant whenever the density gradient in the EQ liquid-vapor tail becomes of significance in this regard.  This happens when the effective width of this tail, which in conventional units is given by $$\label{crw:h_lvt,EQ=}
  h_{lvt,EQ} =\Delta\xi_{lvt,EQ} \left(P_{cr}/\rho_{cr}\right)^{1/2} t,$$ becomes comparable to the laser wavelength $\lambda$ \[actually, $h_{lvt,EQ} {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}(0.1$–0.2)${}\lambda$ would suffice\].  Simple estimates indicate that for multi-nanosecond laser pulses this can often be the case for $\Delta\xi_{lvt,EQ} {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}0.03$, which corresponds to the boiling pressures $p_B {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}0.2$ within our EOS. 

For supercritical isentropes $s_0>s_{cr}$ the difference between the PF and EQ density profiles is significant anyway, and an extra density jump, whose amplitude increases with $s_0$, only reinforces this conclusion. 

Conclusions
===========

The approximation of the phase-flip (PF) hydrodynamics, based on the assumption of instantaneous decay of the metastable liquid/vapor states at the spinodal, provides a possibility to explore the effects of metastability on the dynamic behavior of fluids with a liquid-gas phase transition without invoking more complex kinetic models.  In this work, the first exact solution of the PF hydrodynamics is derived and investigated for the centered rarefaction wave (CRW) by unloading of a uniform layer into vacuum.  Numerical results are obtained for a particular version of the generalized van der Waals EOS in the reduced form. 

Of special interest for practical problems becomes a comparison of the PF solution with the simpler and more familiar EQ solution [@Sokol-Tint_Bialkowski.1998; @Inogamov_Anisimov.1999; @Zhao_Mentrelli.2011], obtained for the same problem by using the fully equilibrium version of EOS in the two-phase region.  One would expect that such a comparison should reveal the maximum possible influence of metastability in the problem considered.  Here, having applied this approach to the classical problem of unloading into vacuum, we obtained the following results. 

We find that the new qualitative feature, emerging in the PF solution, is the appearance of a flow discontinuity, i.e. of a rarefaction shock front.  Nevertheless, so long as the boiling pressure (normalized by $P_{cr}$) on a liquid expansion isentrope stays low $p_B \ll 1$, the difference between the continuous EQ and the discontinuous PF solutions remains, to all purposes, insignificant because a very steep drop in the EQ density profile downstream from the binodal shelf can hardly be distinguished from the step function in the PF profile.  But as the boiling pressure approaches the critical value, metastability becomes an issue from the point of view of absorption and reflection of the laser light.  In particular, for $s_0\leq s_{\ast}$ ($p_B\leq 0.74$ in our EOS) the PF solution predicts a sharp liquid surface at the boundary with vacuum, which for metals can be highly reflective according to the Fresnel formulae.  The EQ solution, on the contrary, exhibits a smooth density slope along the disperse liquid-vapor tail, with the steepest gradient characterized by the length scale $h_{lvt,EQ}$ given by Eq.  (\[crw:h\_lvt,EQ=\]).  With $\Delta\xi_{lvt,EQ} \approx 0.3$ for $s_0=s_{\ast}$ (case 2 in Table \[t:1\]) and $\left(P_{cr}/\rho_{cr}\right)^{1/2} {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}3\times 10^4$ cm/s for typical metals, we calculate $h_{lvt,EQ} {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}0.1$ $\mu$m for times $t {\:\raisebox{.25ex}{$>$}\hspace*{-.75em}
\raisebox{-.93ex}{$\sim$}\:}1$ ns, which would be comparable to the laser wavelength $\lambda$ in most experiments with nanosecond laser pulses.  As a consequence, simulations of laser experiments, based on the EQ EOS, may significantly overestimate the laser absorption at low irradiation intensities near the ablation threshold.  In particular, this could explain the discrepancy between the observed and predicted onset of the ablation regime in Ref.  . 

At the same time, there are many experiments with ultrashort laser pulses where the laser pulse terminates before the target begins to expand, and the key parameters for adequate interpretation of the experimental data become the width, mass and velocity of the plateau layer (the binodal/spinodal shelf) in the rarefaction wave[@Sokol-Tint_Bialkowski.1998; @Basko_Krivokor.2017].  In this regard our analysis confirms and quantifies the earlier conclusion [@Inogamov_Petrov.1999] about the insignificant role of metastability, but under the condition that the expansion isentrope remains undercritical, $s_0< s_{cr}$. For supercritical isentropes $s_0> s_{cr}$, where the rarefaction shock is a condensation jump, the PF approximation practically always predicts significant deviations from the EQ unloading profiles.  This fact should be taken into account by interpretation of certain types of experiments where the EQ EOS is used. 

The author is grateful to Anna Tauschwitz for many stimulating discussions.  This work was supported by the Russian Science Foundation through grant No.  14-11-00699-$\Pi$.
