Nonlinear dispersive regularization of inviscid gas dynamics

Ideal gas dynamics can develop shock-like singularities with discontinuous density. Viscosity typically regularizes such singularities and leads to a shock structure. On the other hand, in 1d, singularities in the Hopf equation can be non-dissipatively smoothed via KdV dispersion. Here, we develop a minimal conservative regularization of 3d ideal adiabatic flow of a gas with polytropic exponent $\gamma$. It is achieved by augmenting the Hamiltonian by a capillarity energy $\beta(\rho) (\nabla \rho)^2$. The simplest capillarity coefficient leading to local conservation laws for mass, momentum, energy and entropy using the standard Poisson brackets is $\beta(\rho) = \beta_*/\rho$ for constant $\beta_*$. This leads to a Korteweg-like stress and nonlinear terms in the momentum equation with third derivatives of $\rho$, which are related to the Bohm potential and Gross quantum pressure. Just like KdV, our equations admit sound waves with a leading cubic dispersion relation, solitary and periodic traveling waves. As with KdV, there are no steady continuous shock-like solutions satisfying the Rankine-Hugoniot conditions. Nevertheless, in 1d, for $\gamma = 2$, numerical solutions show that the gradient catastrophe is averted through the formation of pairs of solitary waves which can display approximate phase-shift scattering. Numerics also indicate recurrent behavior in periodic domains. These observations are related to an equivalence between our regularized equations (in the special case of constant specific entropy potential flow in any dimension) and the defocussing nonlinear Schrodinger equation (cubically nonlinear for $\gamma = 2$), with $\beta_*$ playing the role of $\hbar^2$. Thus, our regularization of gas dynamics may be viewed as a generalization of both the single field KdV&NLS equations to include the adiabatic dynamics of density, velocity, pressure&entropy in any dimension.


I. INTRODUCTION
Gas dynamics has been an active area of research with applications to high-speed flows, aerodynamics and astrophysics. The equations of ideal compressible flow are known to encounter shock-like singularities with discontinuities in density, pressure or velocity 1 . These singularities are often resolved by the inclusion of viscosity. However, as the KdV equation (u t + uu x = εu xxx ) illustrates, such singularities in the one-dimensional (1d) Hopf (or kinematic wave) equation u t + uu x = 0 can also be regularized conservatively via dispersion 2 , as in dispersive shock wave theory (see 1,[3][4][5] and references therein) with applications to undular bores in shallow water and blast waves in Bose-Einstein condensates. In this paper, we develop a minimal conservative regularization of ideal gas dynamics, which we refer to as R-gas dynamics. Somewhat analogous conservative 'rheological' regularizations of vortical singularities in ideal Eulerian hydrodynamics, magnetohydrodynamics and two-fluid plasmas have been developed in [6][7][8] . The current work may be regarded as a way of extending the single-field KdV equation to include the dynamics of density, velocity and pressure and also to dimensions higher than one. There is of course a well-known generalization of KdV to 2d, the Kadomtsev-Petviashvili (KP) equation 9 . However, unlike KP, our regularized equations are rotation-invariant and valid in any dimension. Now, recall 10 that the dispersive regularization term in the KdV equation u t − 6uu x + u 3x = 0 arises from the gradient energy term in the Hamiltonian H = (u 3 + (1/2)u 2 x ) dx, upon use of the Poisson brackets {u(x), u(y)} = ∂ x δ (x − y). In fact, KdV does not conserve mechanical and capillarity energies separately 11,12 . By analogy with this, we obtain our regularized model by augmenting the Hamiltonian of ideal adiabatic flow of a gas with polytropic exponent γ, by a density gradient energy β (ρ)(∇ρ) 2 . Such a term arose in the work of van der Waals and Korteweg [13][14][15][16] in the context of capillarity, but can be important even away from interfaces in any region of rapid density variation, especially when dissipative effects are small, such as in weak shocks, cold atomic gases, superfluids and collisionless plasmas. It has also been used to model liquid-vapor phase transitions and in the thermomechanics of interstitial working 15 . We argue that the simplest choice of capillarity coefficient that leads (using the standard Poisson brackets) to local conservation laws for mass, momentum, energy and entropy (with the standard mass, momentum and entropy densities) is β (ρ) = β * /ρ where β * is a constant. By contrast, the apparently simpler option of taking β (ρ) constant leads, in 1d, to a KdV-like linear dispersive term ρ xxx in the velocity equation, but results in a momentum equation that, unlike KdV 11 , is not in conservation form for the standard momentum density ρu. A consequence of the constitutive law β = β * /ρ is that the ideal momentum flux ρu 2 + p is augmented by a stress −β * (ρ xx − ρ 2 x /ρ) corresponding to a Kortweg-type grade 3 elastic material 15,16 . This leads to new nonlinear terms in the momentum equation with third derivatives of ρ, somewhat reminiscent of KdV. One of the effects of these nonlinear dispersive terms is to allow for 'upstream influence' 17 which is forbidden by the hyperbolic equations of inviscid gas dynamics under supersonic conditions. Interestingly, our regularization term is also related to the quantum mechanical Bohm potential 18 and Gross quantum pressure (p.476 of 19 ) encountered in superfluids. Moreover, unlike KdV, our equations extend in a natural way to any dimension. Remarkably, for potential flow (v = ∇φ ) in the isentropic case (globally constant entropy and p ∝ ρ γ ), the R-gas dynamic equations may be transformed into the nonlinear Schrödinger equation (NLSE) via the Madelung transformation 20 ψ = √ ρ exp iφ /2 β * with β * playing the role ofh 2 . This equivalence, which may be regarded as a conservative analog of the Cole-Hopf transformation for Burgers, applies in any dimension, and results in a defocusing NLSE with |ψ| 2(γ−1) ψ nonlinearity, so that one obtains the celebrated cubic NLSE for γ = 2. The lat-ter is known to admit an infinite number of conservation laws and display recurrence. It is noteworthy that the quantum version of the 1d cubic NLSE (Lieb-Liniger model) has recently been given a hydrodynamical description (generalized hydrodynamics 21,22 ) with infinitely many local conservation laws and has been used to model 1d gases of ultracold Rubidium atoms which retain memory of their initial state 23 .
A brief summary of the paper and its organization follows. We begin in §II by giving the Lagrangian (in terms of Clebsch variables) and Hamiltonian formulations and equations of motion (EOM) of adiabatic R-gas dynamics in 3d. The mass, momentum, energy and entropy equations are all expressed in conservation form. In §III, we specialize to 1d and discuss the special case of constant entropy (isentropic/barotropic) flow in which case the velocity equation also acquires a conservation form. Sound waves are discussed in §IV A and shown to be governed at long wavelengths by a cubic dispersion relation similar to that of the linearized KdV equation. In §IV B, the local conservation laws are used to reduce the determination of steady and traveling wave solutions in 1d to a single quadrature of a generalization of the Ermakov-Pinney equation. A mechanical analogy and phase plane analysis is used to show that the only such non-constant bounded solutions are cavitons (in density) and periodic waves. While these results hold for any value of γ, for γ = 2, closed-form sech 2 and cnoidal wave solutions are obtained, physically interpreted and compared with the corresponding KdV solutions. Aside from overall scales, steady solutions are parametrized by a pair of dimensionless shape parameters: a Mach number and a curvature. A parabolic embedding and a virial theorem for steady flows are given in Appendix C. In §V A, the weak form of the R-gas dynamic equations is given, and in §V B an attempt is made to find a steady shock-like profile by patching half a caviton with a constant solution. However, it is shown that there are no such continuous profiles that satisfy all the Rankine-Hugoniot conditions, though it may be possible to satisfy the mass flux condition alone. To study more general time-dependent solutions of R-gas dynamics and the evolution of initial conditions that could lead to shock-like discontinuities, we set up in §VI, a semi-implicit spectral numerical scheme for the isentropic R-gas dynamic equations with periodic boundary conditions (BCs) in 1d. For γ = 2, our numerical solutions indicate that our regularization evades the gradient catastrophe through the formation of a pair of solitary waves at the top and bottom of a velocity profile with steep negative gradient. Though we do not observe a KdVlike solitary wave train, these solitary waves can suffer collisions and approximately re-emerge with a phase shift. We also observe a rapid decay of energy with mode number and recurrent behavior with Rayleigh quotient fluctuating between bounded limits, indicating an effectively finite number of active Fourier modes. In §VII we use a canonical transformation to reformulate 3d adiabatic R-gas dynamics in terms of a complex scalar field coupled to an entropy field and three Clebsch potentials. For isentropic potential flows, this formulation shows that R-gas dynamics for any γ reduces to a defocusing 3d NLSE. In §VII A, the regularized Bernoulli equation is used to show that steady R-gas dynamic solutions map to solutions of NLSE with harmonic time dependence, with the γ = 2 caviton in 1d corresponding to the dark soliton of the cubic NLSE. In §VII B we relate the conserved quantities and bounded Rayleigh quotient of NLSE to their R-gas dynamic analogues. This connection lends credence to our numerical observations, since the cubic NLSE with periodic BCs in 1d is known to possess an infinity of conserved quantities in involution 24 . We also note in §VII C that the negative pressure γ = 2 isentropic R-gas dynamic equations in 1d are equivalent to the vortex filament and Heisenberg magnetic chain equations. We conclude with a discussion in §VIII.

