text:
# 1 Historical remark and Introduction

## 1 Historical remark and Introduction

In 1943-44, von Neumann became convinced that the calculation of the flows of compressible fluids containing strong shocks could be accomplished only by numerical methods. He conceived the idea of capturing shocks, i.e., of ignoring the presence of a discontinuity. Employing a Lagrangian description of compressible flow, setting heat conduction and viscosity equal to zero, von Neumann replaced space and time derivatives by symmetric difference quotients. Calculations using this scheme were carried out; the approximation resulting from these calculations (see ) showed oscillations on the mesh scale behind the shock. Von Neumann boldly conjectured that the oscillations in velocity represent the heat energy created by the irreversible action of the shock, and that as $`\mathrm{\Delta }x`$ and $`\mathrm{\Delta }t`$ tend to zero, the approximate solutions tend in the weak sense to the discontinuous solution of the equations of compressible flow.

In  it was counterconjectured that von Neumann was wrong in his surmise, i.e., that although the approximate solutions constructed by his method do converge weakly, the weak limit fails to satisfy the law of conservation of energy.

In  J.Goodman and P.Lax investigated von Neumann’s algorithm applied to the scalar equation

$$u_t+uu_x=0$$
(1)
(it is called the Hopf equation ), in the semidiscrete case. Using numerical experimentation and analytical techniques the demonstrated the weak convergence of the oscillatory approximations, and that the weak limit fails to satisfy the scalar equation in question.

Von Neumann’s dream of capturing shocks was realized in his joint work with Richtmyer in 1950, see . Oscillations were eliminated by the judicious use of artificial viscosity; solutions constructed by this method converge uniformly except in a neighborhood of shocks, where they remain bounded and are spread out over a few mesh intervals. The limits appear to satisfy the conservation laws of compressible flow. The conservation of mass and momentum is the consequence of having approximated these equations by difference equations in conservation form; but the von Neumann-Richtmyer difference approximation to the energy equation is not in conservation form.

In the paper  T.Hou and P.Lax compared the results of a von Neumann-Richtmyer calculation with the weak limit of of calculations performed by von Neumann’s original method.

The aim of this paper, at first, it is to introduce the idea of understanding of conservation laws; the second is to proposed a numerical method for calculation of hydrodynamic shocks profile *without using differences schemes*.

We believe that there is no discontinues in the nature of shocks. Discontinues solution of hydrodynamic equations only a rough mathematical model of shocks. When viscosity is taken into account, for example, the shocks are seen to be smeared out, so that the pure mathematical surface of discontinuety are replaced by the thin layers ($`10^7`$$`10^8`$ m) in which pressure, density, temperature and etc. vary rapidly but continuosly (see Fig.1).

The equations of compressible flow in one space dimension can be written in the following Lagrangian form:

$$\begin{array}{c}u_t+p_x=0\hfill \\ v_tu_x=0\hfill \\ e_t+pu_x=0\hfill \end{array}$$
(2)
It is a classical approach. Here $`u`$ is velocity, $`p`$ is pressure, $`v`$ specific volume, and $`e`$ internal energy, connected with $`p`$ and $`v`$ via an equation of state. Here we mostly study the case when the equation of state is given by the $`\gamma `$-law with $`\gamma =1.4`$

$$e=\frac{pv}{\gamma 1}.$$

Substituting that into the (2) we get the following equation

$$\frac{1}{\gamma 1}\left(pv\right)_t+pu_x=0.$$
(3)
The first equation of (2) is conservation of momentum, the second equation conservation of mass, the third the work equation.

We will consider a steady-state shock. Imagine a long pipe containing a fluid initially in equilibrium (thermally and mechanically), into which a piston is pushing from one end with constant speed, as shown in Fig.1. In the presence of dissipation the specific volume, $`v`$, and the fluid velocity, $`u`$ are as shown by the curves. *Our problem is to calculate the exact shape of this curves.* The shock is steady, at least approximately, after it has gone to a sufficiently great distance from the initiating piston. Then $`u,v,e`$, etc. depend on $`x`$ and $`t`$ only through the combination $`y=xct,`$ where $`c`$ is a speed of the shock relative to the original, or Lagrangian, coordinates.

Now, we are going to study the equations (2) in specific sense. Namely, we are going to rewrite mentioned equations in the sense of $`𝐑\epsilon `$–distributions. We give the definition of the special kind of solutions of the some conservation laws in the sense of $`𝐑\epsilon `$-distributions and consider the method for the numerical calculations of the smooth shocks and soliton like solutions of the Hopf equation and equations of compressible flow in the mentioned sense. This method based on orthogonal system of the Hermite functions as a base for calculation of such solutions (i.e. shocks and infinitely narrow solitons). Calculations of profiles of infinitely narrow soliton and shock wave are reduced to the nonlinear system of algebraic equations in $`𝐑^{n+1}`$, $`n>1`$. We proved, using the Schauder fixed point theorem , that the mentioned system has at least one solution in $`𝐑^{n+1}`$. We showed that there is possibility to find out some of the solutions of this system using the Newton iteration method . We considered examples and numerical tests. We also should emphasis that proposed numerical approach do not use a difference scheme (see also ).

Let us consider a bit of theory which we will apply to conservation laws.

## 2 Non-Archimedean field of Laurent series and $`𝐑\epsilon `$–distributions.

The theory of Non-Archimedean fields was considered in the book by A.H.Lightstone and A.Robinson .

###### Definition 2.1

A *Laurent series* is a formal object

$$\underset{n=0}{\overset{\mathrm{}}{}}\xi _{n+k}\epsilon ^{n+k}$$

where $`k`$ is a fixed (i.e., fixed for this Laurent series), each $`\xi _i𝐑`$, and either $`\xi _k0`$ or each $`\xi _i=0`$.

The Laurent series $`_{n=0}^{\mathrm{}}\xi _n\epsilon ^n`$, where $`\xi _0=1`$ and $`\xi _n=0`$ if $`n>0`$, is denoted by 1. It is easy to see that the Laurent series is a field. Let us denote it by $`𝐑\epsilon `$. The norm on the field of Laurent series can define

$$|x|_\nu =e^{\nu (x)}\text{for}\text{each}x𝐑\epsilon $$

(in place of $`e`$ can use any number greater than 1). The function $`\nu (x)`$ is a Non-Archimedean valuation. Define

$$\nu (0)=\mathrm{}\text{and}\nu \left(\underset{n=0}{\overset{\mathrm{}}{}}\xi _{n+k}\epsilon ^{n+k}\right)=k\text{if}\underset{n=0}{\overset{\mathrm{}}{}}\xi _{n+k}\epsilon ^{n+k}0,\xi _k0.$$

The norm $`||_\nu `$ have properties

1. $`|x|_\nu =0\text{if and only if}x=0,`$
2. $`|xy|_\nu =|x|_\nu |y|_\nu ,`$
3. $`|x+y|_\nu \mathrm{max}\{|x|_\nu ,|y|_\nu \}.`$

Here, we propose a general construction of the $`𝐑\epsilon `$–valued generalized functions , . These objects are a natural generalization of Sobolev-Schwartz distributions. We call them as $`𝐑\epsilon `$–distributions.

1. Consider all functions $`f(x,\epsilon )C^{\mathrm{}}(𝐑\times (0,1))`$ such that integrals

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}f(x,\epsilon )\psi (x)𝑑x$$

exist for any $`\epsilon `$ and for all $`\psi (x)`$ from a given class of functions $`𝒳`$ ($`𝒳`$ can be $`C_0^{\mathrm{}}(𝐑),`$ $`𝒮(𝐑)`$ and etc.).
2. Suppose also that $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}f(x,\epsilon )\psi (x)𝑑x`$ is a number $`a_{f,\epsilon }(\psi )`$ from the field of Laurent series $`𝐑\epsilon `$.
3. The two functions $`f(x,\epsilon )`$ and $`g(x,\epsilon )`$ call equivalent with respect to test functions $`𝒳`$ if and only if

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}f(x,\epsilon )\psi (x)𝑑x=a_{f,\epsilon }(\psi )=a_{g,\epsilon }(\psi )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}g(x,\epsilon )\psi (x)𝑑x.$$

The equality means in sense of the field of Laurent series $`𝐑\epsilon `$ for all functions $`\psi 𝒳`$. Classes of equivalent functions call $`𝐑\epsilon `$*–functions*. The expression

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}f(x,\epsilon )\psi (x)𝑑x$$

associates a number from $`𝐑\epsilon `$ with every $`\psi `$. Such a quantity is called a functional. In this case a linear functional map $`𝒳`$ into the Non-Archimedean field $`𝐑\epsilon `$. Call these functionals as $`𝐑\epsilon `$*-distributions*.

Thus,

###### Proposition 2.1

$`𝐑\epsilon `$–function $`f(x,\epsilon )=0`$ if and only if

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}f(x,\epsilon )\psi (x)𝑑x=0𝐑\epsilon $$

for every $`\psi `$ from $`𝒳`$.

The set of all $`𝐑\epsilon `$*–distributions* denote by $`(𝒳)`$

###### Remark 2.1

Recall that the idea of representation of a function $`fL_{\mathrm{𝑙𝑜𝑐}}^1(𝐑)`$ in terms of a linear functional

$$C_0^{\mathrm{}}(𝐑)\psi \underset{\mathrm{}}{\overset{+\mathrm{}}{}}f(x)\psi (x)𝑑x𝐑$$

based on well-known proposition that if $`fL_{\mathrm{𝑙𝑜𝑐}}^1(𝐑)`$ and $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}f(x)\psi (x)𝑑x=0`$ for any $`\psi C_0^{\mathrm{}}(𝐑)`$ then $`f=0`$ almost everywhere.

