text:
# A tomographic approach to quantum nonlocality

## 1 Introduction

In their famous paper of 1935 , Einstein, Podolsky, and Rosen (EPR) questioned the completeness of quantum mechanics (QM). Their argument was based on the premise of no action at a distance (locality) and realism. Several years later Bell  showed that the prediction of QM are incompatible with the premises of local realism (or local hidden variable theories).

Experiments  based on Bell’s result support QM, indicating the failure of local hidden variable theories. Quantitative tests used mainly Bohm’s (dichotomic) version  of the EPR entangled states instead of the original EPR states with continuous degrees of freedom.

In recent years, systems with continuous variable has attracted much atttention in connection with the burgeoning field of quantum information . In such a field EPR entanglement and quantum nonlocality are of practical importance. Nevertheless, the generalization of Bell’s inequalities to quantum systems with continuous variabbles still remains a challenging issue.

In continuous variable bipartite quantum systems, EPR aspects can arise when trying to infer quadratures of one subsystem from those of the other . These can be tested by exploiting the Heisenberg uncertainty principle . Such approach has been employed in Refs.. Few other theoretical proposals use the field quadrature measurements , while a phase space approach to quantum nonlocality has been developed  and tested .

On the other hand quantum tomography  provides a useful tool to reach a complete state knowledge. Such knowledge could then be used for nonlocal tests. Some aspects of EPR problem have been considered in tomographic approach in .

Here, we propose to use the tomographic data to do nonlocal tests. At first instance, we shall develop an inequality based on the tomograms. However, this will result quite difficult to violate. A state does not have to violate all possible Bell’s inequalities to be considered quantum nonlocal; rather, a given state is nonlocal when it violates any Bell’s inequality . Thus the degree of quantum nonlocality that we can uncover crucially depends not only on the given quantum state, but also on the Bell operator . Then, we shall introduce pseudospin operators whose statistics can be inferred from the data characterizing the reconstructed state, thus giving the possibility to use standard Bell’s inequalities. We shall also present some illuminating examples.

## 2 Quantum tomography revisited

A quantum state is described by the density operator $`\rho `$ or by any its phase space representation like the Wigner function $`W(q,p)`$, where $`q`$, $`p`$ are canonical conjugate variables. Quantum tomography is a technique which allows to recover the quantum state from a set of measured probability distributions. The latter, sometimes called tomograms, are line projections of the Wigner function. Such projections can be parametrized through a symplectic transform . Practically, the tomograms can be expressed as

$$w(X,\mu ,\nu )=W(q,p)\delta (Xq\mu p\nu )𝑑q𝑑p,$$
(1)
where $`\mu ,\nu 𝐑`$, are the parameters of the symplectic transformation and $`X`$ is a stochastic variable. It represents the random outcomes of the measurement of the observable $`X=\mu q+\nu p`$, and $`w`$ results the probability distribution associated to such observable. It is worth noting that the above symplectic transformation, once thought as composition of rotations and squeezing, can be realized in optical systems . A remarkable property of the tomograms (1) is their homogeneity.

To get the phase space picture of the quantum state one has to invert the relation (1), obtaining

$$W(q,p)=w(X,\mu ,\nu )e^{i\mu qi\nu p+iX}𝑑X𝑑\mu 𝑑\nu .$$
(2)
This simply results an inverse Fourier transform. Neverthless, it becomes a more involved inverse transform, namely an inverse Radon transform, when reducing the symplectic transformation to a mere rotation, i.e., $`\mu =\mathrm{cos}\theta `$, $`\nu =\mathrm{sin}\theta `$. This is the case of optical homodyne tomography (OHT).

By expanding the density operator in terms of a complete set of operators, it is also be possible to directly relate the quantum state to the tomograms, that is

$$\rho =w(X,\mu ,\nu )K(X,\mu ,\nu )𝑑X𝑑\mu 𝑑\nu ,$$
(3)
where the kernel operator takes the form

$$K(X,\mu ,\nu )=\frac{1}{2\pi }\lambda ^2\mathrm{exp}\left(i\lambda X+i\lambda ^2\frac{\mu \nu }{2}\right)\mathrm{exp}(i\lambda \mu q)\mathrm{exp}(i\lambda \nu p).$$
(4)
Here, $`\lambda `$ can be set equal to unity; this freedom reflects the overcompleteness of the information obtainable by means of all possible tomograms (1).

In OHT it is preferable to use Eq. (3) instead of Eq. (2) since the kernel in the latter case is unbounded, thus requiring the use of an artificial cutoff to directly sample the Wigner function. On the contrary, Eq. (3) allows a direct sample of the density matrix elements on the basis where $`K`$ results bounded .

## 3 Bell-like inequality for tomograms

The above arguments can be generalized to the case of a bipartite two-mode system . In such a case the tomograms will be given by

$`w(X_1,\mu _1,`$ $`\nu _1`$ $`,X_2,\mu _2,\nu _2)={\displaystyle }dq_1dp_1{\displaystyle }W(q_1,p_1,q_2,p_2)\times `$ (5)
$`\times `$ $`\delta (X_1q_1\mu _1p_1\nu _1)\delta (X_2q_2\mu _2p_2\nu _2)dq_2dp_2,`$
where $`\mu _1`$, $`\nu _1`$, $`\mu _2`$, $`\nu _2`$ are parameters. Let us restrict to the case of OHT, and suppose to perform homodyne detection in both subsystems $`1`$ and $`2`$. Then, $`\mu _1=\mathrm{cos}\theta _1`$, $`\nu _1=\mathrm{sin}\theta _1`$, $`\mu _2=\mathrm{cos}\theta _2`$, $`\nu _2=\mathrm{sin}\theta _2`$, where $`\theta _1`$, $`\theta _2`$ are the rotation angles related to the local oscillators phase. As a consequence we have $`w(X_1,\mu _1,\nu _1,X_2,\mu _2,\nu _2)`$ $`w(X_1,\theta _1,X_2,\theta _2)`$.

We now classify the results of the measurements to be $`+`$ if the quadrature result $`X`$ is greater than zero, and $``$ otherwise . Then, we can construct the following probabilities from the tomograms

$`w_{++}(\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑X_1{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}w(X_1,\theta _1,X_2,\theta _2)𝑑X_2,`$ (6)
$`w_+(\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}𝑑X_1{\displaystyle \underset{\mathrm{}}{\overset{0}{}}}w(X_1,\theta _1,X_2,\theta _2)𝑑X_2,`$ (7)
$`w_+(\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{0}{}}}𝑑X_1{\displaystyle \underset{0}{\overset{\mathrm{}}{}}}w(X_1,\theta _1,X_2,\theta _2)𝑑X_2,`$ (8)
$`w_{}(\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle \underset{\mathrm{}}{\overset{0}{}}}𝑑X_1{\displaystyle \underset{\mathrm{}}{\overset{0}{}}}w(X_1,\theta _1,X_2,\theta _2)𝑑X_2.`$ (9)
Then, the Bell’s inequalities can be written in terms of the above probabilities. In particular, the CHSH inequality  can be written as

$$B\left|E(\theta _1,\theta _2)+E(\theta _1,\theta _2^{})+E(\theta _1^{},\theta _2)E(\theta _1^{},\theta _2^{})\right|2,$$
(10)
where now

$$E(\theta _1,\theta _2)=w_{++}(\theta _1,\theta _2)w_+(\theta _1,\theta _2)w_+(\theta _1,\theta _2)+w_{}(\theta _1,\theta _2).$$
(11)
## 4 Reconstructed data and pseudo-Spin operators

As discussed in Section 2, measuring tomograms allows to statistically sample the density matrix elements at least in some basis, e.g. Fock basis. This represents a complete knowledge of the state, also for a bipartite system. Then, such information, say the density matrix elements in Fock basis, can be used to derive the statistics of any hermitian operator, even if this is not directly measurable. The price one has to pay is the large amount of data needs to collect.

Let us now see how these arguments can be fruitful applied to the problem of quantum nonlocality. We introduce the local pseudo-spin operators

$`S_x^{(j)}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(|2n_j2n+1|+|2n+1_j2n|\right),`$
$`S_y^{(j)}`$ $`=`$ $`i{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}\left(|2n_j2n+1||2n+1_j2n|\right),`$ (12)
$`S_z^{(j)}`$ $`=`$ $`{\displaystyle \underset{n=0}{\overset{\mathrm{}}{}}}()^n|n_jn|,`$
where $`j=1,2`$ and $`|n_j`$ are Fock states. The operators (4) obey the commutation relation of spin-$`\frac{1}{2}`$ algebra, namely

$$[S_x^{(j)},S_y^{(j^{})}]=2i\delta _{j,j^{}}\epsilon _{xyz}S_z^{(j)},$$
(13)
with $`\epsilon _{xyz}`$ the totally antisymmetric tensor. Spin tomography approach was considered e.g. in .

In terms of the operators (4) the CHSH inequality  reads in its standard form, that is

$$|(𝐮,𝐯)+(𝐮^{},𝐯)+(𝐮,𝐯^{})(𝐮^{},𝐯^{})|2,$$
(14)
with

$$(𝐮,𝐯)=\left(𝐮𝐒^{(1)}\right)\left(𝐯𝐒^{(2)}\right).$$
(15)
Here $`𝐮`$, $`𝐮^{}`$, $`𝐯`$, $`𝐯^{}`$ are unit vectors in $`𝐑^3`$, while $`𝐒^{(j)}(S_x^{(j)},S_y^{(j)},S_z^{(j)})`$ and the dot indicates the ordinary scalar product in $`𝐑^3`$. Furthermore, the angle brackets in Eq. (15) denote the average over the density matrix elements assumed available from tomographic reconstruction.

Let us now study the possible violations of inequality (14) for the two preceeding examples. For simplicity, we consider the directions $`𝐮`$, $`𝐮^{}`$, $`𝐯`$, $`𝐯^{}`$ coplanar, in the plane $`xz`$, so that

$$(𝐮,𝐯)=\left(S_z^{(1)}\mathrm{cos}\theta _u+S_x^{(1)}\mathrm{sin}\theta _u\right)\left(S_z^{(2)}\mathrm{cos}\theta _v+S_x^{(2)}\mathrm{sin}\theta _v\right),$$
(16)
where $`\theta _u`$ ($`\theta _v`$) is the relative angle betweeen $`𝐮`$ ($`𝐯`$) and $`𝐳`$.

## 5 Examples

Let us going to study the possible violations of the inequalities (10) and (15) in three paradigmatic cases.

### 5.1 Example A

We first consider the two-mode squeezed vacuum state which can be generated in the non-degenerate optical parametric amplifier

$$|\psi =\sqrt{1\lambda ^2}\underset{n=0}{\overset{\mathrm{}}{}}\lambda ^n|n_1|n_2,$$
(17)
where $`\lambda =\mathrm{tanh}s`$ and $`s`$ is the squeezing parameter. Such state can also be described by the following Wigner function

$`W(q_1,p_1,q_2,p_2)`$ $`=`$ $`{\displaystyle \frac{4}{\pi ^2}}\mathrm{exp}(e^{2s}[(q_1q_2)^2+(p_1+p_2)^2].`$ (18)
$``$ $`.e^{2s}[(q_1+q_2)^2+(p_1p_2)^2]).`$
Notice that for $`s\mathrm{}`$ the state (18) becomes precisely the EPR state.

Inserting Eq. (18) into Eq. (5) we obtain

$$w(X_1,\theta _1,X_2,\theta _2)=\frac{2}{\pi }N\mathrm{exp}\left(2aX_1^22aX_2^24aX_1X_2\right),$$
(19)
where

$`a`$ $`=`$ $`{\displaystyle \frac{\mathrm{cosh}(2s)}{\mathrm{cosh}^2(2s)\mathrm{sinh}^2(2s)\mathrm{cos}^2(\theta _1+\theta _2)}},`$
$`b`$ $`=`$ $`{\displaystyle \frac{\mathrm{sinh}(2s)\mathrm{cos}(\theta _1+\theta _2)}{\mathrm{cosh}^2(2s)\mathrm{sinh}^2(2s)\mathrm{cos}^2(\theta _1+\theta _2)}},`$ (20)
$`N`$ $`=`$ $`\sqrt{a^2b^2}.`$
In this case, the tomograms depend only on the sum of the parameters.

Using the position of Eqs. (6) we obtain

$`w_{++}(\theta _1,\theta _2)`$ $`=`$ $`w_{}(\theta _1,\theta _2)={\displaystyle \frac{1}{4\pi }}\left[\pi 2\mathrm{arctan}\left({\displaystyle \frac{b}{N}}\right)\right],`$ (21)
$`w_+(\theta _1,\theta _2)`$ $`=`$ $`w_+(\theta _1,\theta _2)={\displaystyle \frac{1}{4\pi }}\left[\pi +2\mathrm{arctan}\left({\displaystyle \frac{b}{N}}\right)\right],`$ (22)
where the signs in front of $`\pi `$ refers to $`\pi /2(\theta _1+\theta _2)\pi /2`$, and should be adjusted accordingly to the range of $`(\theta _1+\theta _2)`$.

In Fig. 1(a) we show the behavior of the functions (21)-(22) vs $`\theta _1+\theta _2`$ for different values of $`\lambda `$. As we can see they never exceed the value $`1/2`$, thus the quantities $`E`$, in Eq. (11), will be always confined between $`1`$ and $`1`$. In turns, the inequality (10) will be never violated. This results seems in agreement with Bell’s conclusions stating that the original EPR state does not exhibit nonlocality because its Wigner function is positive everywhere, and as such allows a hidden variable description . However, we shall see that is simply a matter of the type of measurements.

For the state (17), Eq. (16) becomes

$$(𝐮,𝐯)=\mathrm{cos}\theta _u\mathrm{cos}\theta _v+\frac{2\lambda }{1+\lambda ^2}\mathrm{sin}\theta _u\mathrm{sin}\theta _v.$$
(23)
Then, in Fig. 1(b) we show $``$ for different values of $`\lambda `$. In this case the EPR state show its nonlocality. Indeed a maximal violation of inequality (14) is achieved for $`\lambda 1`$ ($`s\mathrm{}`$).

### 5.2 Example B

As second example we consider the entangled superposition

$$|\psi =\frac{1}{\sqrt{2}}\left(|0_1|0_2+|n_1|n_2\right).$$
(24)
Since $`|0`$ and $`|n`$ are orthogonal for $`n0`$, the state (24) resembles a spin triplet state. The corresponding Wigner function is