II. HAMILTONIAN AND LAGRANGIAN FORMULATIONS OF 3D R-GAS DYNAMICS
It is well-known 1 that adiabatic dynamics of an ideal gas with constant specific heat ratio γ = c p /c v is governed by the continuity, momentum and internal energy equations with the temperature in energy units given by T = mp/ρ for a molecular mass m. In adiabatic flow, specific entropy (per unit mass) is advected (D t s ≡ ∂ t s + v · ∇s = 0), while the entropy per unit volume is locally conserved, ∂ t (ρs) + ∇ · (ρsv) = 0. Though the terms 'reversibly adiabatic' and 'isentropic' are often used interchangeably, in this paper we use adiabatic for D t s = 0 and isentropic for the special case where s is a global constant. For adiabatic flow, ρ and p may be taken as independent variables with s being a function of them. For a polytropic gas, s = c v log ((p/p)/(ρ/ρ) γ ) wherep,ρ are reference values. These equations follow from the Hamiltonian and Hamilton's equationsḟ = { f , H} using the (non-zero) non-canonical Poisson brackets (PB) 25 where w = ∇ × v is the vorticity. Our conservative regularization involves adding a density gradient term to the Hamiltonian while retaining the same PBs: The density gradient energy, which could arise from capillarity 13,14 , has been chosen ∝ (∇ρ) 2 to ensure positivity, parity conservation and to prevent discontinuities in density, so as to conservatively regularize shock-like discontinuities. It involves the capillarity coefficient β (ρ) = β * /ρ, where β * is a constant with dimensions L 4 T −2 . β * can be taken as λ 2 c 2 where λ is a short-distance cut-off and c a typical speed. This is the simplest form for β (ρ) that ensures the mass, momentum and energy equations are all in conservation form for the ideal mass and momentum densities. It also leads to other nice properties such as a transformation to the NLSE for isentropic potential flow. The continuity and entropy equations following from (4) and (3) are as in the ideal model. The momentum and consequently the velocity equation however, differ due to the presence of a capillary force term β * F: Remarkably, β * F = β * ∇Φ is a gradient, so that for barotropic flow (∇p/ρ = ∇h), it augments the specific enthalpy h → h + β * Φ. Thus, the vorticity evolves exactly as in ideal gas dynamics (in other words, we only regularize the 'potential' part of the velocity and don't deal with vortical singularities as in 6,7 ). Thus, Kelvin's theorem would apply in R-gas dynamics, unchanged. The momentum and velocity equations may be expressed in terms of a regularized stress tensor: The scalar part of σ defines a regularized pressure p * which includes the Gross 'quantum pressure' 19 : The energy equation for the energy density E defined in (4) is given by: The fact that (8) is in local conservation form follows from the PB formulation. Indeed, {H, H} = 0 implies that E t = {E , H} must be a divergence. The internal energy per unit volume is therefore These regularization terms in the pressure, enthalpy and internal energy depend upon density gradients and are therefore not strictly thermodynamic properties of the gas, any more than the regularized stress tensor. They are conservative analogues of the viscous stress tensor which depends on velocity gradients in dissipative gas dynamics. Interestingly, the potential Φ in (5) is also the Bohm potential U 18 that arises as a correction to the classical potential V in the quantum-corrected Hamilton-Jacobi equation for the Schrödinger wavefunction ψ = √ ρe iS/h : Our regularized stress σ also resembles the Korteweg stress σ Kor of Ref. 14,16 . Indeed, if β = β * /ρ, However, though σ Kor i j has the term (β * /ρ)(∂ i ρ)(∂ j ρ) in common with σ i j (6), they are not quite equal. Thus, though the qualitative physical features of our equations may be similar to those of the Korteweg equations, ours additionally possess some remarkable mathematical properties facilitating the analysis in this paper.
Finally, if the flow domain is all of R 3 , then β * can be scaled out by defining R = r/ β * and T = t/ β * , just as we may eliminate the dispersion coefficient in KdV on the whole real line. By contrast, in the presence of a characteristic length scale l, β * cannot be scaled out and l/λ serves as a conservative analogue of the Reynolds number.