Let us consider an example of the $`𝐑\epsilon `$-distribution.

###### Example 2.2

Take $`𝒳=C_0^{\mathrm{}}(𝐑)`$ and $`f(x,\epsilon )=\phi (x/\epsilon )`$, $`\phi (x)C_0^{\mathrm{}}(𝐑)`$ then $`𝐑\epsilon `$–distribution can write in the following form.

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi (x/\epsilon )\psi (x)𝑑x=\epsilon \underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi (x)𝑑x\psi (0)+\epsilon ^2\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x\phi (x)𝑑x\frac{\psi ^{}(0)}{1!}+\mathrm{}.$$

Note that $`\phi (x/\epsilon )`$ converges to the function

$$u(x)=\{\begin{array}{cc}\phi (0),\hfill & \text{if}x=0,\hfill \\ 0,\hfill & \text{if}x0.\hfill \end{array}$$

Last function almost everywhere equals to zero.

Like Sobolev-Schwartz distributions we can *differentiate* $`𝐑\epsilon `$–distributions. For example,

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{d}{dx}\phi (x/\epsilon )\psi (x)𝑑x=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi (x/\epsilon )\frac{d}{dx}\psi (x)𝑑x,$$

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi (x/\epsilon )\frac{d}{dx}\psi (x)𝑑x=\epsilon \underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi (x)𝑑x\psi ^{}(0)\epsilon ^2\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x\phi (x)𝑑x\frac{\psi ^{\prime \prime }(0)}{1!}\mathrm{}.$$

It is evident that $`𝐑\epsilon `$–distributions are more general objects than Sobolev-Schwartz distributions , .

## 3 Conservation laws. Non-Archimedean approach.

A conservation law asserts that the rate of change of the total amount of substance contained in a fixed domain $`G`$ is equal to the flux of that substance across the boundary of $`G`$. Denoting the density of that substance by $`u`$, and the flux by $`f`$, the conservation law is

$$\frac{d}{dt}_Gu(t,x)𝑑x=_Gf\stackrel{}{n}𝑑S.$$

Applying the divergence theorem and taking $`d/dt`$ under the integral sign we obtain

$$_G(u_t+\mathrm{𝐝𝐢𝐯}f)𝑑x=0.$$

Dividing by vol $`(G)`$ and shrinking $`G`$ to a point where all partial derivatives of $`u`$ and $`f`$ are continuous we obtain the differential conservation law

$$u_t(t,x)+\mathrm{𝐝𝐢𝐯}f(u(t,x))=0.$$

Note, that if $`f(u)=u^2/2`$ then we obtained the Hopf equation (1). In general, previous calculations based on the following well known proposition.

###### Proposition 3.1

If $`GL_{\mathrm{𝑙𝑜𝑐}}^1(𝐑)`$ and $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}G(x)\psi (x)𝑑x=0`$ for any $`\psi C_0^{\mathrm{}}(𝐑)`$ then $`G=0`$ almost everywhere.

###### Definition 3.1

Let us consider two sets of the smooth functions, depending on a small parameter $`\epsilon (0,1]`$. Let us take all functions $`v(t,x,\epsilon )`$ which have the type

$$v(t,x,\epsilon )=l_0+\mathrm{\Delta }l\phi \left(\frac{xct}{\epsilon }\right),$$

$`l_0,\mathrm{\Delta }l,c`$ are real numbers, $`\mathrm{\Delta }l0`$ and $`\phi 𝒮(𝐑)`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\phi (y)𝑑y=1`$. We denote this set of functions by $`I`$. We call $`I`$ as a set of infinetely narrow solitons.

###### Definition 3.2

Now, let us take all functions $`w(t,x,\epsilon )`$ which have the type

$$w(t,x,\epsilon )=h_0+\mathrm{\Delta }hH\left(\frac{xat}{\epsilon }\right),$$