$`W(q_1,p_1,q_2,p_2)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi ^2}}(1+{\displaystyle \frac{2^n}{n!}}(q_1ip_1)^n(q_2ip_2)^n+.`$
$`+`$ $`{\displaystyle \frac{2^n}{n!}}(q_1+ip_1)^n(q_2+ip_2)^n+`$
$`+`$ $`.L_n(2q_1^2+2p_1^2)L_n(2q_2^2+2p_2^2))e^{q_1^2q_2^2p_1^2p_2^2}.`$
where $`L_n`$ indicates the Laguerre polynomial of order $`n`$. By means of Eqs. (5) and (5.2) we get the tomogram of the state under consideration

$`w(X_1,\theta _1,X_2,\theta _2)`$ $`=`$ $`{\displaystyle \frac{1}{2\pi }}(1+{\displaystyle \frac{H_n(X_1)H_n(X_2)}{2^{n1}n!}}\mathrm{cos}n(\theta _1\theta _2)+.`$ (26)
$`+`$ $`.{\displaystyle \frac{H_n^2(X_1)H_n^2(X_2)}{2^{2n}(n!)^2}})e^{X_1^2X_2^2}.`$
where $`H_n`$ indicates the Hermite polynomial of order $`n`$. Using the position of Eqs. (6) we obtain probabilities $`w_{\pm \pm }(\theta _1,\theta _2)`$

$`w_{++}(\theta _1,\theta _2)`$ $`=`$ $`w_{}(\theta _1,\theta _2)={\displaystyle \frac{1}{4}}+{\displaystyle \frac{H_{n1}^2(0)}{\pi 2^nn!}}\mathrm{cos}n(\theta _1+\theta _2),`$ (27)
$`w_+(\theta _1,\theta _2)`$ $`=`$ $`w_+(\theta _1,\theta _2)={\displaystyle \frac{1}{4}}{\displaystyle \frac{H_{n1}^2(0)}{\pi 2^nn!}}\mathrm{cos}n(\theta _1+\theta _2).`$ (28)
In Fig. 2(a) we show the behavior of the $`w_{++}`$ vs $`\theta _1+\theta _2`$ for different values of $`n`$. Also in this case they never exceed the value $`1/2`$, thus the quantities $`E`$, in Eq. (11), will be always confined between $`1`$ and $`1`$, and the Bell’s inequality (10) will be satisfied.

For the state (24), Eq. (16) becomes