A. Lagrangian formulation via Clebsch variables
To obtain a Lagrangian for R-gas dynamics, we use the where ε(ρ, s) is the ideal specific internal energy (9). The EOM (5) follow as the Euler-Lagrange equations (EL) for the Bateman-Thellung 28,29 Lagrangian density linear in velocities 30 augmented by the density gradient energy: The EL equations for α and λ imply the advection of s and µ, while that for φ is the continuity equation and that for s and µ are the evolution equations α t + ∇ · (αv) = ρT and λ t + ∇ · (λ v) = 0. The regularization only affects the EL equation for ρ. Upon using p = ρ 2 ∂ ε/∂ ρ, it becomes the time-dependent Bernoulli equation for adiabatic R-gas dynamics: Using these, one obtains (5) for v. There are of course related Lagrangians for the same EOM, e.g., Thus, we may interpret φ , λ and α as Lagrange multipliers enforcing the EOM for ρ, µ and s.

A. Hamiltonian and equations of motion
In what follows, we will primarily be interested in 1d adiabatic R-gas dynamics where ρ, s and p are independent of two of the Cartesian coordinates and v = (u(x,t), 0, 0). The non-zero PBs (3) simplify as w = 0: {u, u} = 0 and The total mass ( ρ dx), entropy ( ρ s dx) and more generally ρΣ(s)dx for any Σ(s) are Casimirs of this algebra. As before, the dynamics is generated by a Hamiltonian that involves a capillary energy where β (ρ) will be chosen by requiring that the momentum equation be in conservation form. The continuity and entropy equations are as in the ideal model: Thus, even with our regularization we continue to have D t p = c 2 s D t ρ where c 2 s = (∂ p/∂ ρ) s = γ p/ρ. The regularized momentum and velocity equations are The simplest way for the former to be in conservation form is for the momentum density to equal ρu and for the regularization term to be a divergence. β = β * /ρ is the simplest capillarity coefficient that ensures this, giving We note that the apparently simpler choice of constant β leads to a KdV-like ρ xxx term in the velocity equation but prevents the momentum equation from being in conservation form.
Our regularization amounts to modifying the pressure p → p * in the momentum and velocity equations: It is instructive to compare our velocity equation with Korteweg's. For capillarity coefficient β = β * /ρ, the 1d Korteweg velocity equation following from (11) is Unlike our force per unit mass β * f which is a gradient, β * f Kor is not. Thus, in the barotropic case of §III B where p x /ρ = h x , our velocity equation (21) (but not Korteweg's) comes into conservation form as in ideal gas dynamics. Finally, our energy equation is also in local conservation form, It takes a compact form in terms of a regularized specific internal energy ε * and enthalpy h * : The internal energy equation may be interpreted as the 1 st law of thermodynamics for adiabatic flow: Evidently, the gas does work against the pressure p * as well as a new type of reversible, non-dissipative work due to the regularization while ensuring that the specific entropy s is constant along the flow.
The action corresponding to (14) possesses three obvious symmetries: (a) constant shift in φ , (b) space translation and (c) time translation, leading to the local conservation laws for mass (ρ), momentum (ρu) and energy ( 1 2 ρu 2 +ρε * ) densities (19,21,24). In addition, under an infinitesimal Galilean boost (t → t, x → x − ct), the fields transform as leading to a change in the Lagrangian (16) by a spatial derivative δ L 2 = ct (p − β * ρ xx ) x . The corresponding Noether charge and flux densities are Here, we sum over χ = ρ, φ , α, s, λ and µ. The resulting conservation law ∂ t j t + ∂ x j x = 0 or (ρ(x −tu)) t + (xρu −tF p ) x = 0 involves an explicitly time-dependent Galilean charge and flux, where F p = ρu 2 + p − β * ρ (log ρ) xx is the regularized momentum flux (21). Thus, G = (x − tu) ρ dx is conserved even though {G, H} = P where P = ρu dx is the total momentum. P, G and H satisfy a 1d Galilei algebra with the total mass M furnishing a central extension: B. Isentropic R-gas dynamics of a polytropic gas Sans entropy sources/sinks and boundaries, one is mainly interested in cases where s =s is initially constant and by (19), independent of time. Thus, we consider isentropic flow where p and ρ satisfy the barotropic relation is a constant that encodes the constant value of entropy and labels isentopes. A feature of isentropic flow is that in addition to the continuity, momentum and energy equations, the velocity equation is also in conservation form: Here ε = Kρ γ−1 and h = γKρ γ−1 are specific internal energy and enthalpy. The corresponding fluxes are In ideal gas dynamics F u = F e /F m , but no such algebraic relation holds when β * = 0. The emergence of a 4 th conservation law in the isentropic case is tied to the global constancy of entropy. These equations follow from the degenerate Landau PBs {ρ, ρ} = {u, u} = 0 and {ρ(x), u(y)} = ∂ y δ (x − y) whose Casimirs include M = ρ dx and u dx. These PBs be- As above, the local conservation laws for mass, momentum, energy and Galilei charge follow from Noether's theorem. However, the conservation law u t +F u x = 0 (30) does not arise from a symmetry via Noether's theorem. This is because C = u dx is a Casimir, it acts trivially on all observables: Alembert-like splitting of the pulse. Dispersion and nonlinearity modify the shape and produce 'forerunners' and 'backrunners'. The evolution is for γ = 2 and ε = 0.1 with periodic BCs using the scheme of §VI with n max = 20 Fourier modes, a time step ∆ = 0.01 and 'nonlinearity strength' δ = 0.1.