$`h_0,\mathrm{\Delta }h,a`$ are real numbers, $`\mathrm{\Delta }h0`$ and $`H(x)={\displaystyle \underset{\mathrm{}}{\overset{x}{}}}\theta (y)𝑑y`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\theta (y)𝑑y=1`$ and $`\theta 𝒮(𝐑)`$. We denote this set of functions by $`J`$. We call $`J`$ as a set of shock waves.

It is natural to consider conservation laws as an integral expressions which contain the time $`t`$ as parameter. Therefore, we introduce the following concept.

###### Definition 3.3

The function $`vI`$ (or $`wJ`$) will be a solution of the Hopf equation up to $`e^l`$, $`l𝐍_0`$ in the sense of $`𝐑\epsilon `$–distributions if for any $`t[0,T]`$

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{v_t(t,x,\epsilon )+v(t,x,\epsilon )v_x(t,x,\epsilon )\right\}\psi (x)𝑑x=\underset{k=l}{\overset{+\mathrm{}}{}}\xi _k\epsilon ^k𝐑\epsilon ,$$
(4)
$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{w_t(t,x,\epsilon )+w(t,x,\epsilon )w_x(t,x,\epsilon )\right\}\psi (x)𝑑x=\underset{k=l}{\overset{+\mathrm{}}{}}\eta _k\epsilon ^k𝐑\epsilon $$
(5)
for every $`\psi 𝒮(𝐑).`$ In case when $`l`$ is equal to $`+\mathrm{}`$ the function $`v(t,x,\epsilon )`$ (or $`w(t,x,\epsilon )`$) exactly satisfies the Hopf equation in the sense of $`𝐑\epsilon `$–distributions.

Certainly, one can consider instead of the Hopf equation some conservation law.

From mathematical point of view, we deal with a infinitely differentiable functions in definitions 3.1 and 3.2, so that we avoid the problem of distribution multiplication. From physical point of view, functions from the set $`I`$ or $`J`$ can describe fast processes. Mathematical models of such processes based on functions from $`I`$ or $`J`$ may give additional information and take in account a short zone where physical system make a jump from one position to another.

Thus, we will consider solutions of the Hopf equation which are infinitely narrow solitons or shock waves. It easy to see that

$$v(t,x,\epsilon )\{\begin{array}{cc}l_0+\mathrm{\Delta }l\phi (0),\hfill & ifx=ct,\hfill \\ l_0,\hfill & ifxct.\hfill \end{array}\text{as}\epsilon 0$$

$$w(t,x,\epsilon )h_0+\mathrm{\Delta }hH(xat),\text{as}\epsilon 0$$

$`H`$ is Heaviside function.

## 4 Method for the numerical calculations of the microscopic profiles of soliton like solutions of the Hopf equation in the sense of $`𝐑\epsilon `$–distributions.

Thus, conservation laws are integral expressions. Therefore, it is natural, that one can interpret the Hopf equation in the sense of the definition 3.3.

We will seek a solution of the Hopf equation in the type of infinitely narrow soliton, i.e. let us $`vI`$. Substitute $`v(t,x,\epsilon )`$ into integral expression (4) using the following formulas

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{}{t}\left\{\phi \left(\frac{xct}{\epsilon }\right)\right\}\psi (x)𝑑x=\underset{k=0}{\overset{+\mathrm{}}{}}c\epsilon ^km_k\frac{1}{k!}\psi ^{(k+1)}(ct),$$
(6)
$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{}{x}\left\{\frac{1}{2}\phi ^2\left(\frac{xct}{\epsilon }\right)\right\}\psi (x)𝑑x=\underset{k=0}{\overset{+\mathrm{}}{}}\epsilon ^kg_k\frac{1}{k!}\psi ^{(k+1)}(ct).$$
(7)
We denote

$$m_k(\phi )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}y^k\phi (y)𝑑y,g_k(\phi )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}y^k\frac{\phi ^2(y)}{2}𝑑y,k=0,1,2,\mathrm{}.$$
(8)
Thus, we obtain

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{v_t+vv_x\right\}\psi 𝑑x=\underset{k=0}{\overset{+\mathrm{}}{}}\left\{\mathrm{\Delta }l(cl_0)m_k(\mathrm{\Delta }l)^2g_k\right\}\epsilon ^k\frac{\psi ^{(k+1)}(ct)}{k!}.$$
(9)
From the last expression we have conditions for the function $`\phi (x)`$. Namely,

$$g_k(\phi )\frac{cl_0}{\mathrm{\Delta }l}m_k(\phi )=0,k=0,1,2\mathrm{}.$$
(10)
From the first ($`k=0`$) we have

$$\frac{cl_0}{\mathrm{\Delta }l}=\frac{g_0}{m_0}=\frac{1}{2}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi ^2(x)𝑑x.$$
(11)
Hence, we can rewrite conditions (10) as follows.

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi ^2(x)𝑑x\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\phi (x)𝑑x=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\phi ^2(x)𝑑x,k=0,1,2\mathrm{}.$$
(12)
Now, let us prove the following lemma.

###### Lemma 4.1

For any non-negative integer $`n`$ exists such function $`\phi 𝒮(𝐑)`$, $`\phi 0`$ which satisfies the following system of non-linear equations:

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\phi (x)𝑑x=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\phi ^2(x)𝑑x/\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi ^2(x)𝑑xk=0,1,2\mathrm{}n.$$
(13)
Proof. First, we will seek function $`\phi (x)`$ in the following type:

$$\phi (x)=c_0h_0(x)+c_1h_1(x)+\mathrm{}+c_nh_n(x),$$
(14)
where

$$h_k(x)=\frac{H_k(x)}{\sqrt{2^kk!}\sqrt[4]{\pi }}e^{x^2/2}\mathrm{are}\mathrm{Hermit}\mathrm{functions}.$$
(15)
Then we substitute the expression (14) into conditions (13). After that we will have nonlinear system of $`n+1`$ equations with $`n+1`$ unknowns ($`c_0,c_1,c_2,\mathrm{},c_n`$). We write this system by the following way.

$$A\stackrel{}{x}=𝒩(\stackrel{}{x}),\stackrel{}{x}=(c_0,c_1,\mathrm{},c_n)$$
(16)
$`A`$ is a matrix with elements

$$A_{kj}=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^kh_j(x)𝑑x=(𝐢)^j𝐢^k\sqrt{2\pi }h_j^{(k)}(0),𝐢=\sqrt{1},k,j=0,1,2,\mathrm{}n$$

$`𝒩`$ is nonlinear map such that

$$𝒩(\stackrel{}{x})=\frac{1}{\stackrel{}{x}^2}\underset{k=0}{\overset{n}{}}(N(k)\stackrel{}{x},\stackrel{}{x})\stackrel{}{e}_k\underset{k=0}{\overset{n}{}}f_k(\stackrel{}{x})\stackrel{}{e}_k$$
(17)
Vector $`\stackrel{}{e}_k=(e_0,e_1,\mathrm{},e_n)`$ such that $`e_k=1`$ and $`e_j=0`$ for all $`jk`$. $`N(k)`$ are matrices with elements

$$N_{ij}(k)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^kh_i(x)h_j(x)𝑑x,i,j,k=0,1,2\mathrm{}n$$
(18)
and functions

$$f_k(\stackrel{}{x})=\frac{(N(k)\stackrel{}{x},\stackrel{}{x})}{\stackrel{}{x}^2}.$$

Note that functions $`f_k(\stackrel{}{x})`$ are continuous everywhere except $`\stackrel{}{x}=0`$ and $`|f_k(\stackrel{}{x})|N(k)`$ due to Cauchy-Bunyakovskii unequality. Matrix $`A`$ is invertible for any $`n`$ because of $`det(A)`$ is a Wronskian for the linear independent system of Hermit functions $`h_0(x)`$, $`h_1(x)`$, … $`h_n(x)`$ and

$$det(A)=(2\pi )^{\frac{(n+1)}{2}}W(h_0(0),h_1(0),\mathrm{}h_n(0)).$$

We can write the system (16) as

$$\stackrel{}{x}=\underset{k=0}{\overset{n}{}}f_k(\stackrel{}{x})A^1\stackrel{}{e}_kF(\stackrel{}{x})\mathrm{or}\stackrel{}{x}=A^1(𝒩(\stackrel{}{x}))F(\stackrel{}{x}).$$
(19)
Let us describe the function $`F:𝐑^{n+1}𝐑^{n+1}`$. It is continuous except $`\stackrel{}{x}=0`$ and bounded. Indeed,

$$F(\stackrel{}{x})A^1\underset{k=0}{\overset{n}{}}N(k),r_n=A^1\underset{k=0}{\overset{n}{}}N(k).$$
(20)
Let us consider function $`𝒩(\stackrel{}{x})`$. It is continuous function everywhere in $`𝐑^{n+1}`$ except $`\stackrel{}{x}=0`$ and, moreover, $`𝒩(𝐑^{n+1}\backslash \{0\})\mathrm{\Pi }_1`$ where $`\mathrm{\Pi }_1=\{\stackrel{}{z}𝐑^{n+1}:z_0=1\}`$ is a plane. Further $`A^1(\mathrm{\Pi }_1)=\mathrm{\Pi }_2`$ where $`\mathrm{\Pi }_2=\{\stackrel{}{y}𝐑^{n+1}:\underset{k=0}{\overset{n}{}}a_{0j}y_j=1\}`$ is another plane.

$$A_{0j}=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}h_j(x)𝑑x=(𝐢)^j\sqrt{2\pi }h_j(0),𝐢=\sqrt{1},j=0,1,2,\mathrm{}n$$

Thus, we can consider the function $`F(\stackrel{}{x})`$ which is defined on the convex compact set $`C_n=\mathrm{\Pi }_2B[0,r_n]`$ such that $`F:C_nC_n`$, where $`B[0,r_n]`$ is a closed ball with radius $`r_n`$. Function $`F`$ is continuous on the $`C_n`$ because of $`\stackrel{}{0}C_n`$. Now we can use J.Schauder theorem.

###### Theorem 4.2 (Schauder fixed-point theorem )

Let $`C`$ be a compact convex subset of a normed space $`E`$. Then each continuous map $`F:CC`$ has at least one fixed point.

Hence, we can conclude that our system (19) and therefore system (16) has at least one solution. Thus, there is a function $`\phi (x)`$ which satisfy to conditions (13) proposed lemma.

###### Remark 4.1

Let us a function $`\phi (x)`$ satisfies lemma condition. If $`\beta 𝐑`$ then the function $`\phi (x+\beta )`$ also satisfies lemma condition. Moreover, if $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\phi ^2(x)𝑑x=\alpha `$ then $`\phi (\alpha x)`$ satisfies lemma condition.

Thus, we can formulate the following result.

###### Theorem 4.3

For any integer $`l`$ there is a infinitely narrow soliton type solution of the Hopf equation (in the sense of the definition 3.3) up to $`e^l`$ with respect to the norm $`||_\nu `$, i.e.