$$(𝐮,𝐯)=\{\begin{array}{cc}\mathrm{cos}(\theta _u\theta _v)\hfill & n=1,\hfill \\ \mathrm{cos}\theta _u\mathrm{cos}\theta _v\hfill & n>1.\hfill \end{array}$$
(29)
Then, in Fig. 2(b) we show $``$ for the case of $`n=1`$. The nonclassical character of the entangled state (24) becomes now manifest.

### 5.3 Example C

As a third, and last, example we consider a pair-coherent state , namely

$$|\psi =N\underset{0}{\overset{2\pi }{}}|re^{i\phi }|re^{i\phi }𝑑\phi ,$$
(30)
where norm $`N`$ is

$$N=\frac{e^{r^2}}{2\pi \sqrt{I_0(2r^2)}}$$
(31)
with $`I_n`$ the modified Besssel function of order $`n`$. The state (30) represents two coherent states, of the same amplitude $`r𝐑`$, having a well defined phase relation although their own phase is random.

The Wigner function corresponding to the state (30) is

$`W(q_1,p_1,q_2,p_2)={\displaystyle \frac{N^2}{\pi ^2}}e^{2r^2}\mathrm{exp}(q_1^2p_1^2q_2^2p_2^2)\times `$ (32)
$`\times `$ $`{\displaystyle \underset{0}{\overset{2\pi }{}}}{\displaystyle \underset{0}{\overset{2\pi }{}}}\mathrm{exp}(\sqrt{2}r(q_1ip_1)(e^{i\phi }+e^{i\phi ^{}})+\sqrt{2}r(q_2ip_2)(e^{i\phi }+e^{i\phi ^{}}))\times `$
$`\times `$ $`\mathrm{exp}\left(r^2[2\mathrm{cos}(\phi \phi ^{})\mathrm{cos}(2\phi )\mathrm{cos}(2\phi ^{})]\right)d\phi d\phi ^{}.`$
The tomogram of the state (30) is

$$w(X_1,\theta _1,X_2,\theta _2)=\frac{N^2}{\pi }e^{2r^2}|I(X_1,\theta _1,X_2,\theta _2)|^2e^{X_1^2X_2^2},$$
(33)
where the integral $`I(X_1,\theta _1,X_2,\theta _2)`$ reads

$`I(X_1,\theta _1,X_2,\theta _2)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{2\pi }{}}}\mathrm{exp}[{\displaystyle \frac{r^2}{2}}(e^{2i(\phi \theta _1)}+e^{2i(\phi +\theta _2)})+.`$ (34)
$`+`$ $`.\sqrt{2}r(X_1e^{i(\phi \theta _1)}+X_2e^{i(\phi +\theta _2)})]d\phi .`$
Using the procedure shown in the Appendix together with some algebra, we can obtain from Eq. (33) the following expressions for the probabilities $`w_{\pm \pm }(\theta _1,\theta _2)`$:

$`w_{\pm \pm }(\theta _1,\theta _2)`$ $`=`$ $`{\displaystyle \frac{N^2}{4}}e^{2r^2}{\displaystyle \underset{0}{\overset{2\pi }{}}}{\displaystyle \underset{0}{\overset{2\pi }{}}}\mathrm{exp}\left(2r^2\mathrm{cos}(\phi _1+\phi _2)\right)\times `$
$`\times `$ $`(1\mathrm{erf}\left({\displaystyle \frac{r}{\sqrt{2}}}[e^{i(\phi _1\phi _0)}+e^{i(\phi _2+\phi _0)}]\right))\times `$
$`\times `$ $`\left(1\mathrm{erf}\left({\displaystyle \frac{r}{\sqrt{2}}}\left[e^{i(\phi _1+\phi _0)}+e^{i(\phi _2\phi _0)}\right]\right)\right)d\phi _1d\phi _2,`$
where $`\mathrm{erf}`$ denotes the error function and $`\phi _0=(\theta _1+\theta _2)/2`$. Notice that the first (second) $`\pm `$ on the l.h.s. corresponds to the first (second) $``$ under the integral on the r.h.s.

The probabilities (5.3) are no longer bounded between the values $`\{0,1/2\}`$, thus, they can lead to violation of Bell’s inequalities. In particular, we show in Fig. 3(a) the dependence of $`B`$ from $`r`$ revealing violation of the CHSH inequality in a small range of the amplitude $`r`$. This result recall that of Ref.  where the state (30) was shown to violate another Bell’s inequality.

For the state (30), Eq. (16) becomes

$$(𝐮,𝐯)=\mathrm{cos}\theta _u\mathrm{cos}\theta _v+r^2\left(1\frac{J_0(2r^2)}{I_0(2r^2)}\right)\mathrm{sin}\theta _u\mathrm{sin}\theta _v,$$
(36)
where $`J_n`$ denotes the Bessel function of order $`n`$. Then, in Fig. 3(b) we show $``$ for $`r=1.05`$. The violation of CHSH inequality becomes now more pronounced.

## 6 Conclusions

In conclusions we have presented a tomographic approach to the quantum nonlocality of a bipartite continuous variable system. At first instance we have proposed a rough use of tomograms for testing violation of local hidden variable theory. In such a case one does not really need of a complete set of tomograms. However, the presence of purely quantum correlations are not easily recognizable. As matter of fact, the EPR state does not show the nonlocal character, while a pair-coherent state does.

Then, we have proposed to exploit the tomographically reconstructed state to get the statistics of pseudo-spin operators. This, practically, maps a continuous variable system into a discrete variable system making possible the use of standard Bell’s inequalities. This seeems a more powerful approach, but the price one ought to pay is the large amount of data (a complete set of tomograms) needed to reconstruct the state. In such a case the EPR state completely shows its nonlocal character. Also the pair-coherent state evidenciates much larger violations .

Our result opens possibilities for testing QM against local hidden variable theories using very efficient detection methods like homodyne detection .

Finally, the present approach could be straightforward extended to multipartite systems, e.g. making possible test of Grenberger-Horne-Zeilinger theorem  in continuous variables version . Moreover, it could be useful to characterize quantum states used in information processing .

## Acknowledgments

VIM and EVS thank the Russian Foundation for Basic Research for partial support under projects nos 00-02-16516 and 01-02-17745.

## Appendix

By refering to Eq. (34) we note that the exponent of the imaginary argument is periodic, and since the integral is taken over a period of the exponent we can shift the angle argument as $`\phi (\theta _1\theta _2)/2\phi `$ so the integral takes the form

$`I(X_1,\theta _1,X_2,\theta _2)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{2\pi }{}}}\mathrm{exp}[{\displaystyle \frac{r^2}{2}}e^{2i\phi _0}(e^{2i\phi }+e^{2i\phi })+.`$ (37)
$`+`$ $`.\sqrt{2}re^{i\phi _0}(X_1e^{i\phi }+X_2e^{i\phi })]d\phi ,`$
where we have set

$$\phi _0=\frac{\theta _1+\theta _2}{2}.$$
(38)
Then, the integral can be rewritten as

$`I(X_1,\theta _1,X_2,\theta _2)`$ $`=`$ $`{\displaystyle \underset{0}{\overset{2\pi }{}}}\mathrm{exp}[{\displaystyle \frac{(\alpha e^{i\phi })^2}{2}}+\sqrt{2}\alpha e^{i\phi }X_1]\times `$ (39)
$`\times `$ $`\mathrm{exp}\left[{\displaystyle \frac{(\alpha e^{i\phi })^2}{2}}+\sqrt{2}\alpha e^{i\phi }X_2\right]d\phi ,`$
where we have set

$$\alpha =re^{i\phi _0}.$$
(40)
Expanding each exponent in the above integral into a series of Hermite polynomials, and taking into account the equality

$$\underset{0}{\overset{2\pi }{}}e^{i(nm)\phi }𝑑\phi =2\pi \delta _{n,m},$$
(41)
we obtain the following expression

$$I(X_1,\theta _1,X_2,\theta _2)=2\pi \underset{n=0}{\overset{\mathrm{}}{}}\frac{H_n(X_1)H_n(X_2)}{2^n(n!)^2}\alpha ^{2n}.$$
(42)
Let us now consider two generic series

$$f(x_1)=\underset{n=0}{\overset{\mathrm{}}{}}f_nx_1^n\mathrm{and}g(x_2)=\underset{m=0}{\overset{\mathrm{}}{}}g_mx_2^m,$$
(43)
whose product yields

$$h(x_1,x_2)=f(x_1)g(x_2)=\underset{n,m=0}{\overset{\mathrm{}}{}}f_ng_mx_1^nx_2^m.$$
(44)
How to express the sum of diagonal elements

$$s(x_1,x_2)=\underset{n=0}{\overset{\mathrm{}}{}}f_ng_nx_1^nx_2^n$$
(45)
in terms of the $`h(x_1,x_2)`$ ? The answer is found by taking into account the equality (41), and it reads

$$s(x_1,x_2)=\frac{1}{2\pi }\underset{0}{\overset{2\pi }{}}h(x_1e^{i\phi },x_2e^{i\phi })𝑑\phi .$$
(46)
The sum in right-hand side of Eq. (42) is the sum of diagonal terms of the product of two following series

$$f_k(\alpha )=\underset{n=0}{\overset{+\mathrm{}}{}}\frac{H_n(X_k)}{n!}\left(\frac{\alpha }{\sqrt{2}}\right)^nk=1,2.$$
(47)
If we apply the formula (46) to this sum we obtain the representation (37) for it which we have started from. To obtain probabilities $`w_{++}`$ etc. we represent $`|I(X_1,\theta _1,X_2,\theta _2)|^2`$ in the following way:

$$|I|^2=\underset{0}{\overset{2\pi }{}}\underset{0}{\overset{2\pi }{}}f_1(\alpha e^{i\phi _1})f_1^{}(\alpha e^{i\phi _1})f_2(\alpha e^{i\phi _2})f_2^{}(\alpha e^{i\phi _2})𝑑\phi _1𝑑\phi _2.$$
(48)
To calculate probability $`w_{++}`$ we need to take the following integral:

$$\underset{0}{\overset{2\pi }{}}\underset{0}{\overset{2\pi }{}}|I(X_1,\theta _1,X_2,\theta _2)|^2e^{X_1^2X_2^2}𝑑\phi _1𝑑\phi _2,$$
(49)
and analogously for other probabilities. Swap integration order in these integrals we obtain the expression (5.3).

## References


format: text

text:
# Two-dimensional Induced Ferromagnetism

## Abstract

Magnetic properties of materials confined to nanometer length scales are providing important information regarding low dimensional physics. Using gadolinium based Langmuir-Blodgett films, we demonstrate that two-dimensional ferromagnetic order can be induced by applying magnetic field along the in-plane (perpendicular to growth) direction. Field dependent exchange coupling is evident in the in-plane magnetization data that exhibit absence of hysteresis loop and show reduction in field required to obtain saturation in measured moment with decreasing temperature.

The role of dimensionality in magnetic ordering has remained an active field of research since the first argument of Bloch in 1930 followed by theoretical work of Mermin and Wagner in 1966 bloch ; mermin . Recent advances in growth and characterization techniques of nanomaterials - materials confined to nanometer length scales in one, two or all three dimensions - has enabled investigation of magnetic ordering in general and ferromagnetism in particular in one-dimensional (1D) kmnat , two-dimensional (2D) 4 ; 5  and zero-dimensional (0D) 6  systems. These studies have given opportunities to refine existing theoretical understanding regarding possible existence of long-range ferromagnetic order in low-dimensional systems ahar ; 8 ; 9 ; 10 .

Pomerantz et al. 10 ; 11 , demonstrated for the first time that one can form literally two-dimensional (2D) magnets using Langmuir-Blodgett (LB) film growth technique physrep . Using LB film growth technique, one can form 2D hexagonal lattice of metallic ions and multilayer stack of these 2D lattices can be kept separated by organic chains physrep ; prb02 . In a typical LB film metal-metal distances within and between the 2D lattices are about 5 Å and 50 Å respectively, making this arrangement ideally suited to form stacks of 2D lattices of single-ion magnets. Although several magnetic ions like manganese lbmag , iron lbiron , cobalt lbco  have been used to form LB films, clear signature of 2D ferromagnetic ordering has not been observed in these systems. Recently, signature of magnetic ordering was observed in gadolinium based LB films at unusually high temperature bohr . But due to doubtful stoichiometry and absence of systematic low temperature magnetic measurements, the nature of ordering could not be established. The magnetic properties of solids based on gadolinium, a lanthanide metal, are primarily determined by the localized 4f moments. One can observe long-range magnetic ordering here provided the exchange coupling is mediated by the hybridized 6s and 5d conduction electrons free ; maiti . On the other hand, one expects to observe paramagnetism henry  in gadolinium compounds due to absence of conduction electrons. However, in a recent systematic angle resolved photoemission measurements of oxygen- induced magnetic surface states of lanthanide metals, it was shown that gadolinium forms GdO instead of nonmetallic sesquioxide Gd2O3 schubler . The remaining one valence electron of ($`5d6s^2`$) hybridized state was found to be responsible for mediating exchange coupling to form magnetic ordering. The temperature dependence of the energy splitting of the oxygen-induced states has confirmed schubler  this possibility.

Here we present evidence of 2D ferromagnetism that can be induced by external magnetic field applied along the metallic planes in gadolinium stearate (GdSt) LB film. The gadolinium ions show normal paramagnetism when the field is applied in the out - of - plane (normal to the 2D gadolinium layers) direction. Several GdSt LB films, having 9 to 101 monolayers (ML), were deposited on 1 mm thick Si(001) substrates using an alternating trough (KSV5000) from a monolayer of stearic acid on Milli-Q (Millipore) water subphase containing 5$`\times 10^4`$ M Gd3+ ions, obtained from dissolved gadolinium acetate. The surface pressure was maintained at 30 mNm-1 during deposition and the dipping speed was 5 mm min-1. The silicon substrates were cleaned and hydrophilized according to RCA cleaning procedure. Grazing incidence x-ray reflectivity and diffuse scattering measurements were performed using a rotating anode x-ray set up, to characterize the structure of deposited LB films physrep .

We have presented magnetization data of 101 ML LB film as a function of T and H measured using a vibrating sample magnetometer (Oxford Instruments). The nature of magnetic ordering was found to be independent of number of monolayers deposited but in 101 ML sample the ordering was evident in the raw data itself despite the presence of diamagnetic signal of Si(001) substrate (refer Fig.1(a)). Magnetization isotherm measurements (M vs. H at a fixed temperature) over all four quadrants including the virgin curve were carried out as a function of magnetic field up to $`\pm `$70 kOe, applied parallel (in - plane) as well as perpendicular (out - of - plane) to the film plane at several temperatures down to 2 K. All these measurements were carried out by cooling the sample from 300 K to the desired temperature of measurement under zero magnetic field. No hysteresis was found in these M vs. H curves. Field-cooled (FC) M vs. T measurements were carried out over 2 to 100 K under 500 Oe in the cooling cycle. We also found that the magnetization values scale with number of monolayers deposited.

Silicon background was subtracted from all the data consistently before performing data analysis. Fig.1(b) depicts M vs. T curves measured with 500 Oe field in two in-plane directions, obtained by rotating the film by 90, and in an out-of-plane direction. Paramagnetic behavior in all three directions is evident. However, at lower temperatures the nature of the out-of-plane data was found to be different from that of in-plane data (Fig.1(b)) indicating different spin response in the in-plane directions as we shall discuss later. In Fig. 1(c) we have shown the specular x-ray reflectivity data of 9 ML GdSt LB film. In this data the presence of both Bragg peaks and Kiessig fringes corresponding to out-of-plane metal-metal distance and total film thickness, respectively gibaud , are evident. For films with large number of layers, Bragg peaks become strong and Kiessig fringes could not be resolved (refer 51 ML data in Fig. 1(c)). The measured 9 ML data matches quite well with the calculated reflectivity obtained from a simple model electron density profile shown in the inset. The organic portion of film has electron density of 0.32 el Å-3 and dips going to the value of 0.17 el Å-3, as observed earlier gibaud . The electron density of the head region takes a value of 0.64 el Å-3 to produce the strong Bragg peaks observed in reflectivity data. This value of electron density in the metal-plane corresponds to a molecular structure where two stearic acid tails ( having 20 Å2 area) are attached to a single gadolinium ion prb02 . We have also shown the calculated reflectivity data and corresponding electron density profile in Fig.1(c) assuming that three tails are attached to single gadolinium ion. The essential difference between two density profiles being the change of electron density in the metal-plane from 0.64 el Å-3 to 0.48 el Å-3. From these curves, we conclude that out of three valence electrons in gadolinium only two electrons participate in bonding. Diffuse scattering data of these films show clearly that the 2D metal-planes are conformal in nature and have logarithmic in-plane correlation, as observed earlier basu . The interfacial roughness comes out to be around 2 Å.

Figure 2(a) shows out-of-plane magnetization data measured at 5K, 10K and 20K temperatures. All the magnetization data plotted against (H / T) collapses to a single curve as expected for paramagnetism or superparamagnetism yusuf . The data were fitted with the expression M = MsBs(g$`\mu _B`$SH/kBT) where Ms (= Ng$`\mu _B`$S/V) is the saturation magnetization and Bs is the Brillouin function defined as ahar

$$B_s(x)=\frac{2S+1}{2S}coth\left(\frac{(2S+1)x}{2S}\right)\frac{1}{2S}coth\left(\frac{x}{2S}\right)$$
(1)
We obtained the value of spin S as 2.75 instead of expected 3.5 of the 4f moment for gadolinium. Ms (=Ng$`\mu _B`$S/V) was found to be 1.29 $`\times 10^5`$emu/mm2. This value corresponds well with the number density of gadolinium 2.53 $`\times 10^6`$ mm-2. Here, each gadolinium ion was assumed to be at the center of a cylinder of 20 Å2 cross section and length 50 Å , as obtained from fitting of the specular reflectivity data shown in Fig.1(c).

In Fig. 2(b) in-plane magnetization data taken at temperatures of 2K, 5K, 10K and 20K are shown. It is interesting to note that the slope of the curves as well as the respective saturation magnetization ($`M_s`$) values decreases as the temperature is increased. As a result, the in-plane magnetization curves do not collapse to a single curve like the out-of-plane data ruling out the existence of normal paramagnetism or superparamagnetism in these 2D planes.

However, like in out-of-plane data, no hysteresis (i.e., zero remanent magnetization and zero coercive field) was observed here. This type of field-induced ferromagnetism has been observed earlier guertin . It should be noted here that field values (Hs) at which saturation of magnetization sets in was found to decrease with decreasing temperature and we obtained values of 10.2, 21.9, 34.6 and 57.2 kOe for sample temperatures of 2, 5, 10 and 20K, respectively. The saturation magnetization was found to exhibit exponential dependence with temperature ($`logM_s=0.660.034T`$) (refer inset of Fig.2(b)). Here Ms is expressed in a unit of $`\mu _B`$/Gd and the projected value of Ms at 0K comes out to be 4.57 $`\mu _B`$/Gd. This value is less than the value of 5.5 obtained from the fitting of out-of-plane data. It is interesting to note that the magnetization curves taken in the out-of-plane direction do not show this saturation behavior, which implies strong anisotropy in GdSt LB film. In fact, effect of strong anisotropy in ferromagnetism and possible existence of even direction dependent Curie temperature has been discussed earlier 26 . We have used an anisotropic exchange to explain the induced ferromagnetism in the in-plane direction – such an exchange was noted to be sufficient to stabilize ferromagnetism in 2D systems mermin . In this formalism it is expected that increase in applied magnetic field in a particular in-plane (xy) direction will increase the effective exchange field provided Jx,J$`{}_{y}{}^{}`$Jz. The existence of saturation in magnetization only in the in-plane direction corroborates this assumption.

We could analyze all the in-plane data by using Brillouin function as used in ferromagnetism ahar

$$M=M_sB_s\left(\frac{S}{k_BT}[g\mu _BH+J_{||}S_j\underset{j}{}cos\theta _{ij}]\right)$$
(2)
where J||=Jx=Jy are in-plane exchange (assumed to be symmetric on the basis of data shown in Fig.1(b)) and $`\theta _{ij}`$ is the in - plane angle between spins Si and Sj.

With increasing applied field the sum of cosines, in Eq. (2), gets maximized. This relationship can be used to write a simplified expression for in-plane magnetization

$$M=M_sB_s(\frac{SH\mu _B}{k_BT}(g+\xi M))$$
(3)
In Fig. 3 we have plotted measured M along with the fitted curves obtained by using Eq.(3) for 2, 5, 10 and 20K data. In this analysis we have used S= 3.5 and only fitting parameter was $`\xi `$, which increases with increasing temperature. The values come out to be 2.45 $`\times 10^4`$, 8.82 $`\times 10^4`$, 2.15 $`\times 10^5`$ and 8.28 $`\times 10^5`$ mm$`{}_{}{}^{2}emu^1`$ at 2, 5, 10 and 20K temperatures respectively. It should be noted here that without invoking this field dependent exchange interaction ($`\xi `$H) the data couldn’t be analyzed. We have tried fitting the data with a Brillouin function with constant exchange in place of $`\xi `$H. The best fit was obtained (refer dashed curve in Fig. 3) with $`\xi `$H =1.27 $`\times 10^6`$ kOe mm$`{}_{}{}^{2}emu^1`$ for the 5K data. Only field dependent (linear) exchange could explain the observed magnetization data (refer Fig. 3.).

It will be interesting to investigate the roles of various models of ferromagnetism for gadolinium maiti ; schubler  and 2D dipole interactions politi  to explain the observed reduction of Hs and increase of Ms as temperature decreases. Lower saturation value of magnetization (Ms) in the in-plane direction can be simply due to non-participation of few gadolinium planes, located near substrate/film and film/air interfaces, in 2D ferromagnetic ordering. However, strong asymmetry arising from lower spin value (2.75) in the out-of-plane direction and lower Ms value with expected spin of 3.5 in the in-plane direction need systematic microscopic measurements to elucidate the effect of anisotropic g-factor for spins 28 ; 29 ; halperin .

In conclusion, we have demonstrated that GdSt LB films can provide an easy-to-form 2D ferromagnetic system. We have shown that out-of-plane magnetization resembles paramagnetic behavior down to a temperature of 2K and a field of up to 70 kOe. However, field induced ferromagnetic ordering was observed in the in-plane direction and the field (Hs) required to get the saturation magnetization decreases with decreasing temperature. The attachment of two stearic acid tails with single gadolinium ion and the presence of field induced ferromagnetic state in GdSt, found here, are consistent with the reported results of recent photoemission studies schubler . The absence of remanence in the in-plane ferromagnetic ordering, observed perhaps for the first time here, provides experimental validity of Mermin-Wagner theorem which has been extended recently bruno  for the metallic two-dimensional systems.


text:
# Stability of Solutions of Hydrodynamic Equations Describing the Scaling Limit of a Massive Piston in an Ideal gas

## 1 Introduction

The time evolution of a system consisting of a piston of mass $`M`$ moving parallel to the $`x`$-axis in a cube containing non-interacting point particles of unit mass has been studied extensively \[CLS1, CLS2, CL, G, GF, GP, H, KBM, LPS, Li, PG\]. After some rescaling of space and time (by the length of the cube) the problem reduces to that of a one dimensional system with $`N_L`$ ($`N_R`$) particles in the interval $`[0,X]`$ (resp., $`[X,1])`$ where $`X(t)`$ is the position of the piston. The left (right) particles move freely between collisions with the wall at $`x=0`$ ($`x=1`$) and the piston at $`x=X(t)`$. At the walls the velocities of the particles get reversed while at $`x=X(t)`$ the outgoing velocity $`v^{}`$ is related to the incoming velocity $`v`$ by the rules of elastic collisions \[CLS2, CL, LPS\],

$$v^{}=\frac{M1}{M+1}v+\frac{2M}{M+1}V$$
(1.1)
where $`V`$ is the incoming velocity of the piston. It follows from (1.1) that $`N_L,N_R`$, as well as the total kinetic energy, $`\frac{1}{2}_{i=1}^Nv_i^2+\frac{1}{2}MV^2`$ ($`N=N_L+N_R`$) are conserved quantities. The dynamics of the system can be reduced to a billiard in a $`(2N+1)`$-dimensional domain (polyhedron), cf. \[CL\]. It was shown in \[LPS, CLS1, CLS2\], under certain quite restrictive conditions on the initial distribution of gas particles, that, in the limit $`N\mathrm{}`$, $`MN^{2/3}`$, the dynamics of the piston and the gas satisfy a closed system of Euler-type hydrodynamical equations (HE) for a time interval $`(0,\tau )`$ in which any particle had at most two collisions with the piston.

The origin of the scaling $`MN^{2/3}`$ is as follows. For $`N`$ particles with velocities of $`O(1)`$ distributed with density of $`O(1)`$ in a parallelopiped of length $`L`$ and crossectional area $`A`$ the number of particles colliding with the piston per unit (unrescaled) time and hence the pressure (from each side) is proportional to $`A`$. To ensure that, on this time scale, the acceleration of the piston stays of $`O(1)`$ as $`L`$, $`A`$ and $`NO(AL)`$ grows to infinity it is necessary to make the mass of the piston grow as $`A`$. For a cube this corresponds to $`MN^{2/3}`$. In the rescaled units the number of collisions experienced by the piston per unit time is of $`O(N)`$ independent of $`A`$ and the HE hold for general $`M`$ as long as $`MN^\alpha `$, $`\alpha (0,1)`$, i.e. when the kinetic energy of the piston becomes negligible compared to that of the gas. The time interval (in the scaled units) during which the derivation of the HE remains valid depends on $`\alpha `$ (getting larger as $`\alpha 1`$), see Remarks 3 and 4 in Sect. 4 of \[CLS2\]. For $`\alpha =2/3`$ this time interval is such that the piston suffers no more than two collisions with any gas particle. It is however not clear from the derivation to what extent those equations may actually approximate the real evolution of the system with large $`N`$, for longer times. This led us to carry out extensive computer simulations of particle systems, with $`MN^{2/3}`$, $`N`$ as large as $`27\times 10^6`$ \[CL\], and initial conditions for which the hydrodynamic equations have trivial solutions $`X(t)=0.5`$ and $`V(t)=0`$ for all $`t>0`$. We found nevertheless that for certain initial velocity distributions the trajectory of the piston diverged from these values. In particular it was observed in these simulations that the particle plus piston system, after experiencing large random fluctuations, quickly converges to a more stable regime, in which the piston and the gas undergo regular slowly damped oscillations lasting a long time. The parameters of those oscillations (the period, the amplitude, and the rate of damping) seem to depend little on many details of the initial distributions.

A possible explanation for this behavior is that the mechanical trajectory of the system with $`N`$ large but finite, being subjected to the intrinsic random fluctuations associated with the discrete nature of the gas particles, behaves as a perturbed solution of the HE, i.e. as if the initial conditions were subjected to a small random perturbation. There large fluctuations in the particle dynamics observed experimentally in \[CL\] may thus be due to the instability of the HE for the initial conditions considered. This would imply that, on the contrary, when the solution of the hydrodynamical equations is stable, then the mechanical evolution of the system should remain close to that solution for a very long time.

Here we present evidence, but no proof, for this conjecture. We first find a class of initial conditions for which the solution of the HE are stable and note that simulations with such distributions in \[CL\] indeed yielded a piston trajectory, which remained close to the solution of the HE with $`V(t)0`$ . Furthermore, we observe that the mechanical piston trajectories, which do not stay close to the solution of the HE with the prescribed initial conditions, appear to follow some perturbed oscillating solutions of the HE equations. These solutions appear to act as an “attractor” for the HE 111Due to an obvious time-reversibility of the hydrodynamical equations, there can be no attractors, strictly speaking, but saddle-type periodic solutions can very well exist.. We cannot rigorously prove the existence of attracting periodic solutions for the HE yet, but we present some approximate constructions of such solutions.

The paper is organized as follows. In Section 2 we state the hydrodynamical equations. In Section 3 we prove global stability for a class of stationary solution. In Section 4, we use a perturbative analysis to find sufficient conditions for linear instability. In Section 5 we investigate a particular family of stationary solutions which include those used in the \[CL\] simulations and show that it contains both linearly stable and unstable ones, alternating in a very intricate manner. In Section 6, an approximate construction of periodic solutions of the HE is presented.

## 2 Hydrodynamical equations

Let $`X(t)(0,1)`$ be the position of the piston at time $`t`$ and $`V(t)`$ its velocity. We denote the continuum density of the gas in $`[0,1]\times \text{I}\text{R}`$ by a function $`p(x,v,t)`$. The HE describing the time evolution of this, continuum fluid plus piston system, are as follows.

* Free motion. Inside the container the density satisfies the standard continuity equation for a noninteracting particle system without external forces:

$$\left(\frac{}{t}+v\frac{}{x}\right)p(x,v,t)=0$$
(2.1)
for all $`x(0,1)`$ except $`x=0`$, $`x=1`$ and $`x=X(t)`$.

Equation (2.1) has a simple solution

$$p(x,v,t)=p(xvs,v,ts)$$
(2.2)
for $`0<s<t`$ such that $`xvr\{0,X(tr),1\}`$ for all $`r(0,s)`$. Equation (2.2) has one advantage over (2.1): it applies to all points $`(x,v)`$, including those where the function $`p`$ is not differentiable, or even continuous.

* Collisions with the walls. At the walls $`x=0`$ and $`x=1`$ we have

$$p(0,v,t)=p(0,v,t)$$
(2.3)
$$p(1,v,t)=p(1,v,t)$$
(2.4)
* Collisions with the piston. At the piston $`x=X(t)`$ we have (this is obtained from (1.1) when $`M\mathrm{}`$)

$`p(X(t)0,v,t)`$ $`=`$ $`p(X(t)0,2V(t)v,t)\mathrm{for}v<V(t)`$
$`p(X(t)+0,v,t)`$ $`=`$ $`p(X(t)+0,2V(t)v,t)\mathrm{for}v>V(t)`$ (2.5)
where $`v`$ represents the velocity after the collision and $`2V(t)v`$ that before the collision; and

$$X(t)=X(0)+_0^tV(s)𝑑s$$
(2.6)
is the (deterministic) position of the piston.

It remains to describe the evolution of $`V(t)`$ (which we take to be left continuous). Suppose the piston’s position at time $`t`$ is $`X`$ and its velocity $`V`$. The piston is affected by the fluid at $`(x,v)`$ exerting pressure on it from the right ($`x=X+0`$ and $`v<V`$) and from the left ($`x=X0`$ and $`v>V`$). Accordingly, we define the density of the fluid in contact with the piston (“density on the piston”) by

$$q(v,t;X,V)=\{\begin{array}{cc}p(X+0,v,t)\hfill & \mathrm{if}v<V\hfill \\ p(X0,v,t)\hfill & \mathrm{if}v>V\hfill \end{array}$$
(2.7)
* Piston’s velocity. The velocity $`V=V(t)`$ of the piston satisfies the equation

$$_{\mathrm{}}^{\mathrm{}}(vV)^2\mathrm{sgn}(vV)q(v,t;X,V)𝑑v=0$$
(2.8)

The origin of eqs. (H1)-(H3) in the particle system is clear. Eq. (H4) is essentially a force balance equation—since the rate of collision of the piston with particles on either side and consequent force on piston is much larger than mass of the piston when $`N/M\mathrm{}`$, $`V`$ adjusts instantaneously to make the forces from the two sides balance exactly. The system of (hydrodynamical) equations (H1)–(H4) is now closed and, given initial conditions, satisfying (2.3)-(2.8) at $`t=0`$, completely determine the functions $`X(t)`$, $`V(t)`$ and $`p(x,v,t)`$ for $`t>0`$. When the initial conditions do not satisfy these equations one has to imagine that they become satisfied instantaneously for $`t=0+`$. The existence and uniqueness of solutions of (H1)–(H4) were proven, under general conditions, in \[CLS1, CLS2\]. We need only to assume that the $`p(x,v,0)`$ is bounded, piecewise differentiable, and either has a compact support in the $`x,v`$ plane or decays fast enough as $`|v|\mathrm{}`$. We also require that $`p(x,v,0)𝑑v>0`$ for all $`x`$.

The HE, like the Vlasov equations for plasmas, are time-reversible, see \[P\] and \[MP\]. They preserve the classical integrals of motion. The mass of the fluid to the left and to the right of the piston as well as the total kinetic energy of the fluid remain constant along any solution.

* Mass conservation

$$_L=_0^{X(t)}p(x,v,t)𝑑v𝑑x,_R=_{X(t)}^1p(x,v,t)𝑑v𝑑x$$
* Energy conservation

$$2E=v^2p(x,v,t)𝑑v𝑑x$$
(2.9)

Equation (2.8) also preserves the total momentum $`vp(x,v,t)𝑑v𝑑x`$, but the latter changes due to reflections at the walls. We note that the piston itself does not contribute to the total momentum and energy of the system in this model, because its mass and energy vanish, when divided by $`N`$, in the limit $`N\mathrm{}`$. (The mass, energy and momentum of the fluid all correspond to the original quantities in the particle system divided by $`N`$).

The HE define a dynamics on the domain $`G:=\{(x,v):0x1\}`$ in which every point $`(x,v)G`$ moves freely with constant velocity and collides elastically with the walls and the piston. Denote by $`(x_t,v_t)`$ the position and velocity of an arbitrary point at time $`t0`$. Then (H1) translates into $`\dot{x}_t=v_t`$ and $`\dot{v}_t=0`$ whenever $`x_t\{0,1,X(t)\}`$, (H2) becomes $`(x_{t+0},v_{t+0})=(x_{t0},v_{t0})`$ whenever $`x_{t0}\{0,1\}`$, and (H3) gives

$$(x_{t+0},v_{t+0})=(x_{t0},2V(t)v_{t0})$$

whenever $`x_{t0}=X(t)`$. Note that the point $`(x_t,v_t)`$ moves in $`G`$ and reflects at the walls and the piston as if those had infinite mass.

The motion of points in $`G`$ is described by a one-parameter family of transformations $`F^t:GG`$ defined by $`F^t(x_0,v_0)=(x_t,v_t)`$ for $`t>0`$. We will also write $`F^t(x_t,v_t)=(x_0,v_0)`$. According to (H1)–(H3), the density $`p(x,v,t)`$ satisfies a simple equation

$$p(x_t,v_t,t)=p(F^t(x_t,v_t),0)=p(x_0,v_0,0)$$

for all $`t0`$. It is easy to see that for each $`t>0`$ the map $`F^t`$ is one-to-one and preserves area, i.e. $`\mathrm{det}|DF^t(x,v)|=1`$. Hence, the family $`F^t`$ describes an incompressible flow on $`G`$ and consequently:

* Incompressibility. For any $`a<b`$ the Lebesgue measures (areas) of the sets

$$\{(x,v):a<p(x,v,t)<b,0<x<X(t)\}$$

and

$$\{(x,v):a<p(x,v,t)<b,X(t)<x<1\}$$

remain constant in time

A particular case in which it is possible to solve equations (H1)–(H4) analytically, is when the initial distribution is stationary. This happens when $`p(x,v,0)`$ satisfies two conditions:

* Uniformity and symmetry. The initial density $`p(x,v,0)=p(x,v)`$ is of the form

$$p(x,v)=\{\begin{array}{cc}p_L(|v|)\hfill & \mathrm{for}x<X_0\hfill \\ p_R(|v|)\hfill & \mathrm{for}x>X_0\hfill \end{array}$$

for all $`v`$ and $`X(0)=X_0`$.
* Pressure balance. The pressure on the piston from both sides is equal:

$$P_L:=2_0^{\mathrm{}}v^2p_L(v)𝑑v=P_R:=2_0^{\mathrm{}}v^2p_R(v)𝑑v$$
(2.10)

Under conditions (S1)–(S2) the equations (H1)–(H4) have a simple solution: the system remains frozen in its initial state:

$$X(t)X_0,V(t)0,p(x,v,t)p(x,v,0)$$
(2.11)
for all $`t>0`$.

We will analyze in the next three sections the stability of this stationary solution. Note that there is no requirement on the form of $`p_L(|v|)`$ or $`p_R(|v|)`$; all that is required is a balance of forces (2.10).

## 3 Globally stable solutions

Here we consider stationary solutions $`p(x,v)`$ satisfying (S1)–(S2) and an additional monotonicity requirement:

$$p_L(|v_1|)p_L(|v_2|)\mathrm{and}p_R(|v_1|)p_R(|v_2|)$$
(3.1)
for all $`|v_1||v_2|`$. We claim that such solutions are globally stable. This criteria is very similar to the stability criteria for the Vlasov equation described by Penrose \[P\] and by Marchioro and Pulvirenti \[MP\].

Before we state our result, we introduce some notation. Denote by $``$ the following special norm on the space of functions on $`G`$

$$f(x,v)g(x,v)=|f(x,v)g(x,v)|(1+v^2)𝑑v𝑑x$$
(3.2)
###### Theorem 3.1

Let $`p(x,v)`$ satisfy (S1), (S2) and (3.1). Then for any $`\epsilon >0`$ there exists $`\delta >0`$ such that if the initial density $`p(x,v,0)`$ satisfies $`p(x,v,0)p(x,v)<\delta `$, and $`X(0)=X_0`$ then

* $`p(x,v,t)p(x,v)<\epsilon `$;
* $`|X(t)X(0)|<\epsilon `$

for all $`t>0`$.

Proof. The proof of claim (i) is identical to the proof of the stability theorem of Marchioro and Pulvirenti \[MP\]. We note that their theorem is stated in the $`L_1`$ norm, but, in fact it was proven in the (3.2) norm. (There is a minor error in \[MP\], the conditions that assure that the $`L_1`$ norm and the (3.2) norm are equivalent are not explicitely assumed.)

It is clear from (2.9) that, given the position of the piston $`X`$ and values of the areas of the level sets, defined in $`(D3)`$, the minimal possible value of the total energy for any phase-space density $`\pi (x,v)`$, is attained when $`\pi (x,v)`$ is uniform in $`x`$ and monotonically decreasing in $`|v|`$ in each compartment.

Consider first the case where $`\pi (x,v)`$ has the same area of the level sets as some $`p(x,v)`$ satisfying $`(S1)`$ and $`(S2)`$. Then the minimum of the energy when the piston position is $`X`$ is attained when

$$\pi (x,v)=p_L(vX/X_0),0<x<X$$

and

$$\pi (x,v)=p_R(v(1X)/(1X_0)),X<x<1$$

The minimal total energy is then

$`2E_{\mathrm{min}}(t)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle _0^X}v^2\pi (x,v)𝑑x𝑑v+{\displaystyle _0^{\mathrm{}}}{\displaystyle _X^1}v^2\pi (x,v)𝑑x𝑑v`$
$`=`$ $`{\displaystyle \frac{X_0^3}{X^2}}{\displaystyle _0^{\mathrm{}}}v^2p_L(v)𝑑v+{\displaystyle \frac{(1X_0)^3}{(1X)^2}}{\displaystyle _0^{\mathrm{}}}v^2p_R(v)𝑑v`$
(we used a change of variable $`u=vX/X_0`$ in the first integral and $`u=v(1X)/(1X_0)`$ in the second one). Using the pressure balance (2.10) and denoting $`P=P_L=P_R`$ gives

$$E_{\mathrm{min}}(t)=\frac{P}{2}\left(\frac{X_0^3}{X^2}+\frac{(1X_0)^3}{(1X)^2}\right)$$

Consider now the above expression as a function of $`X`$. Its minimum is attained at the point where $`dE_{\mathrm{min}}/dX=0`$, i.e.

$$\frac{X_0^3}{X^3}=\frac{(1X_0)^3}{(1X)^3}$$

which is only possible if $`X=X_0`$. Therefore, the state $`X=X_0`$ provides a unique minimum of the total energy function under the incompressibility constraint (D3). This proves the theorem when the $`p(x,v,0)`$ has exactly the same area of the level sets as $`p(x,v)`$. For perturbed initial densities $`p(x,v,0)`$ the above estimates only hold approximately, hence small changes in the system are possible, but large changes are not. This proves the theorem.

Slightly extending the above theorem, let us assume that the initial density $`p(x,v,0)`$ satisfies (S1), (3.1), but not (S2). If the pressures $`P_L`$ and $`P_R`$ in (2.10) differ by a small amount $`\mathrm{\Delta }=|P_LP_R|`$, then one can estimate how far the piston can swing. The piston can move as long as

$$\frac{X_0^3}{X(t)^2}_0^{\mathrm{}}v^2p(X(t)0,v,t)𝑑v+\frac{(1X_0)^3}{(1X(t))^2}_0^{\mathrm{}}v^2p(X(t)+0,v,t)𝑑v2E(0)$$

where $`X_0=X(0)`$, as before, and $`E(0)`$ is the initial total energy:

$$2E(0)=X_0_0^{\mathrm{}}v^2p_L(v)𝑑v+(1X_0)_0^{\mathrm{}}v^2p_R(v)𝑑v$$

By simple calculations one obtains, to the leading order of $`\mathrm{\Delta }`$, the following bound on the piston displacements:

$$|X(t)X_0|\frac{2\mathrm{\Delta }}{3P_L}\left(\frac{1}{X_0}+\frac{1}{1X_0}\right)^1+O(\mathrm{\Delta }^2)$$

## 4 Perturbative analysis

Here we consider initial densities $`p(x,v,0)`$ satisfying (S2) and the following stricter version of (S1):

* Full uniformity and symmetry. The initial density $`p(x,v,0)`$ is uniform in $`x`$ across the entire cylinder, i.e. $`p(x,v,0)=p_0(|v|)`$ for all $`v`$ and $`0<x<1`$.

We also assume that the piston is initially at the midpoint $`X(0)=0.5`$. Of course, under the conditions (S2)–(S3), the hydrodynamical equations (H1)–(H4) have a simple stationary solution (2.11). On the other hand, we no longer assume monotonicity (3.1). We use perturbative analysis to investigate the stability of the stationary solution (2.11).

From now on we denote by $`p_0(v)=p_0(|v|)`$ an initial density satisfying (S2)–(S3) and by $`p(x,v,0)`$ a perturbed initial density, which we write down as

$$p(x,v,0)=p_0(|v|)+\epsilon p_1(x,v,0)$$

where $`\epsilon p_1(x,v,0)`$ is a small perturbation. We will work to first order in $`\epsilon `$, i.e. ignore terms of order $`o(\epsilon )`$. For $`t>0`$, we decompose the density $`p(x,v,t)`$ as

$$p(x,v,t)=p_0(|v|)+\epsilon p_1(x,v,t)$$

We also set $`p_1(x,v,t)=p_L(x,v,t)`$ for $`x<X(t)`$ and $`p_1(x,v,t)=p_R(x,v,t)`$ for $`x>X(t)`$.

According to (2.8), the velocity $`V(t)`$ of the piston is given by

$$_V^{\mathrm{}}(vV)^2[p_0(v)+\epsilon p_L(X,v,t)]𝑑v=_{\mathrm{}}^V(vV)^2[p_0(v)+\epsilon p_R(X,v,t)]𝑑v$$

where $`X=X(t)`$ is the position of the piston. Expanding in $`\epsilon `$ and ignoring terms of order $`o(\epsilon )`$ gives

$$V(t)=\epsilon \frac{_0^{\mathrm{}}v^2p_L(X,v,t)𝑑v_{\mathrm{}}^0v^2p_R(X,v,t)𝑑v}{4_0^{\mathrm{}}vp_0(v)𝑑v}$$

Note that by integrating by parts we obtain (for piecewise smooth $`p_0`$)

$$2_0^{\mathrm{}}vp_0(v)𝑑v=_0^{\mathrm{}}v^2p_0^{}(v)𝑑v$$

We define

$$h(v)=p_0^{}(v)\mathrm{for}v>0$$

and for the sake of completeness, set $`h(v)=h(v)`$. Then

$$V(t)=\epsilon \frac{_0^{\mathrm{}}v^2p_L(X,v,t)𝑑v_{\mathrm{}}^0v^2p_R(X,v,t)𝑑v}{2_0^{\mathrm{}}v^2h(v)𝑑v}$$
(4.1)
The function $`p_0(v)`$ does not have to be differentiable, we interpret $`h(v)`$ here as the generalized derivative of $`p_0(v)`$. Denote by $`_+`$ the integration $`_0^{\mathrm{}}dv`$, and by $`_{}`$ the integration $`_{\mathrm{}}^0dv`$. Then

$$V(t)=\epsilon \frac{v^2p_L(X,v,t)_+v^2p_R(X,v,t)_{}}{2v^2h(v)_+}$$

The density of the gas after interaction with the piston is given by the formulas (2.5), which imply:

$`p(X0,v,t)`$ $`=`$ $`p(X0,v+2V,t)`$
$`=`$ $`p_0(v+2V)+\epsilon p_L(X,v,t)`$
$`=`$ $`p_0(v)+2Vp_0^{}(v)+\epsilon p_L(X,v,t)`$
$`=`$ $`p_0(v)2Vh(v)+\epsilon p_L(X,v,t)`$
Here we assume $`v>0`$ and ignore terms of order $`o(\epsilon )`$. Hence, the “reflection rule” can be written as

$$p_L(X,v,t)=p_L(X,v,t)h(v)\frac{v^2p_L(X,v,t)_+v^2p_R(X,v,t)_{}}{v^2h(v)_+}$$

This expression suggests the introduction of new functions:

$$q_{R,L}(x,v,t)=\frac{p_{L,R}(x,v,t)}{h(v)}$$

and

$$\rho (v)=\frac{v^2h(v)}{v^2h(v)_+}$$

The above expression for $`p_L(X,v,t)`$ can now be written as

$$q_L(X,v,t)=q_L(X,v,t)q_L(X,v,t)\rho (v)_++q_R(X,v,t)\rho (v)_{}$$
(4.2)
Similarly, on the other side of the piston,

$$q_R(X,v,t)=q_R(X,v,t)+q_L(X,v,t)\rho (v)_+q_R(X,v,t)\rho (v)_{}$$
(4.3)
One can interpret these “reflection rules” as follows: the functions $`q_L`$ and $`q_R`$ “exchange” their average values with respect to the “density” $`\rho (v)`$.

Note that $`\rho (v)`$ is normalized, so that $`\rho (v)_+=1`$, but it is not necessarily positive (or even nonnegative). On the other hand when $`\rho (v)0`$, i.e. the unperturbed density $`p_0(|v|)`$ is nonincreasing, thus satisfying (3.1). In this case the stationary solution (2.11) is stable, as we already know by Theorem 3.1. Here we recover this result by our perturbative analysis:

###### Theorem 4.1

The quantity

$$Q=\frac{q^2(x,v,t)\rho (v)}{|v|}𝑑x𝑑v$$

is constant in time, i.e. $`dQ/dt=0`$. Here $`q=q_L`$ for $`x<X`$ and $`q=q_R`$ for $`x>X`$.

Proof. Clearly, $`Q`$ cannot change just due to the free motion of the gas or due to collisions with the walls, so we only need to worry about collisions with the piston. The gas particles colliding with the piston during an infinitesimal interval $`(t,t+dt)`$ lie in the two triangles on the $`xv`$ plane: $`Xvdt<x<X`$ for $`v>0`$ and $`X<x<Xvdt`$ for $`v<0`$. The outgoing particles lie in similar symmetric triangles. Hence, during the interval $`(t,t+dt)`$, the quantity $`Q`$ decreases by (up to the factor of $`dt`$)

$$_0^{\mathrm{}}|v|\frac{q_L^2(X,v,t)\rho (v)}{|v|}𝑑v+_{\mathrm{}}^0|v|\frac{q_R^2(X,v,t)\rho (v)}{|v|}𝑑v$$

$$=q_L^2(X,v,t)\rho (v)_++q_R^2(X,v,t)\rho (v)_{}$$

and it increases by

$$_{\mathrm{}}^0|v|\frac{q_L^2(X,v,t)\rho (v)}{|v|}𝑑v+_0^{\mathrm{}}|v|\frac{q_R^2(X,v,t)\rho (v)}{|v|}𝑑v$$

$$=q_L^2(X,v,t)\rho (v)_{}+q_R^2(X,v,t)\rho (v)_+$$

After substituting (4.2) and (4.3) into the above expressions for $`q_L`$ and $`q_R`$ and some manipulations, all changes in $`Q`$ cancel out and so it stays constant. $`\mathrm{}`$

When $`p_0(|v|)`$ is strictly decreasing, hence $`\rho (v)>0`$, then $`Q`$ is a norm in the space of functions. Thus, the above theorem implies linear stability.

When $`p_0(|v|)`$ is decreasing, but not strictly, then $`\rho (v)0`$, but there may be regions where $`\rho (v)=0`$. They correspond to the intervals where $`p_0^{}=0`$, i.e. where $`p_0`$ is constant. On such intervals, the reflection rules (4.2)–(4.3) for the perturbed density $`p_{L,R}`$ are trivial:

$$p_L(X,v,t)=p_L(X,v,t)\mathrm{and}p_R(X,v,t)=p_R(X,v,t)$$

In this case $`p_L`$ and $`p_R`$ cannot grow either. Therefore, we obtain linear stability for all nonincreasing $`p_0(|v|)`$.

Next we turn to unstable solutions. The stationary solution for an initial density $`p_0(|v|)`$ satisfying (S2)–(S3) is linearly unstable if some small perturbations grow exponentially in time, i.e. $`p_1(x,v,t)\mathrm{\Lambda }^t`$ for some $`p_1(x,v,0)`$ and $`\mathrm{\Lambda }>1`$. This is equivalent to having a positive Lyapunov exponent in the subspace spanned by the function $`p_1`$ and its images. To investigate the existence of such perturbations we first simplify the collision rules (4.2) and (4.3). Consider the following “symmetric” and “antisymmetric” linear combinations of $`q_L`$ and $`q_R`$:

$$q_+(x,v,t)=[q_L(x,v,t)+q_R(1x,v,t)]/2$$

and

$$q_{}(x,v,t)=[q_L(x,v,t)q_R(1x,v,t)]/2$$

They are defined for $`x<1/2`$. The collision rules (4.2)–(4.3) now take form

$$q_+(X,v,t)=q_+(X,v,t)$$
(4.4)
and

$$q_{}(X,v,t)=q_{}(X,v,t)2q_{}(X,v,t)\rho (v)_+$$
(4.5)
Hence $`q_+`$ is simply a periodic function in $`t`$, so it cannot grow to infinity or decrease to zero. In other words, it cannot affect the stability or instability of the hydrodynamical equations. The latter is determined by $`q_{}`$ alone. So we will only consider $`q_{}`$ and omit “$``$” for brevity. Our collision rule then reduces to a single equation

$$q(X,v,t)=q(X,v,t)2q(X,v,t)\rho (v)_+$$
(4.6)
Next we demonstrate, by example, that densities $`\rho _0(|v|)`$ for which the stationary solution (2.11) is unstable do exist.

Example. Let $`p_0`$ be a rectangular function defined by

$`p_0(v)=\{\begin{array}{cc}1\hfill & \mathrm{if}0.5<|v|<1\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}`$ (4.9)
This $`p_0(v)`$ satisfies (S2)–(S3) but not (3.1). We will show that the corresponding stationary solution is linearly unstable.

First, the function $`h=p_0^{}`$ is the sum of two delta functions:

$$h(v)=\delta _{0.5}+\delta _1$$

(and symmetrically for $`v<0`$). It is easy to compute $`\rho `$ directly

$$\rho =\frac{1}{3}\delta _{0.5}+\frac{4}{3}\delta _1$$

Now the reflection rule (4.6) implies

$`q(1)`$ $`=`$ $`{\displaystyle \frac{5}{3}}q(1)+{\displaystyle \frac{2}{3}}q(0.5)`$
$`q(0.5)`$ $`=`$ $`{\displaystyle \frac{8}{3}}q(1)+{\displaystyle \frac{5}{3}}q(0.5)`$
Note that only the values $`p(x,\pm 0.5,t)`$ and $`p(x,\pm 1,t)`$ will evolve in a nontrivial way, as specified above, since $`h(v)=0`$ for all $`v\{1,0.5,0.5,1\}`$.

We now construct a linear subspace of functions $`q=q_{}`$ that stays invariant under the above transformations and in which functions grow exponentially in time (since the $`q_+`$ component of the perturbation is irrelevant, we set it to zero). We can simplify the construction further by assuming that at time $`t=0`$

$$q(x,\pm 1,0)=u_1,q(x,0.5,0)=u_2,q(x,0.5,0)=u_3$$

with some constants $`u_1,u_2,u_3`$ (the choice of indices $`1,2,3`$ is rather arbitrary). We note that the functions $`p_{L,R}(x,v,0)`$ are now piecewise constant and are completely described by the values $`u_1,u_2,u_3`$. The space of such perturbations is three-dimensional.

It is easy to see that at time $`t=1`$ the functions $`p_{L,R}`$ will again be constant on the same intervals, hence they will be described by some other constants $`u_1^{},u_2^{},u_3^{}`$. Our collision rule (4.6) implies that the vectors $`𝐮^{}=(u_1^{},u_2^{},u_3^{})^T`$ and $`𝐮=(u_1,u_2,u_3)^T`$ are related by a linear transformation

$$𝐮^{}=\mathrm{𝐀𝐮}$$

where $`𝐀`$ is a $`3\times 3`$ matrix:

$$𝐀=\frac{1}{3}\left(\begin{array}{ccc}\hfill 5& \hfill 2& \hfill 0\\ \hfill 0& \hfill 0& \hfill 3\\ \hfill 8& \hfill 5& \hfill 0\end{array}\right)$$

After that, the evolution will proceed periodically – the vector $`𝐮`$ will be multiplied by the matrix $`𝐀`$ at times $`t=1,2,3,\mathrm{}`$. The matrix $`𝐀`$ has three real eigenvalues:

$$\lambda _{1,2}=\frac{4\pm \sqrt{7}}{3}\mathrm{and}\lambda _3=1$$

The largest eigenvalue $`\lambda =(4+\sqrt{7})/32.215`$ has the following (unit) eigenvector:

$$𝐮=(0.4472,0.3680,\mathrm{\hspace{0.17em}0.8152})$$

This eigenvector spans a one-dimensional subspace in the space of perturbation densities, which is invariant over the time period $`t=1`$ and in which the corresponding perturbations are expanded by a factor $`|\lambda |2.215`$. Roughly, the perturbations double over one period.

To explore the above periodic growth of perturbations, we note that the piston velocity is given by

$`V(t)`$ $`=`$ $`{\displaystyle \frac{\epsilon }{2}}\left(q_L(X,v,t)\rho (v)_+q_R(X,v,t)\rho (v)_{}\right)`$
$`=`$ $`\epsilon q(X,v,t)_+`$
Hence in our example, during the time interval $`0<t<1`$

$$V=\frac{\epsilon }{3}\left(4u_1u_3\right)=0.9736\epsilon $$

During the next time interval $`1<t<2`$ we have

$$V=\frac{\epsilon }{3}\left(4u_1^{}u_3^{}\right)=2.156\epsilon $$

and so on. Hence, over a unit period of time, the piston velocity grows by a factor of $`|\lambda |=2.215`$ and changes sign – the piston starts its movements back and forth (oscillations) that increase exponentially in time.

We note that the same density (4.9) was studied in \[CL\] where the trajectory of the piston was computed after an initial configuration of gas molecules was selected randomly from the distribution $`p_0(v)`$ given in (4.9). It was found \[CL\] that the piston indeed made oscillations which increased exponentially in time. The piston’s velocity grew as const$`R^t`$ with $`R1.6`$. This experimental estimate is in a reasonable agreement with our calculation of the largest eigenvalue $`2.215`$.

Next, we modify the unstable perturbations $`q`$ found above and make them smooth (rather than piecewise constant) functions of $`v`$.

We will be looking for the function $`q`$ of the form

$$q(x,v,t)=C(v)e^{z(tx/v)}$$

where $`z`$ is a complex constant. Note that due to (2.2) the function $`q`$ (with $`v`$ fixed) can only depend on $`tx/v`$. We chose the exponential form in order to investigate the existence of solutions of the linear equation which grow exponentially with time. Also, for convenience, we introduce the new space coordinate $`y`$ in the following way: for all $`v>0`$ and $`x<0.5`$ we set $`y=x+0.5`$, for $`v<0`$ and $`x<0.5`$ we set $`y=0.5x`$, for $`v>0`$ and $`x>0.5`$ we set $`y=x0.5`$, and for $`v<0`$ and $`x>0.5`$ we set $`y=1.5x`$. The so designed coordinate $`y`$ assumes value zero when a point $`(x,v)G`$ moving under $`F^t`$ reflects off the piston, then grows to 0.5 when the point travels to the wall, and grows further from 0.5 to 1 when the point travels from the wall back to the piston.

In the new coordinate $`y`$, we will be looking for the function $`q`$ of the form

$$q(y,v,t)=C(|v|)e^{z(ty/|v|)}$$

More precisely, let

$$q(y,\pm 1,t)=C(1)e^{z(ty)}$$

$$q(y,\pm 0.5,t)=C(0.5)e^{z(t2y)}$$

Recall that $`p_0(|v|)`$ is the characteristic function of the interval $`[0.5,1]`$.

Now, the reflection rule (4.6) implies

$$C(1)=\frac{5}{3}C(1)e^z+\frac{2}{3}C(0.5)e^{2z}$$

$$C(0.5)=\frac{8}{3}C(1)e^z+\frac{5}{3}C(0.5)e^{2z}$$

We need to find $`z`$ for which the above system of equations has a nontrivial solution. Put $`e^z=\lambda `$ and introduce an auxiliary variable $`D(0.5)=C(0.5)e^z`$. Now the above system can be rewritten as

$`\lambda C(1)`$ $`=`$ $`{\displaystyle \frac{5}{3}}C(1)+{\displaystyle \frac{2}{3}}D(0.5)`$
$`\lambda D(0.5)`$ $`=`$ $`C(0.5)`$
$`\lambda C(0.5)`$ $`=`$ $`{\displaystyle \frac{8}{3}}C(1)+{\displaystyle \frac{5}{3}}D(0.5)`$
Hence, $`\lambda `$ is an eigenvalue of the matrix of coefficients

$$\frac{1}{3}\left(\begin{array}{ccc}\hfill 5& \hfill 2& \hfill 0\\ \hfill 0& \hfill 0& \hfill 3\\ \hfill 8& \hfill 5& \hfill 0\end{array}\right)$$

which is the same matrix $`𝐀`$ that we encountered before. We take its leading eigenvalue $`|\lambda |>1`$ and set

$$z=\mathrm{ln}|\lambda |+i\pi $$

The function $`q`$ now takes form

$$q(y,v,t)=\pm C(|v|)|\lambda |^{ty/|v|}\mathrm{cos}\pi (ty/|v|)$$

where $`C(0.5)`$ and $`C(1)`$ are the coordinates of the leading eigenvector, and we only take the real part, for obvious reasons. Since $`|\lambda |>1`$, we have an exponential growth of perturbations and thus linear instability. This gives us smooth unstable perturbations.

We now generalize the above construction to arbitrary nonmonotonic initial densities $`p_0`$. Let $`p_0(v)`$ satisfy (S2)–(S3) but not (3.1). We will be looking for perturbations of the form

$$q(y,v,t)=C(|v|)e^{z(ty/|v|)}$$
(4.10)
with the same convention on $`y`$ as before. The reflection rule (4.6) leads to (cancelling $`e^{zt}`$)

$$C(v)=C(v)e^{z/v}2_0^{\mathrm{}}C(v)e^{z/v}\rho (v)𝑑v$$

for all $`v>0`$. Denoting

$$D=2C(v)e^{z/v}\rho (v)𝑑v$$

gives immediately

$$C(v)=\frac{D}{1e^{z/v}}$$

Thus, we not only eliminated $`t`$ but determined the function $`C(v)`$ up to a constant factor. The above solution exists if

$$D=2_0^{\mathrm{}}\frac{De^{z/v}\rho (v)}{1e^{z/v}}𝑑v$$

or, cancelling $`D`$,

$$_0^{\mathrm{}}\frac{\rho (v)}{1e^{z/v}}𝑑v=\frac{1}{2}$$
(4.11)
If this equation has a solution $`z`$ with Re$`(z)>0`$, we immediately obtain an unstable perturbation (4.10). Otherwise our construction of unstable perturbations does not work.

Unfortunately, it does not seem to be easy to solve equation (4.11) for particular functions $`\rho (v)`$ or even to determine if it has solutions with a positive real part, as we will demonstrate in the next section.

Next we mention an important property of (4.11). Let us denote

$$F(z):=_0^{\mathrm{}}\frac{\rho (v)}{1e^{z/v}}𝑑v\frac{1}{2}$$
(4.12)
###### Lemma 4.2

$`F(z)+F(z)=0`$ for all $`z\text{C}\mathit{|}\text{ }\mathit{|}\text{ }`$.

Proof.

$`F(z)+F(z)`$ $`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{\rho (v)}{1e^{z/v}}}𝑑v{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{e^{z/v}\rho (v)}{1e^{z/v}}}𝑑v1`$
$`=`$ $`{\displaystyle _0^{\mathrm{}}}{\displaystyle \frac{(1e^{z/v})\rho (v)}{1e^{z/v}}}𝑑v1`$
$`=`$ $`0\mathrm{}`$
As a result, the existence of a solution of (4.11) with Re$`(z)>0`$ is equivalent to that of a solution with Re$`(z)<0`$. The alternative is when all the solutions lie on the imaginary axis Re$`(z)=0`$.

## 5 A special family of densities

Here we investigate a family of rectangular densities

$`p_0(v)=\{\begin{array}{cc}1\hfill & \mathrm{if}r<|v|<1\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array}`$ (5.3)
where $`0<r<1`$ is the parameter of our family. Note that our example (4.9) is a particular case of (5.3) with $`r=1/2`$. It is easy to compute

$$h(v)=p_0^{}(v)=\delta _r(v)\delta _1(v)$$

and

$$\rho (v)=\frac{v^2h(v)}{_0^{\mathrm{}}v^2h(v)𝑑v}=\frac{1}{1r^2}\left[\delta _1(v)r^2\delta _r(v)\right]$$

Since $`h(v)=0`$ for all $`v\{1,r,1,r\}`$, we only consider perturbations $`q(x,v,t)`$ defined for $`v=1,r,1,r`$. The reflection rule (4.6) now gives

$`q(X,1,t)`$ $`=`$ $`\alpha q(X,1,t)+\beta q(X,r,t)`$
$`q(X,r,t)`$ $`=`$ $`\gamma q(X,1,t)+\alpha q(X,r,t)`$
where

$$\alpha =\frac{1+r^2}{1r^2},\beta =\frac{2r^2}{1r^2},\gamma =\frac{2}{1r^2}$$
(5.4)
While we cannot construct unstable perturbations for arbitrary irrational values of $`r`$, it is relatively easy to investigate the case of rational $`r`$. But even in this case, we can only provide partial answers leaving some questions open.

Let $`r=m/n`$ be a rational number with $`1m<n`$. Now (5.4) takes form

$$\alpha =\frac{n^2+m^2}{n^2m^2},\beta =\frac{2m^2}{n^2m^2},\gamma =\frac{2n^2}{n^2m^2}$$
(5.5)
To investigate the evolution of perturbations $`q(x,v,t)`$ as $`t`$ grows, we consider $`n+m`$ points $`P_iG`$, $`1in+m`$, shown on Fig. 1. The points $`P_i`$ are defined as follows:

$$P_i=\{\begin{array}{cc}(1,0.5(i1)/m)\hfill & \mathrm{for}1i<m/2+1\hfill \\ (1,(i1)/m0.5)\hfill & \mathrm{for}m/2+1im\hfill \\ (r,0.5(im1)/n)\hfill & \mathrm{for}m<i<m+n/2+1\hfill \\ (r,(im1)/n0.5)\hfill & \mathrm{for}m+n/2+1im+n\hfill \end{array}$$

It is crucial to observe that the points $`P_i`$ move under the dynamics in a periodic fashion. In a time period $`\mathrm{\Delta }t=1/m`$, the point $`P_i`$ is mapped to $`P_{i+1}`$ for all $`1i<m`$ and all $`m+1i<m+n`$. Also, $`P_m`$ moves to the piston, gets reflected off it and lands on $`P_1`$. Likewise, $`P_{m+n}`$ moves to the piston, gets reflected off it and lands on $`P_{m+1}`$. Therefore, the time shift $`\mathrm{\Delta }t`$ permutes the points $`P_i`$, $`1im+n`$, in two independent cycles. The reason why we combine the two cycles together is that they are linked by the reflection rule, as we will see shortly.

For each $`i`$, denote $`q_i(t)=q(P_i,t)`$. Then we have

$$q_i(t+\mathrm{\Delta }t)=q_{i1}(t)$$

for all $`2im`$ and $`m+2im+n`$. The reflection rule now implies

$`q_1(t+\mathrm{\Delta }t)`$ $`=`$ $`\alpha q_m(t)+\beta q_{m+n}(t)`$
$`q_{m+1}(t+\mathrm{\Delta }t)`$ $`=`$ $`\gamma q_m(t)+\alpha q_{m+n}(t)`$
Thus, the vector $`𝐪(t)=(q_1(t),\mathrm{},q_{n+m}(t))`$ is updated at time $`t+\mathrm{\Delta }t`$ by the rule

$$𝐪(t+\mathrm{\Delta }t)=𝐁𝐪(t)$$

where $`𝐁`$ is an $`(n+m)\times (n+m)`$ matrix,

$$𝐁=\left(\begin{array}{cccccccc}0& \mathrm{}& & \alpha & & \mathrm{}& & \beta \\ 1& \mathrm{}& & & & & & \\ \mathrm{}& \mathrm{}& 0& \mathrm{}& \mathrm{}& & & \mathrm{}\\ & & 1& 0& & \mathrm{}& & \\ 0& \mathrm{}& & \gamma & 0& \mathrm{}& & \alpha \\ & & & & 1& 0& & \\ \mathrm{}& & & \mathrm{}& & \mathrm{}& \mathrm{}& \mathrm{}\\ 0& \mathrm{}& & 0& & \mathrm{}& 1& 0\end{array}\right)$$

We conclude that the existence of unstable perturbations $`q(t)`$ is equivalent to the existence of an eigenvalue $`\lambda `$ of $`𝐁`$ such that $`|\lambda |>1`$. The characteristic polynomial of the matrix $`𝐁`$ is

$$P(\lambda )=\lambda ^{m+n}+\alpha \lambda ^n\alpha \lambda ^m1$$
(5.6)
where $`\alpha =(n^2+m^2)/(n^2m^2)`$ is defined in (5.5).

Remark. Interestingly, the equation (4.11) can be reduced to $`P(\lambda )=0`$ as well. Indeed, it is easy to see that

$$_0^{\mathrm{}}\frac{\rho (v)}{1e^{z/v}}𝑑v=\frac{1}{1r^2}\left[\frac{1}{1e^z}\frac{r^2}{1e^{z/r}}\right]$$

Now the substitution $`\lambda =e^{z/m}`$ and some algebraic manipulations show that Eq. (4.11) is equivalent to $`P(\lambda )=0`$.

It is easy to see that if $`\lambda `$ is a root of $`P(\lambda )`$, then so is $`1/\lambda `$ (this reciprocability also follows from Lemma 4.2). Thus, the existence of unstable perturbations is equivalent to the existence of eigenvalues of $`𝐁`$ that do not lie on the unit circle $`|\lambda |=1`$.

If an eigenvalue $`|\lambda |>1`$ of $`𝐁`$ exists, then the perturbations in the corresponding eigenspace grow by the factor of $`|\lambda |`$ over the time period $`\mathrm{\Delta }t=1/m`$. Hence, the expansion factor per unit time would be $`\mathrm{\Lambda }=|\lambda |^m`$.

###### Theorem 5.1

Let $`r=m/n`$ be a rational number with an even denominator $`n`$ (hence, $`m`$ is odd). Then there is a unique eigenvalue of $`𝐁`$ such that $`\lambda <1`$. This eigenvalue has multiplicity one. The expansion factor per unit time $`\mathrm{\Lambda }_r=|\lambda |^m`$ depends on $`r`$ continuously, and we have, asymptotically,

$$\mathrm{\Lambda }_r=1+\mathrm{const}r^{3/2}+𝒪(r^2)\mathrm{as}r0$$
(5.7)
and

$$\mathrm{\Lambda }_r\frac{\mathrm{const}}{1r}\mathrm{as}r1$$
(5.8)
Proof. One can easily check that, under the conditions of the theorem, $`P(1)>0`$ and $`P(\mathrm{})<0`$, hence a root $`\lambda <1`$ exists. Next,

$$P^{}(\lambda )=[(n+m)\lambda ^n+\alpha n\lambda ^{nm}\alpha m]\lambda ^{m1}$$

and so $`P^{}(1)<0`$ and $`P^{}(\mathrm{})>0`$. Now let $`Q(\lambda )=(n+m)\lambda ^n+\alpha n\lambda ^{nm}\alpha m`$, then

$$Q^{}(\lambda )=[(n+m)\lambda ^m+\alpha (nm)]n\lambda ^{nm1}$$

and so clearly $`Q^{}(\lambda )<0`$ for all $`\lambda <1`$. Putting these facts together proves the uniqueness and the simplicity of the root $`\lambda <0`$.

The equation $`P(\lambda )=0`$ can be rewritten in terms of $`\mathrm{\Lambda }_r=|\lambda |^m`$ as follows:

$$\mathrm{\Lambda }_r^{1+1/r}+\frac{1+r^2}{1r^2}\mathrm{\Lambda }_r^{1/r}+\frac{1+r^2}{1r^2}\mathrm{\Lambda }_r1=0$$
(5.9)
Now the continuity of $`\mathrm{\Lambda }_r`$, as a function of $`r`$, is obvious. Note that our argument is only valid when $`r=m/n`$ with an even $`n`$ and an odd $`m`$, because this parity condition dictates the signs in (5.9).

To prove (5.7), one can substitute $`\mathrm{\Lambda }_r=1+\epsilon `$ into (5.9) and expand all the terms in Taylor series, the calculation is then straightforward and we omit it. The proof of (5.8) is similar. $`\mathrm{}`$

Figure 2 presents the graph of the Lyapunov exponent $`\mathrm{log}\mathrm{\Lambda }_r`$ as a function of $`r`$.

###### Lemma 5.2

Let $`z`$ be a solution of (4.11) such that $`|e^z|1`$ and $`e^z\text{I}\text{R}`$. Then $`dF/dz0`$ (in fact, $`dF/dz`$ is a real negative number).

Proof. A direct calculation shows that

$$\frac{dF}{dz}=\frac{(1+r+r^2)(1+e^z)^2+r^3(1e^z)^2}{4r^2(1+r)(1e^z)^2}$$

which proves the lemma. $`\mathrm{}`$

For any $`r=m/n`$ with an even $`n`$ and odd $`m`$, the corresponding solution $`e^z=\mathrm{\Lambda }_r`$ satisfies the conditions of the above lemma. Hence, this solution changes continuously with $`r`$, and so we get the following

###### Corollary 5.3

For every $`r=m/n`$ with an even $`n`$ and an odd $`m`$ there is an interval $`(r\epsilon ,r+\epsilon )`$ in which all parameter values have unstable perturbations.

Therefore, unstable perturbations exist for an open and dense set of parameter values $`0<r<1`$. One would naturally wonder if all $`r`$’s have unstable perturbations. The answer is, surprisingly, negative:

Fact 5.4 For the density (5.3) with $`r=1/3`$, there are no solutions $`z`$ of (4.9) with $`Re(z)>0`$, hence there are not solutions of the linearized equation which grow exponentially with time.

Proof. The characteristic equation

$$\lambda ^4+\frac{5}{4}\lambda ^3\frac{5}{4}\lambda 1=0$$

has two real roots ($`\lambda =\pm 1`$) and two complex roots. The complex roots are, on the one hand, conjugate and, on the other, satisfy the reciprocability rule: $`P(\lambda )=0`$ if and only if $`P(1/\lambda )=0`$. Hence, they must belong to the unit circle $`|\lambda |=1`$. $`\mathrm{}`$

It is interesting to know if other rational parameter values $`r=m/n`$ with odd $`n`$ are also stable. We have examined the values $`r=1/n`$ for small odd values of $`n=5,7,\mathrm{},31`$ numerically (by using MATLAB) and always found that all the roots of $`P(\lambda )`$ belonged to the unit circle. Therefore, we conjecture that the values $`r=1/n`$ with odd $`n`$ are stable.

On the other hand, the values $`r=m/n`$ with an odd $`n`$ but $`m>1`$ appear to be unstable. For $`r=2/3`$, $`2/5`$, $`3/5`$, $`3/7`$ we found, again numerically (by using MATLAB), roots $`\lambda `$ such that $`|\lambda |>1`$. All those roots are complex, for example for $`r=2/3`$ they are $`\lambda =0.3778\pm 1.7173i`$. It remains to determine theoretically whether all rational values $`r=m/n`$ with $`m>1`$ are unstable, we leave this question open.

Fact 5.4 seems to disagree with Theorem 5.8. Indeed, let $`p_0(v)`$ be the rectangular density (5.3) corresponding to $`r=1/3`$ and $`p(x,v,0)=p_0(v)+\epsilon p_1(x,v,0)`$ an arbitrary perturbation with an infinitesimally small $`\epsilon `$. According to Fact 5.4, this perturbation cannot grow exponentially in time. On the other hand, let us approximate $`1/3`$ by a rational number $`r=m/n`$ with an even $`n`$. Denote by $`p_0^{}(v)`$ the corresponding rectangular density (5.3) for the chosen $`r=m/n`$. Then we have

$$p(x,v,0)=p_0^{}(v)+\epsilon p_2(x,v,0)\mathrm{with}p_2=p_1+(p_0p_0^{})/\epsilon $$
(5.10)
Hence, if $`|r1/3|<\epsilon `$, then $`(p_0p_0^{})/\epsilon `$ is of order one (in the $`L^1`$ metric), and $`p(x,v,0)`$ becomes an $`\epsilon `$-perturbation of the density $`p_0^{}(v)`$. As such, it “must” grow exponentially in time according to Theorem 5.8. This apparent disagreement requires an explanation, which we provide next.

We recall that smooth unstable perturbations are given by the general formula (4.10). For the rectangular density (5.3), the velocity $`v`$ in this formula only takes two values, $`|v|=r`$ and $`|v|=1`$, hence the factor $`C(|v|)`$ takes two values, as well, and so plays little role. For simplicity, we set $`|v|=1`$ and ignore the constant factor $`C(|v|)=C(1)`$. Now the (real part of) unstable perturbations is described by

$$q(y,1,t)=\mathrm{Re}e^{z(ty)}=e^{(\mathrm{Re}z)(ty)}\mathrm{cos}[(\mathrm{Im}z)(ty)]$$
(5.11)
A similar formula holds for $`|v|=r`$, and we omit it. Now recall that for any rational $`r=m/n`$ we have $`e^{z/m}=\lambda `$, where $`\lambda <1`$ is the eigenvalue of $`𝐁`$ described by Theorem 5.8. Therefore, Re$`z=m\mathrm{log}|\lambda |=\mathrm{log}\mathrm{\Lambda }`$ and Im$`z=\pm m\pi `$.

We see that the real part of $`z`$ changes continuously with $`r=m/n`$ but the imaginary part does not. In particular, when $`r=m/n`$ is close to $`1/3`$ and $`n`$ is even, both $`m`$ and $`n`$ have to be large, so that $`m\mathrm{}`$ and $`|\mathrm{Im}z|\mathrm{}`$ as $`r1/3`$. In terms of the perturbation (5.11), the growth of $`|\mathrm{Im}z|`$, as $`r`$ approaches $`1/3`$, implies that the function $`q(y,1,t)`$ becomes highly oscillatory, and so does the corresponding initial unstable perturbation $`p(x,v,0)=h(v)q(x,v,0)`$. Thus, the linear subspace of unstable perturbations (along which exponential growth takes place) becomes nearly orthogonal to any given function, in particular to $`p_2(x,v,0)`$ defined in (5.10).

This explains the above “disagreement”. The density $`p_2`$ does grow exponentially in time for any $`r=m/n`$ with an even $`n`$, but, as $`r1/3`$, the projection of $`p_2`$ onto the unstable subspace corresponding to the positive Lyapunov exponent $`\mathrm{log}\mathrm{\Lambda }_r>0`$ becomes small and vanishes in the limit, hence the exponential growth is not visible during a long initial interval of time. In the limit $`r1/3`$ that “initial interval” becomes infinite and the instability evaporates.

One can also reverse this line of argument. Indeed, when $`\epsilon `$ is not infinitesimally small but finite, the representation (5.10) implies that any perturbation $`p(x,v,0)`$ of the rectangular density (5.3) for any $`0<r<1`$ will eventually grow exponentially fast in time (because any $`r(0,1)`$ can be approximated by rational numbers $`m/n`$ with even $`n`$). We checked this conclusion experimentally and found that it was indeed correct.

Figure 3: Initial rectangular density (5.3) perturbed by a “bump”.

To investigate the instability experimentally, we solved the hydrodynamical equations (H1)–(H4) numerically starting with a perturbed rectangular density (5.3) shown on Fig. 3. The initial density $`p(x,v,0)`$ takes the value one on the black region and zero elsewhere. The small ‘bump” on the top left edge of the upper rectangle represents the perturbation. The area of the bump in our experiments was less than 1% relative to the total area of each black rectangle.

Figure 4: Piston’s trajectory for a perturbed rectangular density with $`r=1/3`$.

Figure 5: Piston’s trajectory for a perturbed rectangular density with $`r=1/4`$.

Figures 4 and 5 show typical trajectories of the piston $`X(t)`$ for $`r=1/3`$ and $`r=1/4`$, respectively. One can see that by the time $`t=15`$ the piston’s motion shows signs of exponential instability in both cases, but there is a notable difference. For $`r=1/4`$ the piston just swings back and forth with a monotonically increasing amplitude, as time goes on. For $`r=1/3`$, the amplitude of oscillations grows slower but the frequency increases quickly. The higher frequency of oscillations of $`X(t)`$ for $`r=1/3`$ probably reflects the oscillatory structure of unstable perturbations for rational $`r=m/n`$ approximating $`1/3`$.

## 6 Periodic solutions of the hydrodynamical equations

Here we discuss the long-term behavior of our system in the unstable regime.

In our previous work \[CL\] we reported the results of computer simulations of the piston dynamics in an ideal gas with many (up to 27 million) particles. The initial configuration of particles was selected randomly with the average density (4.9), see \[CL\] for details. A typical trajectory of the piston $`X(t)`$ found in our experiments is shown here on Fig. 6. One can see that during the initial interval of time $`0<t<8`$ the piston moves back and forth with an exponentially increasing amplitude, which is consistent with our analysis in Section 4, where the density (4.9) was proven to be unstable.

Figure 6: Piston’s trajectory in the mechanical model with $`10^6`$ particles.

Later on, however, at times $`8<t<15`$, the amplitude of the piston’s oscillations decreases to a certain constant value (nearly a half of its maximum attained at $`t=8`$). Then the piston’s oscillations become very stable and continue almost unchanged for a very long time, up to $`t=50`$ or $`100`$, with a very slowly decreasing amplitude.

On the other hand, we have solved the hydrodynamical equations (H1)–(H4) numerically, starting with the same initial density (4.9) perturbed by a bump shown on Fig. 3. Figure 7 presents the resulting trajectory of the piston. One can see that it behaves almost identically to the simulated trajectory of the piston shown on Fig. 6. Thus, not only the initial instability, but also the long term behavior of the simulated piston trajectory match those of perturbed solutions of the hydrodynamical equations.

Figure 7: Piston’s trajectory from the HE for a perturbed rectangular density with $`r=1/2`$.

The behavior shown on Fig. 7 persists when various perturbations of the initial density (4.9) are applied. It seems like there is a periodic cycle or an invariant manifold of quasi-periodic solutions of (H1)–(H4) that acts as an attractor. Of course, due to the time-reversibility of the hydrodynamical equations there can be no attractors in the strict sense. It is more likely that there is an invariant manifold of periodic or quasi-periodic solutions that acts as a saddle point in the phase space: typical trajectories approach that manifold temporarily and then slowly move away. We cannot rigorously prove the existence of periodic or quasi-periodic solutions, but we construct such solutions by using perturbative analysis.

We will be looking for solution of the hydrodynamical equations (H1)–(H4) such that the piston makes harmonic oscillations

$$X(t)=\frac{1}{2}+\epsilon \mathrm{cos}\omega t,\dot{X}(t)=\epsilon \omega \mathrm{sin}\omega t,$$
(6.12)
with some fixed $`\omega >0`$ and small $`\epsilon >0`$. We will approximate such solutions up to the first order in $`\epsilon `$, i.e. ignoring terms of higher order.

The construction is done in two steps. First, we assume that the piston moves as prescribed by (6.12) and consider the motion of a fluid point bouncing against the moving piston $`X(t)`$ and the fixed wall $`x=0`$. Second, we define the density $`p(x,v,t)`$ which, coupled with the piston’s oscillations (6.12), satisfies equations (H1)–(H4).

Let the piston move according to Eqs. (6.12). Then fluid points in the left compartment $`0<x<X(t)`$ bounce against the wall $`x=0`$ and the piston, the latter simply acts on them as a moving wall. It is known that the phase space of gas particles bouncing against a periodically moving wall necessarily contains many invariant curves. Moreover, the region corresponding to high velocities $`|v|>v_0`$ is densely filled by such invariant curves, the larger $`v_0`$ the higher the density of invariant curves. This fact is a consequence of KAM theory, it was first proved by R. Douady in his thesis \[Do\] and later independently by S. Laederich and M. Levi \[LL\]. We describe these invariant curves approximately, up to the first order in $`\epsilon `$, by equation

$$v+\epsilon F(t,V)=V+𝒪(\epsilon ^2)$$
(6.13)
where $`v`$ denotes the velocity of the particle when it kicks the piston, $`t`$ the the collision time, and $`V`$ is the parameter of the curve. In fact, we will construct invariant curves for all $`V>V_0`$ with some $`V_0>0.`$ Here we only consider particles to left to the piston, the particles to the right of the piston are completely symmetric.

Let us consider successive collisions of a gas particle with the piston. Denote by $`v_n>0`$ the velocity of the particle before its $`n`$-th collision and by $`t_n`$ the time of that collision. Then the law of elastic impact reads

$$v_{n+1}=v_n2\dot{X}(t_n)=v_n+2\epsilon \omega \mathrm{sin}\omega t_n$$
(6.14)
Let $`t_{n+1/2}`$ denote the time at which the particle bounces off the wall $`x=0`$ between its $`n`$-th and $`(n+1)`$-st collisions with the piston. Obviously, $`t_{n+1/2}=t_n+X(t_n)/v_{n+1}`$ and $`v_{n+1}(t_{n+1}t_{n+1/2})=X(t_{n+1}).`$ Since we are interested in knowing $`v_n`$ up to terms $`𝒪(\epsilon ),`$ it is sufficient to find $`t_n`$ up to terms $`𝒪(1).`$ This is easy:

$$t_{n+1}=t_n+\frac{1}{v_{n+1}}+𝒪(\epsilon )=t_n+\frac{1}{V}+𝒪(\epsilon )$$
(6.15)
where we used (6.13).

Now let us look for an invariant curve of the form

$$v+\epsilon F(t,V)=V+𝒪(\epsilon ^2)$$

We have to impose the constraint

$$v_{n+1}+\epsilon F(t_{n+1},V)=v_n+\epsilon F(t_n,V)+𝒪(\epsilon ^2)$$
(6.16)
By equations (6.14), (6.15) and (6.16) we get

$$v_n+2\epsilon \omega \mathrm{sin}\omega t_n+\epsilon F(t_n+1/V+𝒪(\epsilon ),V)=v_n+\epsilon F(t_n,V)$$

Cancelling $`v_n`$ and $`\epsilon `$ and removing the index $`n`$ gives a general equation for an invariant curve:

$$2\omega \mathrm{sin}\omega t+F(t+1/V,V)F(t,V)=0$$
(6.17)
We construct solutions of this equation in the form

$$F(t,V)=a\mathrm{cos}\omega t+b\mathrm{sin}\omega t$$
(6.18)
where $`a`$ and $`b`$ depend on $`V`$. By substituting this expression into (6.17) we find that (6.17) can only hold if

$`a(\mathrm{cos}(\omega /V)1)+b\mathrm{sin}(\omega /V)`$ $`=`$ $`0`$
$`a\mathrm{sin}(\omega /V)b(\mathrm{cos}(\omega /V)1)`$ $`=`$ $`2\omega `$
The solution of the above system is

$`a`$ $`=`$ $`{\displaystyle \frac{\omega \mathrm{sin}(\omega /V)}{1\mathrm{cos}(\omega /V)}}`$
$`b`$ $`=`$ $`\omega `$
Remark. Notice that $`a`$ (and hence the invariant curve) is not defined for $`V=\frac{\omega }{2\pi k}`$, $`k=\pm 1,\pm 2,\mathrm{}`$ To avoid these singularities, we will not use invariant curves corresponding to $`V\frac{\omega }{2\pi }`$. In particular, the density $`p(x,v,t)`$ that we define below will be constant for $`|v|\frac{\omega }{2\pi }`$.

Thus, for any $`V>V_0>\frac{\omega }{2\pi }`$ we can define an invariant curve $`u(x,t;V)`$ in the phase space of gas particles, where $`V`$ is the parameter of the curve and $`u(x,t;V)`$ is the velocity of the particle on the curve at point $`x`$ at time $`t`$. The curve is made by two branches: the upper branch $`u^+`$ and the lower branch $`u^{}.`$ Obviously, we have

$`u^+(X(t),t;V)`$ $`=`$ $`V\epsilon F(t,V)+𝒪(\epsilon ^2)`$
$`u^{}(X(t),t;V)`$ $`=`$ $`[V\epsilon F(t,V)2\dot{X}(t)]+𝒪(\epsilon ^2)`$
$`u^+(t,0,V)`$ $`=`$ $`u^{}(t,0,V)`$
Note that the last equation here is equivalent to (6.16).

Now we define a density $`p(x,v,t)`$ so that its value on each invariant curve $`u(x,t;V)`$, $`|V|>V_0`$, is a constant denoted by $`\rho (V)`$. Between the curves $`u^+(x,t;V_0)`$ and $`u^{}(x,t;V_0)`$ we set the density to a constant equal to one. Therefore

$$p(x,u^+(t,x,V),t)=p(x,u^{}(t,x,V),t)=\rho (V)\mathrm{if}V>V_0$$

and

$$p(x,v,t)1\mathrm{if}u^{}(x,t;V_0)<v<u^+(x,t;V_0)$$

The function $`\rho (V)`$ and the “cutoff” value $`V_0>\frac{\omega }{2\pi }`$ will be specified below.

Example. Let us set $`\rho (V)0`$ for $`V>V_0`$, i.e.

$$p(x,v,t)=\{\begin{array}{cc}1\hfill & \mathrm{for}u^{}(x,t;V_0)<u<u^+(x,t;V_0)\hfill \\ 0\hfill & \mathrm{elsewhere}\hfill \end{array}$$

In order to compute the pressure on the piston we only need to know the density $`p(x,v,t)`$ at the point $`x=X(t)`$, i.e. we need to know the function

$$v(t,V):=u^+(X(t),t,V)=V+\epsilon F(t,V)+𝒪(\epsilon ^2)$$

In our example the density on the piston (on the left hand side) is $`1`$ up to $`v^+=V_0\epsilon F(V_0,t)`$. The density on the piston on the right hand side is $`1`$ up to a similar invariant curve, which is phase shifted by $`\mathrm{\Delta }t=\pi /\omega `$. Therefore the density on the right hand side is $`1`$ down to $`v^{}=V_0\epsilon F(t,V_0).`$ Since $`F(t+\pi /\omega )=F(t)`$ by (6.18), the velocity of the piston is exactly the average of $`v^+`$ and $`v^{}`$ and therefore is $`\epsilon F(t,V_0)`$.

Thus our density and the piston satisfy the hydrodynamical equations (H1)–(H4) if $`\dot{X}=\epsilon F(t,V_0)`$, which gives

$$\epsilon \omega \mathrm{sin}\omega t=\epsilon \omega \mathrm{sin}\omega t\frac{\epsilon \omega \mathrm{sin}(\omega /V_0)\mathrm{cos}\omega t}{1\mathrm{cos}(\omega /V_0)}$$

In our example, the only possible choice is $`V_0=\frac{\omega }{\pi }.`$

Now let us consider the case of a generic function $`\rho (V).`$ The pressure on the piston on the left hand side is equal to

$$P_L=_{\dot{X}}^{\mathrm{}}p_L(v)(v\dot{X})^2𝑑v=_0^{\mathrm{}}p_L(v)(v^22v\dot{X})𝑑v+𝒪(\epsilon ^2)$$

where $`p_L(v)=p(X(t)0,v,t)`$ is the density on the piston (we have used the fact that $`\dot{X}=𝒪(\epsilon )`$). Recall that the density is $`p_L(v)=\rho (V)=\rho (v+\epsilon F(t,v))+𝒪(\epsilon ^2).`$ ¿From now on we neglect terms of order $`𝒪(\epsilon ^2)`$. Then we get

$$P_L=_0^{\mathrm{}}(v^22v\dot{X})\rho (v+\epsilon F(t,v))𝑑v$$

The pressure on the right hand side is given, by analogy,

$$P_R=_{\mathrm{}}^0p_R(v)(v^22v\dot{X})𝑑v=_0^{\mathrm{}}p_R(v)(v^2+2v\dot{X})𝑑v$$

Note that for $`v>0`$ we have $`p_R(v)=\rho (V)=\rho (v+\epsilon F(v,t+\pi ))=\rho (v\epsilon F(t,v))`$. Therefore

$$P_R=_0^{\mathrm{}}(v^2+2v\dot{X})\rho (v\epsilon F(t,v))𝑑v$$

We now conclude that $`P_L=P_R`$ iff

$$\dot{X}=\epsilon \frac{_0^{\mathrm{}}\rho ^{}(v)F(t,v)v^2𝑑v}{_0^{\mathrm{}}\rho (v)\mathrm{\hspace{0.17em}2}v𝑑v}=\epsilon \frac{_0^{\mathrm{}}\rho ^{}(v)F(t,v)v^2𝑑v}{_0^{\mathrm{}}\rho ^{}(v)v^2𝑑v}$$

which is analogous to our early formula (4.1).

Using (6.18) and the subsequent equations we find

$$\dot{X}=\epsilon \omega \mathrm{sin}\omega t\epsilon \frac{_0^{\mathrm{}}\rho ^{}(v)\frac{\mathrm{sin}(\omega /v)}{1\mathrm{cos}(\omega /v)}v^2𝑑v}{_0^{\mathrm{}}\rho ^{}(v)v^2𝑑v}\omega \mathrm{cos}\omega t$$

Our density, coupled with the piston oscillations (6.12), satisfies the hydrodynamical equations (H1)–(H4) if and only if $`\dot{X}=\epsilon \omega \mathrm{sin}\omega t`$. This implies

$$_{\omega /2\pi }^{\mathrm{}}𝑑v\rho ^{}(v)v^2\frac{\mathrm{sin}(\omega /v)}{1\mathrm{cos}(\omega /v)}=0$$
(6.19)
where we have imposed $`\rho ^{}=0`$ for $`v\omega /2\pi .`$

Interestingly, (6.19) seems to be related to our early equation (4.11). Precisely, let $`z`$ in (4.11) be a purely imaginary number, $`z=\omega i`$. Also note that $`\rho (v)`$ in (4.11) is just proportional to $`v^2\rho ^{}(v)`$ here. Then (6.19) becomes equivalent to Im$`F(z)=0`$, with $`F(z)`$ defined by (4.12). In other words, (6.19) expresses the “imaginary part” of the equation (4.11). We already observed in the previous section that Im$`z`$ characterized the frequency of oscillations of unstable perturbations, and here $`\omega =`$Im$`z`$ is the frequency of oscillations of the piston. We note that for $`z=\omega i`$ one always has Re$`F(z)=0`$, as it follows from (4.12), hence in our case (6.19) is equivalent to $`F(z)=0`$.

Next, we note that $`\rho ^{}(v)`$ must be supported on both the intervals $`(\omega /2\pi ,\omega /\pi ]`$ and $`[\omega /\pi ,+\mathrm{})`$. In fact the fraction in Eq. (6.19) is negative for $`v(\omega /2\pi ,\omega /\pi )`$ and positive for $`v>\omega /\pi .`$

Obviously Eq. (6.19) can be satisfied in many different ways, all of them leading to different solutions of the hydrodynamical equations.

Lastly, we want to examine a special case when the initial density $`p(x,v,0)=\rho (v)`$ is a monotonic function in $`|v|`$, i.e. when the hydrodynamic equations are stable. Then there are some quantitative restrictions on the period of oscillations of the piston. The period $`T=2\pi /\omega `$ can be bounded from below by a function of the average kinetic energy $`K=K/M,`$ where

$$K=_0^{\mathrm{}}\rho (v)\frac{v^2}{2}𝑑v\mathrm{and}M=_0^{\mathrm{}}\rho (v)𝑑v$$

###### Proposition 6.1

Let $`\rho ^{}0`$ be supported on the interval $`[\omega /2\pi ,\mathrm{})`$ and satisfy Eq. (6.19). Then the period of oscillations $`T`$ is bounded by

$$T\sqrt{\frac{2}{3K}}$$

The equality holds when $`\rho ^{}=(\pi /\omega )\delta (v\omega /\pi )`$, i.e. when $`\rho (|v|)`$ is constant on $`(0,\omega /\pi )`$ and $`0`$ elsewhere.

Even though the proposition is stated for monotonic densities only, let us apply it to our unstable density (4.9). In this case $`K=7/24`$ and $`T=\sqrt{\frac{16}{7}}1.51186`$. The experimentally determined period of oscillations of the piston is $`T1.62`$, see \[CL\]. We also simulated the piston trajectory with other unstable densities (5.3) with $`r0`$ and observed that the period of oscillations approached its lower bound 2 given by the above proposition.

Proof of Proposition 6.1. Consider a function

$$G(v):=v\frac{\mathrm{sin}(\omega /v)}{1\mathrm{cos}(\omega /v)}$$

Then equation (6.19) reads

$$C:=_{\omega /2\pi }^{\mathrm{}}\rho ^{}vGdv=0$$

Note that in the interval $`[\omega /2\pi ,\mathrm{})`$ the function $`G(v)`$ is strictly increasing and that $`G(\omega /2\pi )=\mathrm{}`$, $`G(\omega /\pi )=0`$, and $`G(\mathrm{})=\mathrm{}`$.

Introducing a new function $`R(v)=\rho ^{}(v)v0`$ and integrating by parts yields

$$M=_{\omega /2\pi }^{\mathrm{}}R(v)𝑑v,K=_{\omega /2\pi }^{\mathrm{}}R(v)\frac{v^2}{6}𝑑v,C=_{\omega /2\pi }^{\mathrm{}}R(v)G(v)𝑑v$$

It is useful to replace $`v`$ by a new variable $`u=G(v)`$, $`\mathrm{}<u<\mathrm{}`$. Since $`G`$ is strictly increasing, we can write

$$M=_{\mathrm{}}^{\mathrm{}}S(u)𝑑u,K=_{\mathrm{}}^{\mathrm{}}S(u)\eta (u)𝑑u,C=_{\mathrm{}}^{\mathrm{}}S(u)u𝑑u$$

where $`S(u)=R(G^1(u))/G^{}(G^1(u))`$ and $`\eta (u)=(G^1(u))^2/6`$ (here $`G^1`$ denotes the inverse of the function $`G`$).

We have to solve the following variational problem: minimize $`K/M`$ under the constraint $`C=0.`$ As we shall see in the sequel the function $`\eta `$ turns out to be convex. This easily implies Proposition 6.1.

Also, the convexity of $`\eta `$ implies that the solution of the variational problem $`\overline{S}`$ is a delta-function centered at $`u=0,`$ i.e. at $`v=\omega /\pi `$ (so that $`C`$ vanished).

So it only remains to prove that $`\eta `$ is a convex function of $`u.`$ By direct computation we get

$$6\eta ^{}=\frac{d}{du}(G^1(u))^2=2G^1(u)\frac{dG^1(u)}{du}=2v/G^{}(v)$$

Hence it is sufficient to prove that the function $`G^{}(v)/v`$ is strictly decreasing in the interval $`v>\omega /2\pi .`$ Without loss of generality we set $`\omega =1`$, then

$$G(v)=\frac{v\mathrm{sin}(1/v)}{1\mathrm{cos}(1/v)}$$

Consider a new function

$$H(v):=\frac{G^{}(v)}{v}=\frac{1+v\mathrm{sin}(1/v)}{v^2(1\mathrm{cos}(1/v))}$$

then

$$H(v)^{}=\frac{v(1+v^2)\mathrm{sin}(1/v)+v\mathrm{cos}(1/v)(1+v\mathrm{sin}(1/v))}{v^4(1\mathrm{cos}(1/v))^2}$$

The denominator of $`H^{}`$ being positive, we only need to prove that the numerator of $`H^{}`$ is negative in the interval $`v>1/2\pi .`$

If we replace $`v`$ by $`1/x`$ and multiply by the numerator by $`x^2`$ we find the expression

$`h(x)`$ $`=`$ $`x+(x^21)\mathrm{sin}(x)+\mathrm{cos}(x)(x+\mathrm{sin}x)`$
$`=`$ $`(\mathrm{cos}x1)(\mathrm{sin}x+x)+x^2\mathrm{sin}x`$
We need to show that $`h(x)<0`$ in the interval $`x(0,2\pi ).`$ First of all, $`h(0)=h(2\pi )=0`$ and for any $`x(\pi ,2\pi )`$ the expression is clearly negative.

It only remains to prove that $`h(x)`$ is negative in $`(0,\pi ].`$ By computing the Taylor expansion of $`h`$ about $`h=0`$ one finds

$$h(x)=\underset{k=3}{\overset{+\mathrm{}}{}}(1)^k\frac{2^{2k}4k^2}{(2k+1)!}x^{2k+1}$$

It is easy to prove that for any $`x(0,\pi ]`$ this is an alternating series, the absolute values of its terms being strictly decreasing.

The first few terms of the above expansion are

$$h(x)=\frac{x^7}{180}\left(1\frac{2}{21}x^2+\frac{1}{240}x^4\frac{19}{166320}x^6\right)+𝒪(x^{15})$$

Therefore

$$\frac{x^7}{180}<h(x)<\frac{x^7}{180}\left(1\frac{2}{21}x^2\right)$$

which implies that $`h(x)<0`$ for any $`x(0,\pi ].`$ $`\mathrm{}`$

Acknowledgements. We thank E. Presutti and Ya. Sinai for their contributions to the early stages of this work, C. Marchioro, M. Pulvirenti, and N. Simanyi for useful discussions. We also thank Weinan E, Rick Falk and Michael Vogelius for helpful suggestions about the numerical solution of the hydrodynamical equations. E. Caglioti was partially supported by MIUR GNFM and by the European network HYKE, contract HPRN-CT-2002-00282. N. Chernov was partially supported by NSF grant DMS-0098788. J. Lebowitz was partially supported by NSF grant DMR-9813268 and by Air Force grant F49620-01-0154. This work was started when the authors stayed at the Institute for Advanced Study with partial support by NSF grant DMS-9729992.