B. Steady and traveling waves in one-dimension
Traveling waves are those where ρ, u, p and s are functions only of (x − ct), where c is the velocity of the wave. The entropy equation s t + us x = 0 becomes (u − c) s = 0. Thus, either s =s is a constant in space and time or u = c. In the former case, we have isentropic flow. In the latter, s can be an arbitrary function of (x − ct), but the fluid is at rest ('aerostatic') in a frame moving at velocity c. We will focus on the first possibility and look at steady solutions, subsequently 'boosting' them to get traveling waves.

Isentropic steady solutions
For steady flow (c = 0) the mass, momentum and velocity fluxes (31) are constant: Moreover, the steady continuity equation uρ x + u x ρ = 0 implies that the constant energy flux of Eqn. (31) is not independent: F e = F m F u . Eliminating u = F m /ρ we get two expressions for ρ xx : Taking a linear combination allows us to eliminate the ρ γ term and arrive at the second order equation In Appendix C, a different linear combination that eliminates the ρ 2 x /ρ term is considered, leading to additional results. The current choice makes it easier to treat all values of γ in a uniform manner. Interpreting x and ρ as time and position, this describes a Newtonian particle of mass β * moving in a (linear + harmonic + logarithmic) potential V on the positive half-line subject also to a 'velocity-dependent' force ∝ ρ 2 x /ρ. This ensures that the motion is 'time-reversal' (x → −x) invariant. The qualitative nature of trajectories is elucidated via a ρ-ρ x phase plane analysis in Appendix A. There are only two types of non-constant bounded solutions for ρ(x): solitary waves of depression (cavitons) and periodic waves. The latter correspond to closed trajectories around an elliptic fixed point (O-point) in the phase portrait while cavitons correspond to the homoclinic separatrix orbit that encircles an Opoint and begins and ends at a hyperbolic X-point to its right. The location of these fixed points are determined (for any γ) by the roots of the quadratic V (ρ) = 0 whose discriminant In the generic non-aerostatic situation (i.e. u ≡ 0 or equivalently F m = 0), the only cases when we get non-constant bounded solutions for ρ are (a) F p , F u > 0: both periodic solutions and cavitons and (b) F u < 0: only periodic solutions.
Remark: Eqn. (40) is a generalization of the Ermakov-Pinney equation 31,32 which corresponds to V (ρ) = ρ 4 /2, γ = 2 and β * = 1. This leads to an alternate approach to understanding (40), since the transformation z 2 = 1/ρ converts it into a Newton equation with sextic potential and no velocity dependent force for any γ: Reduction to quadrature: Subtracting the two equations in (39), we get a first order ordinary differential equation (ODE) for ρ: The 'potential energy' U is related to the potential V via ρ γ+1 U (ρ) = V (ρ). This allows us to reduce the determination of the steady density profile ρ(x) to quadrature: For integer γ ≥ 1 (43) is a hyperelliptic integral though it reduces to an elliptic integral when γ = 2 (see §IV B 3). For other values of γ, steady solutions may be found via the parabolic embedding of Appendix C.

Nondimensionalizing the steady equation
To integrate (42), it is convenient to replace the four constants (F m , F p , F u , K) with two dimensionless shape parameters (κ • , M • ) and two dimensional ones (ρ • , c • ) that set scales. These parameters are adapted to the solutions one seeks to find: ρ • is the density at a point x • where ρ x = 0. For a caviton, ρ • can be the asymptotic or trough density while for a periodic wave, it can be the trough or crest density (or the trough density of an unbounded solution with the same K). This choice will simplify the expressions for the constant fluxes (38) and where p • and c • are the pressure and sound speed at x • . We may use c 2 • to trade β * for a regularization length Here, measures the curvature of the density profile at x • = λ • ξ • . A virtue of this nondimensionalization is that unlike the four parameters in (42), only two parameters κ • and M • appear in (45). Eqn. (45) describes zero energy trajectoriesρ(ξ ) of a unit mass particle in the potential −T (ρ). Evidently, the allowed values ofρ must lie between adjacent positive zeros of T with T > 0 in between (Fig. 2). To obtain (45), we used the following expressions for the constant fluxes and entropy: Conversely, we may invert (46) by first determining ρ • by solving the algebraic equation following from (42). The remaining new parameters follow from (46) and (38): For γ = 2, T (ρ) (45) becomes a cubic with rootsρ = 1 and The density of periodic and solitary waves must lie between adjacent positive roots with T > 0 in between. Interestingly, it can be shown that if for γ = 2, all three roots of T are positive, then the same holds for any 1 < γ < 2. So some qualitative features of solutions for γ = 2 are valid more generally. For γ = 2, the nature of solutions on the κ • -M • plane (Fig. 3a) changes when the two relevant roots coalesce, i.e., when the discriminant of the cubic T vanishes: Interestingly, when we reinstate dimensions, the periodic solutions from regions (a), (b) and (c) of the κ • -M • plane (Fig 3a) are physically identical. They differ by the choice of nondimensionalizing density ρ • which could be any one of the roots of the cubic in (47). Thus, the parameters (κ • , M • , c • , ρ • ) generically furnish a 3-fold cover (redundancy) of the original space of constants ((F m ) 2 , F u , F p , K). For solitary waves, it degenerates into a double cover: the two families of cavitons ((i) and (ii)) in Fig. 3a differ via the choice of ρ • as trough or asymptotic density. Moreover, M • is the Mach number at the trough in (i) while it is the asymptotic Mach number in (ii). In a caviton, the flow goes from asymptotically subsonic to supersonic at x = 0. A caviton is like a pair of normal shocks joined at the trough. Finally, the map between the two sets of parameters becomes a 1-fold cover for the aerostatic cavitons at κ • = M • = 0, since their trough densities vanish and ρ • can only be chosen as the asymptotic density.
In light of the above remarks, we now restrict attention to the yellow triangular region (c) where κ • > 0 and 0 Here, the roots of T are 0 <ρ − <ρ + < 1 and (45) is reduced to quadrature: Here, F(φ |m) = φ 0 (1 − m sin 2 θ ) −1/2 dθ is the incomplete elliptic integral of the first kind with m the square of the elliptic modulus k (see 17.4.62 of 33 ). Inverting, we writeρ(ξ ) in terms of the Jacobi cn(u, m) function: (52) Here ξ − = ξ (ρ − ). The periodic wave extends from a trough atρ − to a crest atρ + with amplitude and wavelength Here K(m) = F( π 2 |m) is the complete elliptic integral of the 1 st kind. When we approach the left boundary κ • → 0 + with 0 < M • < 1, the wavelength diverges (K(1/2κ • ) ∼ √ κ • log κ • for small κ • ) and the periodic waves turn into cavitons which extend from a trough densityρ =ρ − = M 2 • < 1 to an asymptotic densityρ =ρ + = 1: (54) On the other hand, when we approach the upper boundary M • = 1 − √ 2κ • , A → 0 and we get constant solutions. By contrast, on the lower boundary (M • = 0, 0 < κ • < 1 2 ) we continue to have periodic solutions except thatρ vanishes at the trough (ρ − = 0) while the crest densityρ + = 1 − 2κ • : The dimensional ρ, u and p, are got by reinstating the con- (57) for isentropic flow. Reversing the sign of M • reverses the flow direction leaving p and ρ unaltered. Moreover, sinceρ ≥ 0 and u = F m /ρ, the flow is unidirectional with positive velocity solitary waves being waves of elevation in u and vice-versa. A caviton is superficially a bifurcation of the constant solution ρ(x) ≡ ρ(±∞). However, though the caviton and constant solutions have the same constant specific entropy, they have different values of mass and energy (per unit length). For instance, the energy density (18) of an aerostatic caviton is less than that of the constant state: (58) As a consequence, the uniform state cannot, by any isentropic time-dependent motion with fixed BCs at ±∞, in the absence of sources and sinks, evolve via R-gas dynamics to the caviton, or vice versa, since the two states have different invariants. However, one could imagine creating, say an aerostatic caviton, by starting with a uniform state and introducing a point sink at the origin, to suck fluid out. A symmetrical pair of expansion waves would travel to infinity on both sides, without affecting the conditions at infinity. When the density reaches 0 at the origin, we stop the sink. The pressure gradient will then be balanced by the regularizing density gradient force, and the solution should tend to the aerostatic caviton as t → ∞. Since temperature T = Km(γ − 1)ρ γ−1 , the caviton corresponds to a region of width λ where the temperature drops. Loosely, the regularizing force is like Pauli's exchange repulsion, capable of maintaining a depression in density with variable but isentropic temperature and pressure distributions.