$$v(t,x,\epsilon )=l_0+\mathrm{\Delta }l\phi \left(\frac{xct}{\epsilon }\right),$$
(21)
$`l_0,\mathrm{\Delta }l,c`$ are real numbers, $`\mathrm{\Delta }l0`$ and $`\phi 𝒮(𝐑)`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\phi (y)𝑑y=1`$. Moreover,

$$\frac{cl_0}{\mathrm{\Delta }l}=\frac{1}{2}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\phi ^2(x)𝑑x.$$
(22)
For example, calculations in case $`l=7`$ give the “profile” $`\phi (x)`$ (see Fig. 3) for the infinitely narrow soliton $`v(t,x,\epsilon )=\phi \left(\frac{xct}{\epsilon }\right)`$:

$$\phi (x)=\left\{\frac{c_0}{\sqrt[4]{\pi }}+\frac{c_2(4x^22)}{\sqrt{2^22!}\sqrt[4]{\pi }}+\frac{c_4(16x^448x^2+12)}{\sqrt{2^44!}\sqrt[4]{\pi }}\right\}e^{x^2/2},$$
(23)
where $`c_0=0.66583`$, $`c_2=0.23404`$, $`c_4=0.05028`$, $`c=0.25032`$ ($`c`$ is a velocity of the soliton). Numbers $`c_0`$, $`c_2`$, $`c_4`$ and $`c`$ were found approximately by iteration method using the following sequence.

$$\stackrel{}{x}_{m+1}=A^1(𝒩(\stackrel{}{x}_m)),m=0,1,2,\mathrm{}.$$
(24)
Matrix $`A`$ and a nonlinear $`𝒩`$ were introduced in the lemma proof.

Calculations of soliton-like profiles $`\phi (x)`$ for the Hopf equation in case $`l=13,`$ $`15,`$ $`17,`$ $`19,`$ $`21`$ give us pictures (Fig. 3, 5, 5, 7, 7).

For the $`l`$ greater than $`21`$ matrix $`A`$ is close to singular and calculations can be inaccurate.

## 5 Calculations of the microscopic profiles of the shock wave solutions of the Hopf equation in the sense of $`𝐑\epsilon `$–distributions.

A solution of the Hopf equation in this case we will seek in the set $`J`$. Namely,

$$w(t,x,\epsilon )=h_0+\mathrm{\Delta }hK\left(\frac{xat}{\epsilon }\right),$$

$`h_0,\mathrm{\Delta }h,a`$ are real numbers, $`\mathrm{\Delta }h0`$ and

$$K(x)=\underset{\mathrm{}}{\overset{x}{}}\theta (y)𝑑y,\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\theta (y)𝑑y=1,\theta 𝒮(𝐑).$$

Substitute $`w(t,x,\epsilon )`$ into the integral expression (5) using the following formulas

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\frac{}{t}\left\{K\left(\frac{xat}{\epsilon }\right)\right\}\psi (x)𝑑x=\underset{k=0}{\overset{+\mathrm{}}{}}(a)\epsilon ^km_k\frac{\psi ^{(k)}(at)}{k!},$$
(25)
$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}K\left(\frac{xat}{\epsilon }\right)\frac{}{x}\left\{K\left(\frac{xat}{\epsilon }\right)\right\}\psi (x)𝑑x=\underset{k=0}{\overset{+\mathrm{}}{}}\epsilon ^kr_k\frac{\psi ^{(k)}(at)}{k!}.$$
(26)
We denote by

$$m_k(\theta )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}y^k\theta (y)𝑑y,r_k(\theta )=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)\left(\underset{\mathrm{}}{\overset{x}{}}\theta (y)𝑑y\right)𝑑x,k=0,1,2,\mathrm{}.$$
(27)
Thus, we get

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{w_t+ww_x\right\}\psi 𝑑x=\underset{k=0}{\overset{+\mathrm{}}{}}\left\{(\mathrm{\Delta }h)^2r_k\mathrm{\Delta }h(ah_0)m_k\right\}\epsilon ^k\frac{\psi ^{(k)}(at)}{k!}.$$
(28)
From the last expression we have conditions for the function $`\theta (x)`$

$$r_k(\theta )\frac{ah_0}{\mathrm{\Delta }h}m_k(\theta )=0,k=0,1,2\mathrm{}.$$
(29)
From the first ($`k=0`$) we have

$$\frac{ah_0}{\mathrm{\Delta }h}=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\theta (x)\left(\underset{\mathrm{}}{\overset{x}{}}\theta (y)𝑑y\right)𝑑x=\frac{1}{2}.$$
(30)
Therefore, we can rewrite (29) as

$$\frac{1}{2}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)𝑑x=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)\left(\underset{\mathrm{}}{\overset{x}{}}\theta (y)𝑑y\right)𝑑xk=0,1,2\mathrm{}.$$
(31)
The same method one can prove that there is such function $`\theta (x)𝒮(𝐑)`$ which satisfies the following conditions

$$\frac{1}{2}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)𝑑x=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)\left(\underset{\mathrm{}}{\overset{x}{}}\theta (y)𝑑y\right)𝑑xk=0,1,2\mathrm{}n.$$
(32)
Thus, we can formulate next result.

###### Theorem 5.1

For any integer $`l`$ there is a shock wave type solution of the Hopf equation (in the sense of the definition 3.3) up to $`e^l`$ with respect to the norm $`||_\nu `$.

$$w(t,x,\epsilon )=h_0+\mathrm{\Delta }hK\left(\frac{xat}{\epsilon }\right),$$
(33)
$`h_0,\mathrm{\Delta }h,a`$ are real numbers, $`\mathrm{\Delta }h0`$ and $`K(x)={\displaystyle \underset{\mathrm{}}{\overset{x}{}}}\theta (y)𝑑y`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\theta (y)𝑑y=1`$ and $`\theta 𝒮(𝐑).`$ Moreover,

$$\frac{ah_0}{\mathrm{\Delta }h}=\frac{1}{2}.$$
(34)
Note that the condition (34) is *Rankine — Hugoniot condition* for the velocity of a shock wave.

As in previous section we seek function $`\theta (x)`$ in the following type:

$$\phi (x)=a_0h_0(x)+a_1h_1(x)+\mathrm{}+a_nh_n(x),$$
(35)
where $`h_k(x)`$ are Hermite functions. Calculations in case $`l=7`$ give the following “profile” ($`K(x)`$) for the shock wave $`w(t,x,\epsilon )=K\left({\displaystyle \frac{xat}{\epsilon }}\right)`$ (where $`h_0=0,\mathrm{\Delta }h=1`$).

$$K(x)=\underset{\mathrm{}}{\overset{x}{}}\left\{\frac{c_0}{\sqrt[4]{\pi }}+\frac{c_2(4\tau ^22)}{\sqrt{2^22!}\sqrt[4]{\pi }}+\frac{c_4(16\tau ^448\tau ^2+12)}{\sqrt{2^44!}\sqrt[4]{\pi }}\right\}e^{\tau ^2/2}𝑑\tau $$
(36)
where $`c_0=0.79617`$, $`c_2=0.53004`$, $`c_4=0.17923`$, $`c=1/2`$ is a velocity of the shock wave (see Fig. 9). Numbers $`c_0`$, $`c_2`$, $`c_4`$ were found approximately.

Note that the function $`K(x)`$ is not unique. There is a different function $`K_1(x)`$ which satisfies mentioned above conditions. It has the following type

$$\begin{array}{c}K_1(x)=\underset{\mathrm{}}{\overset{x}{}}\left\{\frac{c_0}{\sqrt[4]{\pi }}+\frac{c_12\tau }{\sqrt{2^11!}\sqrt[4]{\pi }}+\frac{c_2(4\tau ^22)}{\sqrt{2^22!}\sqrt[4]{\pi }}\right\}e^{\tau ^2/2}𝑑\tau +\hfill \\ +\underset{\mathrm{}}{\overset{x}{}}\left\{\frac{c_3(8\tau ^312\tau )}{\sqrt{2^33!}\sqrt[4]{\pi }}+\frac{c_4(16\tau ^448\tau ^2+12)}{\sqrt{2^44!}\sqrt[4]{\pi }}\right\}e^{\tau ^2/2}𝑑\tau \hfill \end{array}$$
(37)
where $`c_0=0.18357`$, $`c_1=0.73567`$, $`c_2=0.74733`$, $`c_3=0.15327`$ $`c_4=0.29539`$, $`c=1/2`$ is a velocity of the shock wave (see Fig. 9). Coefficients $`c_0`$, $`c_1`$, $`c_2`$, $`c_3`$, $`c_4`$ were found approximately by the Newton iteration method.

Taking in account the *Rankine — Hugoniot condition* (34) we also have graphs (Fig. 11, 11) as a shock profiles.

Here we describe how it is possible to find coefficients $`c_0,c_1,\mathrm{},c_n`$ in this case by the Newton iteration method for the following system of nonlinear equations.

$$P(\stackrel{}{c})=A\stackrel{}{c}2\underset{k=0}{\overset{n}{}}(S(k)\stackrel{}{c},\stackrel{}{c})\stackrel{}{e}_k=0,\stackrel{}{c}=(c_0,c_1,\mathrm{},c_n)$$
(38)
Vector $`\stackrel{}{e}_k=(e_0,e_1,\mathrm{},e_n)`$ such that $`e_k=1`$ and $`e_j=0`$ for all $`jk`$. $`S(k)`$ are matrices with elements

$$S_{ij}(k)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^kh_i(x)\underset{\mathrm{}}{\overset{x}{}}h_j(y)𝑑y𝑑x,i,j,k=0,1,2\mathrm{}n$$
(39)
Matrix $`A`$ have elements

$$A_{ij}=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^ih_j(x)𝑑x,i,j=0,1,2\mathrm{}n$$
(40)
We can write the formula for the Newton iteration method .

$$\stackrel{}{x}_{m+1}=\stackrel{}{x}_m\left[P^{}(\stackrel{}{x}_m)\right]^1\left[P(\stackrel{}{x}_m)\right],$$
(41)
where $`\left[P^{}(\stackrel{}{x})\right]`$ is a linear map depending on the vector $`\stackrel{}{x}`$.