Traveling waves of isentropic R-gas dynamics
Here we generalize the above steady solutions to waves traveling at speed c: (ρ, u, p)(x,t) = (ρ, u, p)(x − ct). The continuity, velocity and momentum equations (30) are readily integrated, giving the constant fluxes: Eliminating u = c + F m /ρ and taking a linear combination leads us to a 2 nd order nonlinear ODE for ρ: Proceeding as in §IV B 1 and §IV B 2, this may be reduced to the nonlinear first order equation Comparing with (45), we see that the passage from steady to traveling waves involves only a shift in the square of the Mach number M 2 Here,c is the speed of the traveling wave in units of the sound speed c • at the point where ρ (x − ct) = 0. Thus, for each value of β * , we have a 5-parameter (κ • , M • , ρ • , c • ,c) space of traveling cavitons and periodic waves. The dimensionful ρ, u and p for any value of (M • −c) are given by Traveling cavitons for γ = 2: In these nondimensional variables, traveling cavitons have the profilẽ Here,ρ ± are got by replacing M 2 • → (M • −c) 2 in (49). Their wavelength is given by the same formula (53). Our cnoidal waves are very similar to those of KdV (u t − 6uu x + u xxx = 0): (65) Here ξ = x − ct and f (ξ 3 ) = f 3 and f 1,2,3 are the roots of f 3 + 1 2 c f 2 + A f + B with A and B constants of integration.

A. Weak form of R-gas dynamic equations
The R-gas dynamic equations ( (19) and (21)) involve u x , p x , ρ x , ρ xx and ρ xxx . Thus, classical solutions need to be C 1 in u and p and C 3 in ρ. However, by multiplying the conservation equations by C ∞ test functions φ , ψ and ζ and integrating by parts, we arrive at a weak form of the equations: (ρu) t ψ − p + ρu 2 + β * ρ 2 Thus, for weak solutions, it suffices that ρ be C 1 and u and p be merely continuous.

B. Steady shock-like profile from half a caviton
Here, we try to use the steady solutions of §IV B to model the structure of a normal shock propagating to the right in the lab frame. As in Fig. 4, in the shock frame, the shock is assumed to be located around x = 0. The undisturbed preshock medium is to the right (x > 0) while the 'disturbed' post-shock medium is to the left (x < 0) 1 . As x → ±∞ the variables approach the asymptotic values ρ ±∞ , u ±∞ and p ±∞ with ρ +∞ < ρ −∞ . The Rankine-Hugoniot (RH) conditions are obtained by equating the conserved fluxes F m (19), F p (20) and F e (24) at x = ±∞: In our γ = 2 cavitons (54), the flow is subsonic at x = ±∞ and supersonic at x = 0. We exploit this in trying to find a shocklike steady solution by patching half a caviton with a constant solution. Thus, we seek a steady solution where ρ(x) ≡ ρ ∞ in the pre-shock medium and is the left half of a caviton for x < 0. The half-caviton profile is got from (54) by taking the reference location x • = −∞ so that ρ • = ρ −∞ and κ • = 0: Here On the other hand, ρ +∞ must correspond to a constant solution with fluxes F m,p,e ∞ . It can take one of two density values corresponding to the X/O point (A3): We now attempt to patch these two solutions requiring ρ and ρ x to be continuous at x = 0. Since ρ has a local minimum at the trough of a caviton, ρ x (0 − ) = 0 = ρ x (0 + ) of the undisturbed medium. However, a difficulty arises in trying to ensure that ρ is continuous at x = 0. Indeed, suppose we impose the RH conditions, and F e −∞ = F e ∞ , then both the pre-and post-shock regions correspond to a common phase portrait. We observe (see also Fig. 9a) that the caviton trough density ρ(0 − ) = ρ −∞ M 2 −∞ is strictly less than both ρ X,O ∞ . Thus, the post-shock semi-caviton solution cannot be continuously extended into the pre-shock region.
Stated differently, in the above patched shock construction, if ρ +∞ is chosen to be equal to ρ(0 − ) in order to make ρ continuous, then the RH conditions are violated. Let us illustrate this with a γ = 2 aerostatic example. We take the pre-shock region to be vacuum (ρ = u = p ≡ 0 for x > 0) and try to patch this at x = 0 with the left half of the aerostatic caviton [ρ(x) = ρ −∞ tanh 2 (x/2λ −∞ ) and u ≡ 0] of (56). This caviton corresponds to the values F m −∞ = F e −∞ = κ −∞ = M −∞ = 0 and has trough density ρ(0 − ) = 0 with trough density gradient ρ x (0 − ) = 0 as well. Since the caviton is isentropic, its trough pressure p(0 − ) = Kρ(0 − ) 2 = 0. Thus, ρ and p are both C 1 at x = 0, while u ≡ 0 so that this is a weak solution in the sense of §V A. However, it violates the RH conditions: while F m ≡ 0 is continuous, F p and F u are not. In fact, in the pre-shock vacuum region F p = F u = 0 while in the post-shock region they are non-zero, as evaluating them at x = −∞ shows: (70) In the post-shock region K − = F p −∞ /ρ 2 −∞ = 0 whereas in the pre-shock vacuum K is arbitrary since p = ρ = 0 for x > 0. In conclusion, there are no continuous steady shock-like solutions in the shock frame that satisfy the RH conditions. To see how initial conditions (ICs) that would lead to shocks in the ideal model are regularized, we turn to a numerical approach.

A. Spectral method with nonlinear terms isolated
In this section, we discuss the numerical solution of the isentropic R-gas dynamic initial value problem (IVP). It is convenient to work with the nondimensional variablesρ,û, s, etc. of §IV A. The continuity equation isρˆt + (ρû)x = 0 while for isentropic flow,ŝ is a global constant which may be taken to vanish by adding a constant to entropy. Thus, we can eliminatep =ρ γ in the velocity and momentum equations, both of which are in conservation form (30): The energy equation is ∂ˆtÊ +F ê x = 0 where the energy density and flux arê These equations follow from the PB {ρ(x),û(ŷ)} = ∂ŷδ (x−ŷ) and the Hamiltonian We will consider ICs that are fluctuations around a uniform state. For stability of the numerics, we separate the linear and nonlinear terms in the equations and treat the former implicitly and the latter explicitly. Introducing the book-keeping parameter δ (which will also enter through the ICs and may eventually be set to 1), we writê (74) We will consider flow in the domain −π ≤x ≤ π with periodic BCs and thus expandρ andũ as (75) Since ρ dx is conserved, ρ 0 can be taken time-independent. Furthermore, we choose the constantρ used to nondimensionalize ρ as the (conserved) average density, so that ρ 0 = 0. We also suppose that the 'background' flow velocity u is independent of positionx. Since û dx is conserved and û dx = 2π(u(t) + δ u 0 (t)), we may absorb δ u 0 (t) into u(t) and thereby take u 0 = 0. Next, we write the continuity equation with nonlinear terms isolated: The velocity equation (71) in conservation form is Separating out the linear part we get In (76) and (78) the linear terms are at O(δ 0 ) while the nonlinearities involve δ to higher powers, depending on the value of γ. Interestingly, for γ = 2, the pressure gradient doesn't contribute to the nonlinear part of the acceleration: Expanding F m = ∑ n F m n e inx and F u = ∑ n F u n e inx , the EOM in Fourier space becomė When nonlinearities (F m , F u ) are ignored and we assume (ρ n (t), u n (t) ∝ e −iω n t ), one finds the dispersion relation (ω n − nu) 2 = n 2 (1 + ε 2 n 2 ) familiar from §IV A. To deal with the fully nonlinear evolution given by (80), we use a semi-implicit numerical scheme outlined in Appendix D.
We find that at early times, whereû has a negative slope, its gradient increases and decreases where the slope is positive. Without the regularization (ε = 0), the higher Fourier modes can then get activated and the velocity and density profiles become highly oscillatory with steep gradients. Moreover, amplitudes begin to grow and the code eventually ceases to conserve energy and momentum. By contrast, in the presence of the regularization (say ε = 0.2), we find that the above gradient catastrophe is averted and energy is conserved (to about 3 parts in 1000) while momentum is conserved to machine precision. In fact, we find that the real and imaginary parts Q r and Q i of the next conserved quantity (97) are also conserved. Interestingly, when u develops a steep negative gradient, a pair of solitary waves emerge at the top (wave of elevation) and bottom (wave of depression) of the u profile and the gradient catastrophe is avoided (see Fig. 5). This mechanism by which the incipient shocklike discontinuity is regularized is to be contrasted with KdV, where an entire train of solitary waves can form 34 . However, as with KdV, our solitary waves can suffer a head-on collision and pass through each other. After the collision, they re-emerge with roughly similar shapes and a phase shift. Fig.  6 shows the space-time trajectories of the centers of several of these solitary waves, showing their collisions. We also find that higher modes u n and ρ n (n 9) grow from zero but soon saturate and remain a few orders of magnitude below the first few modes (see Fig. 7). This justifies truncating the Fourier series at n max = 16. The modes also display an approximate periodicity in time. This suggests recurrent motion [35][36][37] . This behavior is also captured in Fig. 8a where we plot the Rayleigh quotient or mean square mode number as a function of time. It is found to fluctuate between bounded limits indicating that effectively only a finite number of modes participate in the dynamics. Another interesting statistic is the spectral distribution of energy (E n ) and its dependence on time. Fig. 8b, shows the time evolution of log E n for n = 2, 6, 10, 16 for the ICρ = cos 2x andũ = 0 and demonstrates that the energy in the higher modes remains small. Moreover, each mode ρ n and u n oscillates between an upper and lower bound and shows approximate periodicity with differing periods. In Fig. 8c, log E n vs n is plotted for a few values oft. Unlike the power law decay n −5/3 in the inertial range of fully developed turbulence, here we see that E n drops exponentially with n. In particular, there is no equipartition of energy among the modes.
Specializing to isentropic potential flow where s =s is constant, p = K(γ − 1)ρ γ (29) and v = ∇φ , (83) simplifies to Using the above PBs for ψ, one finds that ψ satisfies the defocusing nonlinear Schrödinger equation Interestingly, in 1d, the R-gas dynamic form of the focusing cubic (γ = 2) NLSE had been obtained in the context of the Heisenberg magnetic chain 39,40 . However, as noted in 40 , the Heisenberg chain leads to negative pressure! Returning to (85), we see that the real and imaginary parts of the NLSE correspond to the Bernoulli and continuity equations. The ∇ 2 ψ term leads to the divergence of the mass flux, v 2 and regularization terms in the isentropic R-gas dynamic equations (86) Evidently, β * plays the role ofh 2 . The nonlinear term (γK/2)|ψ| 2(γ−1) ψ corresponds to the isentropic pressure p = (γ − 1)Kρ γ whose positivity implies we get the defocusing/repulsive NLSE. Thus, our regularization term 2β * (∇ 2 √ ρ)/ √ ρ is like a quantum correction to the classical isentropic pressure. For γ = 2 we get the cubic NLSE or Gross-Pitaevskii equation (without an external trapping potential). Note that 1d isentropic flow on the line with v = (u(x), 0, 0) is always potential flow: u = φ x . So the above transformation takes 1d isentropic R-gas dynamics (30) to the defocusing 1d NLSE, which for γ = 2 and periodic BCs admits infinitely many conserved quantities in involution 24 . This explains our numerical observations of approximate phase shift scattering of solitary waves and recurrence.
A. NLSE interpretation of steady R-gas dynamic cavitons and cnoidal waves It turns out that steady solutions of 1d isentropic R-gas dynamics ( §IV B) correspond to NLSE wavefunctions ψ with harmonic time dependence. For γ = 2, our aerostatic caviton corresponds to the dark soliton of NLSE. More generally, aerostatic steady solutions correspond to ψ of the form ρ(x) exp(−iF u t/2 β * ) where F u is the constant velocity flux (31). Finally, non-aerostatic cnoidal waves correspond to interesting asymptotically plane wave NLSE wavefunctions modulated by a periodic cnoidal amplitude. Here, we consider the cavitons and cn waves in increasing order of complexity. Aerostatic caviton: The simplest caviton solution of §IV B 3 is aerostatic: Here, ρ • , λ • and c 2 • are positive constants. The corresponding specific enthalpy and velocity flux are Thus, the Bernoulli equation φ t + F u = 0 (31) is satisfied provided we take the velocity potential φ = −c 2 • t to be timedependent. The resulting ψ is the dark soliton solution of the defocusing NLSE (see §6.6 of 2 ) Non-aerostatic caviton: More generally, for a non-aerostatic caviton (for 0 < M • < 1 and ρ • , λ • and c 2 • positive constants) In this case, the constant velocity flux is F u = c 2 • (2 + M 2 • )/2. The resulting velocity potential is Thus, ψ is asymptotically a plane wave with phase speed c • (2 + M 2 • )/2M • . This ψ may be regarded as a highfrequency carrier wave modulated by a localized signal. Aerostatic snoidal waves: The simplest steady periodic solutions is the aerostatic snoidal wave (55): Here, F u = c 2 • (1 − κ • ) and φ = −F u t. The resulting ψ is a snoidal wave with harmonic time dependence: Non-aerostatic cnoidal waves: Finally, in the general case, we have from (52), where 0 < κ • < 1/2 and 0 < M • < 1 − √ 2κ • (lower triangular region of Fig. 3a). Hereρ ± (κ • , M • ) are as in (49). Furthermore, ρ • , c 2 • and λ • are positive constants that set scales. As before, φ = −F u t + x 0 u(x ) dx , where the constant velocity flux F u depends on c • , κ • and M • , but is independent of ρ • and λ • . Though an explicit formula for φ (x,t) is not easily obtainable, it is evident that for large x, φ grows linearly in x with a subleading oscillatory contribution. Thus, ψ has a purely harmonic time-dependence exp −iF u t/2 β * , a periodic cnoidal modulus |ψ| = ρ(x) and an argument that asymptotically grows linearly in x. Thus, asymptotically, ψ is a plane wave modulated by the periodic amplitude ρ(x).