$$\left[P^{}(\stackrel{}{x})\right][\stackrel{}{h}]=A\stackrel{}{h}2\left\{\underset{k=0}{\overset{n}{}}(S(k)\stackrel{}{x},\stackrel{}{h})\stackrel{}{e}_k+\underset{k=0}{\overset{n}{}}(S^T(k)\stackrel{}{x},\stackrel{}{h})\stackrel{}{e}_k\right\}$$
(42)
Calculations of shock profiles $`K(x)`$ for the Hopf equation in case $`l=8,`$ $`9,`$ $`10,`$ $`11,`$ $`12,`$ $`13`$ give us the following pictures (Fig. 13, 13, 15, 15, 17, 17). Here, we show only two different types of the shock type solutions of the Hopf equation. We can find more solutions if we take a different initial data for the Newton iteration method.

###### Remark 5.1

It is not easy to see that there is exist function

$$\theta (x)=\underset{n=1}{\overset{\mathrm{}}{}}a_nh_n(x),\stackrel{}{a}=(a_0,a_1,\mathrm{},a_n,\mathrm{})l_2,$$
(43)
such that

$$\frac{1}{2}\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)𝑑x=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\theta (x)\left(\underset{\mathrm{}}{\overset{x}{}}\theta (y)𝑑y\right)𝑑x,k=0,1,2\mathrm{}.$$
(44)
We think that it is true.

## 6 Calculations of the microscopic profiles of the shock wave solutions of the equations of compressible flow in the sense of $`𝐑\epsilon `$–distributions.

Now, we are going to study the equations of compressible flow

$$u_t+p_x=0$$
(45)
$$v_tu_x=0$$
(46)
$$\frac{1}{\gamma 1}\left(pv\right)_t+pu_x=0$$
(47)
in specific sense. Namely, we rewrite mentioned equations in the sense of $`𝐑\epsilon `$–distributions.

We will seek for a solution of the equations in the following form

$$u(t,x,\epsilon )=u_0+\mathrm{\Delta }uU\left(\frac{xat}{\epsilon }\right),$$
(48)
$`u_0,\mathrm{\Delta }u,v`$ are real numbers, $`\mathrm{\Delta }u0`$ and $`U(x)={\displaystyle \underset{\mathrm{}}{\overset{x}{}}}\stackrel{~}{U}(y)𝑑y`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\stackrel{~}{U}(y)𝑑y=1`$ and $`\stackrel{~}{U}𝒮(𝐑).`$

$$p(t,x,\epsilon )=p_0+\mathrm{\Delta }pP\left(\frac{xat}{\epsilon }\right),$$
(49)
$`p_0,\mathrm{\Delta }p,a`$ are real numbers, $`\mathrm{\Delta }p0`$ and $`P(x)={\displaystyle \underset{\mathrm{}}{\overset{x}{}}}\stackrel{~}{P}(y)𝑑y`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\stackrel{~}{P}(y)𝑑y=1`$ and $`\stackrel{~}{P}𝒮(𝐑).`$

$$v(t,x,\epsilon )=v_0+\mathrm{\Delta }vV\left(\frac{xat}{\epsilon }\right),$$
(50)
$`v_0,\mathrm{\Delta }v,a`$ are real numbers, $`\mathrm{\Delta }v0`$ and $`V(x)={\displaystyle \underset{\mathrm{}}{\overset{x}{}}}\stackrel{~}{V}(y)𝑑y`$, $`{\displaystyle \underset{\mathrm{}}{\overset{+\mathrm{}}{}}}\stackrel{~}{V}(y)𝑑y=1`$ and $`\stackrel{~}{V}𝒮(𝐑).`$ Note that $`a`$ is a velocity of the shock waves.

In the other hand, we suppose

$$\stackrel{~}{U}(x)=a_0h_0(x)+a_1h_1(x)+\mathrm{}+a_nh_n(x),\stackrel{}{a}=(a_0,a_1,\mathrm{},a_n),$$
(51)
$$\stackrel{~}{P}(x)=b_0h_0(x)+b_1h_1(x)+\mathrm{}+b_nh_n(x),\stackrel{}{b}=(b_0,b_1,\mathrm{},b_n),$$
(52)
$$\stackrel{~}{V}(x)=c_0h_0(x)+c_1h_1(x)+\mathrm{}+c_nh_n(x),\stackrel{}{c}=(c_0,c_1,\mathrm{},c_n),$$
(53)
where $`h_k(x)`$ are Hermite functions.

We understand the solution of the system in sense of $`𝐑\epsilon `$–distributions.

###### Definition 6.1

Functions $`uJ`$, $`pJ`$ and $`vJ`$ is a solution of the system (45), (46), (47) up to $`e^l`$, $`l𝐍_0`$ in the sense of $`𝐑\epsilon `$–distributions if for any $`t[0,T]`$

$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{u_t+p_x\right\}\psi (x)𝑑x=\underset{k=l}{\overset{+\mathrm{}}{}}\xi _k\epsilon ^k𝐑\epsilon ,$$
(54)
$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{v_tu_x\right\}\psi (x)𝑑x=\underset{k=l}{\overset{+\mathrm{}}{}}\zeta _k\epsilon ^k𝐑\epsilon $$
(55)
$$\underset{\mathrm{}}{\overset{+\mathrm{}}{}}\left\{\frac{1}{\gamma 1}\left\{p_tv+pv_t\right\}+p_tu_x\right\}\psi (x)𝑑x=\underset{k=l}{\overset{+\mathrm{}}{}}\eta _k\epsilon ^k𝐑\epsilon $$
(56)
for every $`\psi 𝒮(𝐑)`$.

In case when $`l`$ is equal to $`+\mathrm{}`$ functions $`u(t,x,\epsilon )`$, $`p(t,x,\epsilon )`$ and $`v(t,x,\epsilon )`$ exactly satisfies the system (45), (46), (47) in the sense of $`𝐑\epsilon `$–distributions.

Substituting $`u`$, $`p`$ and $`v`$ into (54), (55) we get the following relations for the moments.

$$a\mathrm{\Delta }um_k(\stackrel{~}{U})+\mathrm{\Delta }pm_k(\stackrel{~}{P})=0,k=0,1,2,\mathrm{}n,$$
(57)
$$a\mathrm{\Delta }vm_k(\stackrel{~}{V})\mathrm{\Delta }um_k(\stackrel{~}{U})=0,k=0,1,2,\mathrm{}n.$$
(58)
We denote as usual by

$$m_k(\stackrel{~}{U})=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\stackrel{~}{U}(x)𝑑x,m_k(\stackrel{~}{U}U)=\underset{\mathrm{}}{\overset{+\mathrm{}}{}}x^k\stackrel{~}{U}(x)\left(\underset{\mathrm{}}{\overset{x}{}}\stackrel{~}{U}(y)𝑑y\right)𝑑x,$$
(59)
where $`k=0,1,2,\mathrm{}.`$

It is easy to find $`a`$ from (57) and (57) when k=0. Indeed,

$$a=\frac{\mathrm{\Delta }p}{\mathrm{\Delta }u},a=\frac{\mathrm{\Delta }u}{\mathrm{\Delta }v}$$
(60)
These are *Rankine — Hugoniot conditions.* Indeed, from the last expressions we have

$$a^2=\frac{\mathrm{\Delta }p}{\mathrm{\Delta }v}.$$
(61)
See J. von Neumann and R.D. Richtmyer  formula (21). We also conclude that $`\mathrm{\Delta }u<0`$, $`\mathrm{\Delta }p<0`$ and $`\mathrm{\Delta }v>0`$.

Because of $`\mathrm{\Delta }u`$, $`\mathrm{\Delta }p`$ and $`\mathrm{\Delta }v`$ some real numbers, therefore, all three vectors with coordinates $`m_k(\stackrel{~}{U})`$, $`m_k(\stackrel{~}{P})`$ and $`m_k(\stackrel{~}{V})`$, $`k=0,1,2,\mathrm{}n`$, respectively should be collinear. However,

$$m_0(\stackrel{~}{U})=m_0(\stackrel{~}{P})=m_0(\stackrel{~}{P})=1.$$

Hence, $`\stackrel{}{a}=\stackrel{}{b}=\stackrel{}{c}`$.

Substituting $`u`$, $`p`$ and $`v`$ into (56) and taking in account the last equalities we get the following relations for the moments $`(k=0,1,\mathrm{}n):`$

$$\left\{p_0\mathrm{\Delta }u\frac{a(\mathrm{\Delta }pv_0+\mathrm{\Delta }vp_0)}{\gamma 1}\right\}m_k(\stackrel{~}{U})+\left\{\mathrm{\Delta }p\mathrm{\Delta }u\frac{2a\mathrm{\Delta }p\mathrm{\Delta }v}{\gamma 1}\right\}m_k(U\stackrel{~}{U})=0.$$
(62)
When $`k=0`$ we will have *Rankine — Hugoniot conditions* for our system