B. Conserved quantities and Rayleigh quotient of NLSE and R-gas dynamics
The cubic 1d NLSE admits an infinite tower of conserved quantities. The first three are N = |ψ| 2 dx, P NLSE = ℑ(ψ * ψ x ) dx, and These correspond to the mass M = N, momentum P = 2 β * P NLSE and energy H = 2β 3/4 * E NLSE of R-gas dynamics. The next conserved quantity of NLSE (with periodic BCs) is 24 ℜQ and ℑQ correspond to the following in R-gas dynamics: In §VI A we use the conservation of Q r and Q i to test our numerical scheme. In Ref. 35 , it was shown that for the cubic 1d NLSE, the Rayleigh quotient or mean square mode number for periodic BCs, is related to the number of active degrees of freedom and is bounded in both focussing and defocussing cases. This has a simple interpretation in R-gas dynamics, since Consequently, 1d isentropic R-gas dynamic potential flows are recurrent in the sense discussed in 35 .

C. Connection to vortex filament and Heisenberg chain equations
Intriguingly, the negative pressure γ = 2 isentropic R-gas dynamic equations in 1d for ρ and u (30) are equivalent to the vortex filament equationξ = G ξ × ξ where G is a constant. Here the curve ξ (x,t) represents a vortex filament with tangent vector ξ . The equivalence is most transparent when the vortex filament equation is expressed in the Frenet-Serret frame, with curvature κ = √ ρ and torsion τ = u 39 . The FS frame (ξ , n, b) is an orthonormal basis along the filament, where n = ξ /κ is the unit normal and b = ξ × n the binormal. The FS equations describe how the frame changes along the filament. ξ = κn, n = −κξ + τb and b = −τn (100) We may use the vortex-filament and FS equations to find evolution equations for the FS frame: Eqn. (101) leads to evolution equations for κ (usingṅ · n = 0 as n · n = 1) and τ = n · b: (102) Taking G = 1/2, κ = √ ρ and τ = u, (102) reduces to the continuity and velocity equations for γ = 2 isentropic R-gas dynamics with β * = 1/2, but with a negative p = −ρ 2 /8 40 . Furthermore, it is well-known 41 that the vortex filament equation is related to the Heisenberg magnetic chain equatioṅ S = G S × S with S = ξ . An open question is to give a geometric or magnetic chain interpretation of the positive pressure R-gas dynamic equations as well as those for general γ. It is noteworthy that negative pressures (relative to atmospheric pressure) can arise in real flows, for instance in the presence of strong currents 42 .