$$\left\{p_0\mathrm{\Delta }u\frac{a(\mathrm{\Delta }pv_0+\mathrm{\Delta }vp_0)}{\gamma 1}\right\}+\frac{1}{2}\left\{\mathrm{\Delta }p\mathrm{\Delta }u\frac{2a\mathrm{\Delta }p\mathrm{\Delta }v}{\gamma 1}\right\}=0.$$
(63)
Using the equalities (60), (63), we will have

$$\gamma p_0+\frac{\mathrm{\Delta }p}{\mathrm{\Delta }v}v_0+\frac{1}{2}(\gamma +1)\mathrm{\Delta }p=0\text{or}\frac{p_1}{p_0}=\frac{(\gamma +1)v_0(\gamma 1)v_1}{(\gamma +1)v_1(\gamma 1)v_0}.$$
(64)
Then taking in account the condition (61), we get

$$\left\{v_0+\frac{1}{2}(\gamma +1)\mathrm{\Delta }v\right\}a^2=\gamma p_0.$$
(65)
Hence,

$$a=\left(\frac{2}{(\gamma +1)\frac{\mathrm{\Delta }v}{v_0}+2}\right)^{1/2}\left(\gamma \frac{p_0}{v_0}\right)^{1/2}.$$
(66)
Thus, we get the same formula for the shock velocity as in the paper of J. von Neumann and R.D. Richtmyer  (formula (72)).

Thus, from (62) follows that

$$m_k(\stackrel{~}{U})=2m_k(\stackrel{~}{U}U),k=0,1,2,\mathrm{}n.$$

This system we already know how to solve by the Newton iteration method. See conditions (32) and solutions in this case. Thus, we can formulate the following result.

###### Theorem 6.1

For any integer $`l`$ there is a solution of the system of equations (45), (46), (47) in the sense of the definition 6.1 if

$$a=\frac{\mathrm{\Delta }p}{\mathrm{\Delta }u},a=\frac{\mathrm{\Delta }u}{\mathrm{\Delta }v},$$

$$\left\{p_0\mathrm{\Delta }u\frac{a}{\gamma 1}(\mathrm{\Delta }pv_0+\mathrm{\Delta }vp_0)\right\}+\frac{1}{2}\left\{\mathrm{\Delta }p\mathrm{\Delta }u\frac{2a}{\gamma 1}\mathrm{\Delta }p\mathrm{\Delta }v\right\}=0.$$

If we assume that $`v_0`$ and $`\mathrm{\Delta }v`$ are known then we can calculate all constants $`\mathrm{\Delta }u`$, $`\mathrm{\Delta }p`$, $`a`$ and profiles of the shocks. Let us take $`v_0=0.2`$, $`\mathrm{\Delta }v=0.8`$. We denote by $`v_1=v_0+\mathrm{\Delta }v`$, $`v_1`$ is a volume before the shock formation. The quantity $`\eta =v_1/v_0`$ is a measure of the shock strength. Denote by $`p_1=p_0+\mathrm{\Delta }p`$, $`v_1`$ is a normal pressure (it is known, $`p_1=1atm`$) before the shock formation. From the (64) we can get

$$\mathrm{\Delta }p=\frac{2\gamma p_1}{\gamma 2v_0/\mathrm{\Delta }v1}.$$

Therefore, we can calculate step by step

$$a=\sqrt{\frac{\mathrm{\Delta }p}{\mathrm{\Delta }v}},\mathrm{\Delta }u=a\mathrm{\Delta }v.$$

Shock profiles of the considered system (45), (46), (47) one can find on pictures (Fig. 19, 19). We considered only one type profile of the shock, the case when $`l=13`$, $`v_0=0.2`$, $`\mathrm{\Delta }v=0.8`$, $`p_1=1`$ and $`\gamma =1.4`$.

Finally, $`\mathrm{\Delta }p=28p_1`$, $`a=5.9`$, $`\mathrm{\Delta }u=4.7`$, $`\stackrel{}{a}=(0.47894,`$ $`0.70727,`$ $`0.31783,`$ $`0.4706,`$ $`0.53605,`$ $`0.20792,`$ $`0.4099,`$ $`0.057108,`$ $`0.16832,`$ $`0.0088166,`$ $`0.029395)`$.

In conclusion we should emphasis that our calculation method looks like the Fourier method for linear differential equations but applied to the nonlinear equations. Our method allowed to obtain all known formulas for the shocks characteristics and, in addition, find a microscopic behaviour of shocks in the thin layer. According to the our model the pressure in the thin layer can be negative. It is possible that the concept of “pressure” in the thin layer (where the jump of $`p(t,x,\epsilon )`$ took place) one should understand in special sense. We think that the role of $`\epsilon `$ in the calculations can play so-called the average length of free movement of gas molecules. From our point of view the phenomena of the formation and development of characteristics of the shock wave should describe in terms of non-Archemedean distance or non-Archemedean geometry.

In addition, we can use Laguerre functions or harmonic functions in our calculations instead of Hermite functions.

We hope that our approach will be useful for the problems of nonlinear Optics and problems of Quantum Mechanics.

Research is partially supported by Belarussian Fundamental Research Foundation Grant No F99M-082.


format: text

text:
# 𝜃 dependence of SU(𝑁) gauge theories

## 1 Introduction

It is rather well established that four-dimensional SU($`N`$) gauge theories must have a nontrivial dependence on the angle $`\theta `$ that appears in the Euclidean Lagrangian as

$$_\theta =\frac{1}{4}F_{\mu \nu }^a(x)F_{\mu \nu }^a(x)i\theta q(x)$$
(1)
where $`q(x)`$ is the topological charge density

$$q(x)=\frac{g^2}{64\pi ^2}ϵ_{\mu \nu \rho \sigma }F_{\mu \nu }^a(x)F_{\rho \sigma }^a(x).$$
(2)
Indeed, the most plausible explanation of how the solution of the so-called U(1)A problem can be compatible with the $`1/N`$ expansion (performed keeping $`g^2N`$ fixed ) suggests that in the pure gauge theory the ground-state energy depends on $`\theta `$ to leading order in $`1/N`$ . Numerical nonperturbative studies, based on Monte Carlo simulations of the lattice formulation of SU($`N`$) gauge theories (see e.g. Ref.  for recent reviews), have given nonzero estimates for the topological susceptibility $`\chi `$ at $`\theta =0`$, which is the zero-momentum two-point correlation function of $`q(x)`$, and therefore the second derivative of the ground state energy with respect to $`\theta `$ computed at $`\theta =0`$.

The $`\theta `$ dependence of the theory is described by the partition function

$$Z(\theta )=[dA]\mathrm{exp}\left(d^4x_\theta \right)\mathrm{exp}[VF(\theta )]$$
(3)
where $`F(\theta )`$ is the free energy (ground state energy) per unit space-time volume. Some information on the behavior of $`F(\theta )`$ comes from large-$`N`$ arguments. Witten argued  that in the large-$`N`$ limit $`F(\theta )`$ is a multibranched function of the type

$$F(\theta )=N^2\mathrm{min}_kH\left(\frac{\theta +2\pi k}{N}\right)$$
(4)
which is periodic in $`\theta `$, but not smooth since at some value of $`\theta `$ there is a jump between two different branches. More recently, the conjecture was refined  leading to a rather simple expression for $`\mathrm{\Delta }F(\theta )F(\theta )F(0)`$ in the large-$`N`$ limit, that is

$$\mathrm{\Delta }F(\theta )=𝒜\mathrm{min}_k(\theta +2\pi k)^2+O\left(1/N\right).$$
(5)
In particular, for sufficiently small values of $`\theta `$, i.e. $`\theta <\pi `$,

$$\mathrm{\Delta }F(\theta )=𝒜\theta ^2+O\left(1/N\right).$$
(6)
Thus possible $`O(\theta ^4)`$ terms are expected to be depressed by powers of $`1/N`$. This conjecture has been supported using arguments based on a duality of large-$`N`$ gauge theories with string theory . It has also been discussed in a field-theoretical framework in Ref. .

In this work we study the $`\theta `$ dependence of SU($`N`$) gauge theories, using numerical simulations of their Wilson lattice formulation. Numerical Monte Carlo studies of the $`\theta `$-dependence are made very difficult by the complex nature of the $`\theta `$ term. In fact the lattice action corresponding to the Lagrangian (1) cannot be directly simulated for $`\theta 0`$, not to mention the difficulty to define a topological charge density on the lattice. Here we restrict ourselves to the region of small $`\theta `$ values, where one may expand $`F(\theta )`$ around $`\theta =0`$

$$\mathrm{\Delta }F(\theta )=\underset{i=1}{}A_{2i}\theta ^{2i}.$$
(7)
The coefficient $`A_2`$ gives the topological susceptibility, i.e.