VIII. DISCUSSION
It is a significant feature of our attempt to conservatively regularize singularities in gas dynamics that it has led us (in the case of isentropic potential flow) to the defocusing NLSE. Heuristically, the defocusing interaction tends to amplify linear dispersive effects and thereby prevent blowups. For γ = 2, this connection to the cubic NLSE should provide powerful tools [including the inverse scattering transform in 1d (see §9.10 of 2 , 24 and 43 )] to deal with the initial-boundary value problems in various dimensions as well as alternative numerical schemes. Moreover, the bound on the NLSE Rayleigh quotient obtained in 35 (see also §VII B) generalizes to 2d and 3d as well as to nonlinearities other than cubic (γ = 2). This would have implications for recurrence in more general Rgas dynamic isentropic potential flows, even in the absence of integrability. The techniques of dispersive shock wave theory 1,4,44 could provide additional tools to address R-gas dynamic flows. In 45,48 , a classification of semilinear PDEs [perturbations of linear equations by lower order nonlinear terms] into subcritical, critical and supercritical, based on conserved quantities (mass, energy etc.), scaling symmetries and regularity of initial data is described. Though the equations of R-gas dynamics given in §II are not semilinear, in the special case of isentropic potential flow, the transformation to NLSE makes them semilinear. It is thus interesting to examine the implications of this classification for R-gas dynamics. For example, the critical scaling regularity 45 of NLSE in d spatial dimensions is s c = d/2 − 1/(γ − 1). Thus, if the initial data is such that the number of particles and energy are finite (so that ψ, ∇ψ ∈ L 2 and ψ ∈ H 1 ), then according to the scaling heuristic, the NLSE is subcritical for any γ > 1 in 1d and 2d and also for 1 < γ < 3 in 3d.
In §V B we argued that 1d R-gas dynamics does not admit any smooth or continuous shock-like steady solutions. In fact, we found that if we try to patch half a caviton at its trough density with a constant state, then the mass, momentum and energy fluxes in the pre-shock region cannot all match their values in the post-shock region, so that the Rankine-Hugoniot conditions are violated. We conjecture that this absence of steady shock-like solutions is a general feature of a wide class of conservatively regularized gas dynamic models. Loosely, this is like d'Alembert's 'theorem' that continuous solutions of Euler's equations cannot ever lead to drag, although possibly to lift. On the other hand, inclusion of viscosity does permit drag as well as steady shock-like solutions as in the Burgers equation 1,2 . Allowing for non-steady solutions, we find that in R-gas dynamics, the gradient catastrophe is averted through the formation of a pair of solitary waves (see §VI B). It would be interesting to see if this mechanism is observed in any physical system, say, one where dissipative effects are small as in nonlinear optics, weak shocks, cold atomic gases or superfluids. For further discussion on steepening gradients and criteria for detecting wave breaking in dispersive hydrodynamics, see 46,47 .
Though we have formulated R-gas dynamics in 3d, our analytic and numerical solutions have been restricted to 1d. It would be interesting to extend this work to higher dimensional problems such as oblique shocks and the Sedov-Taylor spherical blast wave problem. The mechanical and thermodynamic stability of our traveling caviton and periodic wave so-lutions is also of interest. One also wishes to examine whether the capillarity energy considered here arises from kinetic theory in a suitable scaling limit of small Knudsen number as for the Korteweg equation 16,49 . Finally, our Hamiltonian and Lagrangian formulations of R-gas dynamics can be used as starting points in formulating the quantum theory. The transformation to NLSE provides another approach to quantization for isentropic potential flow especially when γ = 2. Though we have focused on the conservatively regularized model, a more complete and realistic treatment would have to include viscous dissipation just as in the KdV-Burgers equation.
Our attempt to generalize KdV to include the adiabatic dynamics of density, velocity and pressure has led to a interesting link between KdV and NLSE that is quite different from the well known ones (see E.g. 2 ). In fact, we may view R-gas dynamics as a natural generalization of both. While it extends the KdV idea of a minimal conservative dispersive regularization to adiabatic gas dynamics in any dimension and shares with it the cnoidal and sech 2 solutions, it also reduces to the defocusing cubic NLSE for isentropic potential flow of a gas with adiabatic index two. Thus, the cubic 1d NLSE is a very special member of a larger class of R-gas dynamic equations that make sense in any dimension and for nonlinearities other than cubic while also allowing for adiabatically evolving entropy and vorticity distributions. Isentropic aerostatic steady solutions: In the aerostatic limit (u ≡ 0), both the fluxes F e and F m vanish, though their ratio F e /F m = F u is finite. Eqn. (40) for steady solutions becomes β * ρ xx = −V (ρ) + (γ + 1)β * 2 In this limit, the small-ρ logarithmic barrier in V (A1) is absent, and the singularity along ρ = 0 becomes 'naked'. One of the FPs in Eqn. (A3) tends to (0, 0), while the other one tends to (γF p /(γ − 1)F u , 0). Table I(b) summarizes the nature of physical fixed points and bounded solutions for various possible signs of F p and F u . Interestingly, for F p = 0 and F u < 0, there is a family of solitary waves of elevation, though with negative pressure.
The Euler-Lagrange equation reduces to (40). Thus, we have Hamiltonian and Lagrangian formulations for both the full Rgas dynamic field equations and their reduction to the space of steady solutions.
Appendix C: Parabolic embedding and Lagrange-Jacobi identity for steady flow For γ = 2 the quadrature in (43) cannot be done using elliptic functions, but could be done numerically. An alternative approach is to take a linear combination of the equations in (39) to obtain a form of the steady equation for ρ without the velocity dependent term but with a generally non-integral power of ρ: If we introduce a pseudo-time τ, then steady solutions can be obtained via a parabolic embedding in a nonlinear heat equation with a source: By prescribing suitable BCs and starting with an arbitrary initial condition, the solution of this PDE should relax to the stable solutions of (C1). Eqn.
(C1) may also be used to derive additional virial/Lagrange-Jacobi-type identities by multiplying it by ρ and using (42) and repeating the process. The first two such identities are These integral identities can provide valuable checks on any numerics used to obtain steady solutions.