$$A_2=\frac{1}{2}\chi ,\chi \chi _2=\frac{1}{V}Q^2_{\theta =0},$$
(8)
where $`Q`$ is the topological charge $`Q=d^4xq(x)`$. Higher order coefficients of the expansion (7) can be related to higher moments of the probability distribution $`p(Q)`$ of the topological charge. For instance,

$`A_4={\displaystyle \frac{1}{24}}\chi _4,`$
$`\chi _4={\displaystyle \frac{1}{V}}\left[Q^4_{\theta =0}3\left(Q^2_{\theta =0}\right)^2\right].`$ (9)
A convenient parametrization of the free energy is achieved by introducing the dimensionless function

$$f(\theta )=\sigma ^2\mathrm{\Delta }F(\theta ),$$
(10)
where $`\sigma `$ is the string tension at $`\theta =0`$. The coefficients of the $`\theta `$-expansion of $`f(\theta )`$ are dimensionless scaling quantities. Let us write $`f(\theta )`$ in the form

$$f(\theta )=\frac{1}{2}C\theta ^2s(\theta )$$
(11)
where

$`C={\displaystyle \frac{2A_2}{\sigma ^2}}={\displaystyle \frac{\chi }{\sigma ^2}},`$ (12)
$`s(\theta )=1+b_2\theta ^2+b_4\theta ^4+\mathrm{},`$ (13)
$`b_2={\displaystyle \frac{A_4}{A_2}}={\displaystyle \frac{\chi _4}{12\chi }},`$ (14)
The function $`s(\theta )`$ parametrizes the corrections to a simple Gaussian behavior. Witten’s conjecture implies that $`C`$ has a finite nonzero large-$`N`$ limit $`C_{\mathrm{}}`$ with $`O(1/N^2)`$ corrections, and $`b_{2i}=O(1/N^2)`$. In particular one expects that $`b_2`$ is small, and that it should rapidly decrease with increasing $`N`$. We recall that the nonzero value of $`\chi _{\mathrm{}}`$ is essential to provide an explanation to the U(1)A problem, and can be related to the $`\eta ^{}`$ mass  through the relation

$$\chi _{\mathrm{}}=2𝒜=\frac{f_\pi ^2m_\eta ^{}^2}{4N_f}+O(1/N).$$
(15)
The quantity $`b_2`$ also lends itself to a physical interpretation, being related to the $`\eta ^{}\eta ^{}`$ elastic scattering amplitude .

We present results for four-dimensional SU($`N`$) gauge theories with $`N=3,4,6`$. In particular, we determine the ratio $`C=\chi /\sigma ^2`$ and the coefficient $`b_2`$ of the $`O(\theta ^2)`$ term in the function $`s(\theta )`$. Once verified that the results for $`N=3,4,6`$ are consistent with the expected behavior of the type $`A+B/N^2`$, the large-$`N`$ limit is determined by a fit. This allows us to obtain a rather accurate estimate of $`C_{\mathrm{}}`$. The estimated values of $`b_2`$ turn out to be very small, of the order of a few per cent, and appear to decrease rapidly with $`N`$, supporting Witten’s conjecture. Our results provide a first quantitative study of the function $`s(\theta )`$, cf. Eq. (11), related to the corrections to the $`\theta ^2`$ term in the free energy $`F(\theta )`$, i.e. to a simple Gaussian behavior.

The paper is organized as follows. In Sec. 2 we describe the Monte Carlo simulations that we have performed. In Sec. 3 we present and discuss the results.

## 2 Numerical simulations

We consider the Wilson formulation of lattice gauge theories:

$$S=N\beta \underset{x,\mu >\nu }{}\mathrm{Tr}[U_\mu (x)U_\nu (x+\mu )U_\mu ^{}(x+\nu )U_\nu ^{}(x)+\mathrm{h}.\mathrm{c}.],$$
(16)
where $`U_\mu (x)`$ SU($`N`$) are link variables. In our simulations we employed the Cabibbo-Marinari algorithm  to upgrade SU($`N`$) matrices by updating their SU(2) subgroups (we selected $`N(N1)/2`$ subgroups and each matrix upgrading consists of $`N(N1)/2`$ SU(2) updatings). This was done by alternating microcanonical over-relaxation and heat-bath steps, typically in a 4:1 ratio. Tables 1, 2 and 3 contain some information on our MC runs respectively for $`N=3,4,6`$: The coupling values111 In order to compare Monte Carlo results for various values of $`N`$, it is useful to introduce the rescaled coupling $`\gamma 1/(g_0^2N)=\beta /(2N^2)`$, which is kept fixed when taking the large-$`N`$ limit of the lattice theory. $`\beta `$ or $`\gamma \beta /(2N^2)`$, the corresponding value of the string tension (taken from Refs. ), lattice sizes and the number of sweeps. The values of $`\beta ,\gamma `$ were chosen to lie in the weak-coupling region, i.e., beyond the first order phase transition in the case $`N=6`$, and beyond the crossover region characterized by a peak of the specific heat for $`N=3,4`$; see Ref.  for a more detailed discussion of this point. Computing quantities related to topology using lattice simulation techniques is not a simple task. In a lattice theory the fields are defined on a discretized set, therefore the topological properties are strictly trivial. One relies on the fact that the physical topological properties are recovered in the continuum limit. Various techniques have been proposed and employed to associate a topological charge to a lattice configuration, see e.g. the review . Here we employ a cooling technique. The topological charge has been measured on cooled configurations, using the standard twisted double plaquette operator. As is well known, the sum over the whole lattice of the twisted double plaquette, measured on cooled configurations, takes on values $`Q_tk\alpha `$, where $`k`$ is an integer and $`\alpha 1`$. We determine the typical value of $`\alpha `$ by minimizing the average deviation of the twisted double plaquette from integer multiples of $`\alpha `$, and then assign to $`Q`$ the integer closest to $`Q_t/\alpha `$. This method eliminates the need for expensive, protracted cooling; usually fewer than 20 steps suffice (when using $`N(N1)/2`$ subgroups) in order to observe a substantial convergence of the results, which become indistinguishable from results obtained using either improved operators for $`Q`$ or more cooling steps. We measured the topological charge $`Q`$ typically every 100 sweeps.

As already observed in Ref. , a severe form of critical slowing down affects the measurement of $`Q`$, posing a serious limitation for numerical studies of the topological properties in the continuum limit, especially at large values of $`N`$. In Tables 1, 2 and 3 we report the estimates of the autocorrelation time $`\tau _Q`$ obtained by a blocking analysis of the data, where the quoted errors on the estimates of $`\tau _Q`$ are just indicative, and attempt to quantify the uncertainty in extracting $`\tau _Q`$ from the blocking procedure outlined in Ref. . These results are suggestive of an interesting phenomenon: As shown in Fig. 1, the data turn out to be well reproduced by an exponential behavior of the type

$$\mathrm{ln}\tau _QN\left(e_N\xi _\sigma +c_N\right)$$
(17)
where $`\xi _\sigma `$ is the length scale associated with the string tension, i.e. $`\xi _\sigma \sigma ^{1/2}`$, and $`e_N`$, $`c_N`$ are constants. A fit to the data for $`N=3,4,6`$, using the function (17), gives respectively $`e_3=0.40(1)`$ and $`c_3=0.42(4)`$, $`e_4=0.42(1)`$ and $`c_4=0.43(4)`$, $`e_6=0.48(3)`$ and $`c_6=0.55(9)`$, showing that $`\mathrm{ln}\tau _Q`$ is approximately proportional to $`N`$ (at fixed $`\xi _\sigma `$). A qualitative explanation of this severe form of critical slowing down may be that topological modes give rise to sizeable free-energy barriers separating different regions of the configuration space. As a consequence, the evolution in the configuration space may present a long-time relaxation due to transitions between different topological charge sectors. We do not have a priori arguments in favor of the exponential behavior (17), which is just a guess arising naturally after observing the plot in Fig. 1. Note also that this dramatic effect is not observed in plaquette-plaquette or Polyakov line correlations, suggesting an approximate decoupling between topological modes and nontopological ones, such as those determining the confining properties.

## 3 Results

In this section we present and discuss the Monte Carlo results concerning the topological susceptibility $`\chi `$ and the function $`s(\theta )`$, cf. Eq. (11). In Tables 1, 2 and 3 we report $`\chi `$, the ratio $`C=\chi /\sigma ^2`$, and the coefficient $`b_2`$, cf. Eq. (14). Data for different lattice sizes allow us to check for finite size effects: their comparison shows that they are under control. The ratio $`C`$ is plotted in Fig. 2 versus $`\sigma `$, to evidentiate possible scaling corrections, which are expected to be $`O(a^2)`$ apart from logarithms. Such corrections are indeed clearly observed for $`N=3`$ and $`N=4`$, and one needs to allow for them in the fitting procedure to extrapolate to the continuum limit. For this purpose we assume a linear behavior

$$C=C_{\mathrm{cont}}+u\sigma ,$$
(18)
which takes into account of the leading $`O(a^2)`$ scaling correction and it is also supported by the data. We obtain

$`C_{\mathrm{cont}}=0.0282(12),u=0.077(24)`$ $`\mathrm{for}N=3,`$ (19)
$`C_{\mathrm{cont}}=0.0257(10),u=0.049(15)`$ $`\mathrm{for}N=4.`$ (20)
For $`N=6`$ no evidence of scaling corrections is observed. Therefore, in this case as final estimate we take the weighted average of the data for the largest values of $`\beta ,\gamma `$, i.e. for $`\gamma =0.348`$ and lattice size $`L=12,16`$ and for $`\gamma =0.350`$, obtaining

$$C_{\mathrm{cont}}=0.0236(10)\mathrm{for}N=6.$$
(21)
Results for smaller $`\gamma `$ are used to control scaling corrections: they are perfectly consistent, indicating that scaling corrections are at most of the size of the error 222 We mention that a fit of the $`N=6`$ data to a constant gives $`C_{\mathrm{cont}}=0.0232(5)`$ with $`\chi ^2/d.o.f.0.5`$, while a linear fit to $`a+b\sigma `$ gives $`C_{\mathrm{cont}}=0.0240(25)`$ and $`u=0.01(3)`$..

These results of $`C_{\mathrm{cont}}`$ are plotted in Fig. 3 versus $`1/N^2`$, which is the expected order of the asymptotic corrections in the large-$`N`$ limit. They are clearly consistent with a behavior of the type

$$C_{\mathrm{cont}}(N)=C_{\mathrm{}}+\frac{B}{N^2}.$$
(22)
Assuming such a behavior, a fit to the results for $`N=3,4,6`$ gives

$$C_{\mathrm{}}=0.0221(14),B=0.055(18),$$
(23)
which provides a rather accurate estimate of the large-$`N`$ limit of $`C`$; this limit is clearly seen to be nonzero. Note that, extrapolating using only the $`N=4,6`$ results, we obtain a perfectly consistent result, i.e. $`C_{\mathrm{cont}}=0.0219(20)`$ and $`B=0.06(4)`$.

In order to compare with the numbers appearing in the literature, we translate the above results in physical units assuming the rather standard value $`\sqrt{\sigma }=440\mathrm{MeV}`$, for $`N3`$. We obtain $`\chi _{\mathrm{}}^{1/4}=170(3)\mathrm{MeV}`$, and $`\chi ^{1/4}=180(2)\mathrm{MeV}`$ for $`N=3`$. These numbers have to be compared with the r.h.s. of Eq. (15), which, using the actual values of $`f_\pi `$, $`m_\eta ^{}`$, and $`N_f=3`$, gives $`(191\mathrm{MeV})^4`$; using the formula refined by Veneziano  for which $`m_\eta ^{}^2m_\eta ^{}^2+m_\eta ^22m_K^2`$, one obtains $`(180\mathrm{MeV})^4`$. Our results for $`N=3`$ are in substantial agreement with the estimates obtained so far using various methods, see, e.g., Refs. . Ref.  presents an estimate of $`C_{\mathrm{}}`$ obtained by extrapolating results for $`N=2,3,4,5`$ using the formula (22); they obtain $`C_{\mathrm{}}=0.0200(43)`$ from $`C_{\mathrm{cont}}(N=2)=0.0545(25)`$, $`C_{\mathrm{cont}}(N=3)=0.0355(33)`$, $`C_{\mathrm{cont}}(N=4)=0.0224(39)`$ and $`C_{\mathrm{cont}}(N=5)=0.0224(49)`$. Our work improves these results essentially for two reasons: we performed simulations with much higher statistics, and we do not rely on $`N=2`$ data for our extrapolation, which may be too small a value of $`N`$. The large-$`N`$ extrapolation of Ref.  is in agreement with ours, although we note some differences in the finite $`N`$ results. The extrapolation of Ref.  strongly depends on the $`N=2`$ result, while their $`N=2`$ result is not fitted by our extrapolation (23), which would give approximately 0.036 for $`N=2`$. Thus, at this stage the agreement of the large-$`N`$ extrapolations appears to be rather accidental. The result for $`N=2`$ of Ref.  would suggest that $`N=2`$ is outside the region where the behavior (22) sets in.

Let us now turn to the coefficient $`b_2`$ in the expansion (13) of the function $`s(\theta )`$, which determines the correction of $`F(\theta )`$ from the simple $`\theta ^2`$ behavior. Note that a reasonably accurate computation of $`\chi _4`$, which is necessary to determine $`b_2`$, requires very high statistics due to a large cancellation between the two terms in its definition (9). Sufficiently accurate results have been obtained only in the cases $`N=3,4`$. The results are listed in Tables 1 and 2 for $`N=3,4`$ respectively. For $`N=6`$ they were always consistent with zero but with a large error, essentially due to the much smaller statistics available, taking also into account the critical slowing down phenomenon discussed before. We report only the most precise result obtained by a high-statistics simulation: $`b_2=0.01(2)`$ for $`\gamma =0.348`$ and $`L=12`$. The results for $`N=3`$ and $`N=4`$ are plotted in Fig. 4. They show scaling within the errors. Assuming that the scaling corrections are negligible with respect to the typical errors of the data, as final estimate we consider the weighted average of the data, and as error we quote their typical error (in this computation we discard the data for $`\beta =6.0`$ and $`L=12`$ because of possible finite size effects). We obtain 333 Fitting the data to a constant, one obtains $`b_2=0.0226(34)`$ for $`N=3`$ and $`b_2=0.0125(36)`$ for $`N=4`$, with acceptable $`\chi ^2`$. These errors are substantially smaller than those shown in Eqs. (24) and (25). We believe that the final errors reported in the text are more reliable, since scaling corrections appear to be negligible only within the typical errors of the data. Nevertheless, we note that the $`1/N^2`$ behavior is verified even within the smaller errors of the fits.

$`b_2=0.023(7),`$ $`\mathrm{for}N=3,`$ (24)
$`b_2=0.013(7),`$ $`\mathrm{for}N=4.`$ (25)
For $`N=6`$ we report the estimate $`b_2=0.01(2)`$ obtained at $`\gamma =0.348`$. These results show that $`b_2`$ is very small and are consistent with the hypothesis $`b_2=O(1/N^2)`$. One may also assume $`b_2b_{2,2}/N^2`$ and fit the above results, obtaining $`b_{2,2}=0.21(5)`$.

As already mentioned, we used a cooling technique to determine the topological charge of a lattice configuration. Several numerical studies have shown that cooling methods provide quite reliable results. Nevertheless, some of the arguments on which cooling relies are not rigorous. Therefore it would be worthwhile to confirm our results using alternative methods. A promising and more rigorous approach is the one based on the index theorem for Ginsparg-Wilson fermions , using, e.g., the overlap lattice fermion formulation .

In conclusion, we have computed the $`\theta `$ dependence of the ground-state energy of four-dimensional SU($`N`$) gauge theories around $`\theta =0`$, by evaluating numerically the coefficients of its expansion to $`O(\theta ^4)`$. We have verified that the coefficient of the $`\theta ^2`$ term, which is the topological susceptibility, has a nonzero large-$`N`$ limit, in agreement with the picture proposed by Witten  to provide an explanation to the U(1)A problem. Corrections to the leading $`\theta ^2`$ term turn out to be very small. Indeed, the coefficient $`b_2`$ of the $`\theta ^4`$ term is a few per cent for $`N=3`$, and appears to decrease with increasing $`N`$, consistently with Witten’s conjecture  implying $`b_2=O(1/N^2)`$. Clearly, higher-order terms in the expansion in powers of $`\theta ^2`$ might not be suppressed by powers of $`N`$ and therefore might spoil the simple $`\theta ^2`$ behaviour predicted by Witten for large $`N`$. The study of these terms requires the measurement of higher moments of the topological charge probability distribution and are beyond the scope of this work. Thus, for $`N3`$ the simple Gaussian form

$$f(\theta )\frac{1}{2}C\theta ^2$$
(26)
is expected to provide a good approximation of the dependence on $`\theta `$ of the scaling free energy, up to $`O(\theta ^6)`$ terms, and therefore for a relatively large range of values of $`\theta `$.

###### Acknowledgments.

We thank Maurizio Davini for his valuable and indispensable technical support. H. P. would like to thank the Theory Group in Pisa for their hospitality during various stages of this work. We thank G. Veneziano for helpful comments.