text:
# On a slow drift of a massive piston in an ideal gas that remains at mechanical equilibrium

## 1 Introduction

Consider an isolated cylinder filled with an ideal gas and divided into two compartments by a large piston which is free to move along the axis of the cylinder, Fig. 1. The piston interacts with the gas atoms via elastic collisions. Assume that the gas in each compartment separately is at equilibrium with temperature and density $`T_{},n_{}`$ and $`T_+,n_+`$, respectively. Let the gases exert equal pressures on the piston, i.e. let

$$P_{}=n_{}k_BT_{}=n_+k_BT_+=P_+$$

Then the system is at the so called mechanical equilibrium, and according to the laws of thermodynamics this state should be (macroscopically) stable.

However, the system as a whole is not in a true equilibrium state, unless $`T_{}=T_+`$, hence microscopically it is not stable yet and should find ways to evolve to a true, thermal equilibrium, in which $`T_{}=T_+`$. This was predicted earlier by Landau and Lifshitz \[LL\], Feynman \[F\] and others. Recently Gruber, Piasecki and Frachebourg and others \[GP, GF, P\] derived heuristically, by means of kinetic theory, exact formulas describing the slow drift of the piston and the slow heat transfer from the hotter gas to the cooler gas.

Our goal is to derive rigorously the main formulas of \[GP, GF, P\] describing the slow drift of the piston and the slow heat transfer between the gases.

The piston model trivially reduces to a one-dimensional system by the projection onto the axis of the cylinder. Then one obtains an ideal gas on an interval, and the piston itself becomes a heavy point particle. The motion of a heavy particle in an infinite ideal gas of light particles is a classical example of Brownian motion studied by van Kampen \[vK\] and many others \[L, H, DGL, GF, GP, P\].

We consider a one-dimensional ideal gas on the entire line, without boundaries. The heavy piston is initially placed at the origin $`X(0)=0`$ and is at rest $`V(0)=0`$. The initial configuration of gas atoms and their velocities is chosen at random as a realization of a (two-dimensional) Poisson process on the $`(x,v)`$-plane with density $`p(x,v)`$. This means that for any domain $`D\text{I}\text{R}^2`$ the number $`N_D`$ of gas particles $`(x,v)D`$ at time $`t=0`$ is a Poisson random variable with parameter

$$\lambda _D=_Dp(x,v)𝑑x𝑑v$$

The system evolves according to the rules of elastic collisions. Denote the mass of the piston by $`M`$ and the mass of an atom by $`m`$. Since atoms have identical masses, their collisions can be ignored. When an atom with velocity $`v`$ collides with the piston, whose velocity is $`V`$, their velocities after the collision are given by

$$V^{}=\frac{Mm}{M+m}V+\frac{2m}{M+m}v$$
(1.1)
$$v^{}=\frac{Mm}{M+m}v+\frac{2M}{M+m}V$$
(1.2)
These rules preserve the total kinetic energy and the total momentum. Between collisions, all the particles and the piston move with constant velocities.

The position $`X(t)`$ and velocity $`V(t)=dX(t)/dt`$ of the piston make a random process whose characteristics are determined by the initial gas density $`p(x,v)`$. It is natural to assume that $`p(x,v)`$ is symmetric in $`v`$ and spatially homogeneous, i.e. $`p(x,v)=p(v)`$ and $`p(v)=p(v)`$. We can set $`p(v)=nf(v)`$, where $`f(v)=f(v)`$ is a probability density and $`n>0`$ is a (constant) spatial density. Then one can approximate $`X(t)`$ and $`V(t)`$ by certain Gaussian stochastic processes:

###### Theorem 1.1 (Holley \[H\])

Let the density $`f(v)`$ have a finite fourth moment $`v^4f(v)𝑑v<\mathrm{}`$. Then for every finite $`t_0<\mathrm{}`$, the function $`V(t)\sqrt{M}`$ on the interval $`[0,t_0]`$ converges, in distribution, as $`M,n\mathrm{}`$ and $`M/n\mathrm{const}`$, to an Ornstein-Uhlenbeck velocity process $`𝒱_t`$, while $`X(t)\sqrt{M}`$ converges to an Ornstein-Uhlenbeck position process $`𝒳_t`$.

An Ornstein-Uhlenbeck process $`(𝒳_t,𝒱_t)`$ is defined by \[Ne\]

$$d𝒳_t=𝒱_tdt,d𝒱_t=a𝒱_tdt+\sqrt{D}d𝒲_t$$

where $`a>0`$, $`D>0`$ are constants and $`𝒲_t`$ a Wiener process. The Ornstein-Uhlenbeck position process $`𝒳_t`$ converges in an appropriate limit (e.g. $`a\mathrm{}`$, $`a^2/D=`$const) to a Wiener process.

We note that Dürr et al \[DGL\] extended the above theorem to arbitrary dimensions.

Our paper concerns with another physically interesting situation, where the initial density is not spatially homogeneous, but a spatial homogeneity is assumed separately for the gases to the right and to the left of the piston. So we assume that

$$p(x,v)=\{\begin{array}{cc}p_+(v)& \mathrm{for}x>0\\ p_{}(v)& \mathrm{for}x<0\end{array}$$

and we also assume symmetry: $`p_{}(v)=p_{}(v)`$ and $`p_+(v)=p_+(v)`$. We can set $`p_\pm (v)=n_\pm f_\pm (v)`$, where $`f_\pm (v)=f_\pm (v)`$ are probability densities and $`n_\pm >0`$ are (constant) spatial densities.

In addition, we want to exclude macroscopic motion of the piston in either direction. The macroscopic velocity of the piston vanishes if and only if

$$n_{}_0^{\mathrm{}}v^2f_{}(v)𝑑v=n_+_0^{\mathrm{}}v^2f_+(v)𝑑v$$
(1.3)
as it was shown heuristically in \[LPS\] and under some conditions rigorously in \[CLS\]. The physical interpretation of Eq. (1.3) is the pressure balance. If we define the pressures of the gases by

$$P_\pm =mn_\pm _{\mathrm{}}^{\mathrm{}}v^2f_\pm (v)𝑑v$$

(which would be a proper thermodynamical pressure if the velocity distributions were Maxwellian), then the condition (1.3) means exactly that $`P_{}=P_+`$.

To simplify some technical considerations we assume that the initial densities of the gases satisfy the following velocity cutoff:

$$p_\pm (v)=0,\mathrm{if}|v|v_{\mathrm{min}}\mathrm{or}|v|v_{\mathrm{max}}$$
(1.4)
for some $`0<v_{\mathrm{min}}<v_{\mathrm{max}}<\mathrm{}`$. Hence, the initial velocities of atoms are bounded away from zero and infinity. Under these conditions, our arguments are rigorous. We also discuss in Section 4 how to relax these assumptions, leaving this work for the future.

## 2 Markov approximation

Our cutoff condition (1.4) has an important implication. As long as the speed of the piston remains small enough, every gas atom collides with the piston at most once. Indeed, let $`|V(t)|<V_{\mathrm{max}}`$, where

$$V_{\mathrm{max}}=\frac{Mm}{3M+m}v_{\mathrm{min}}$$

Note that $`V_{\mathrm{max}}`$ is close to $`v_{\mathrm{min}}/3`$ when $`Mm`$. Then it follows directly from (1.4) and (1.2) that the atoms’ velocities after collisions are at least

$$\frac{Mm}{M+m}v_{\mathrm{min}}\frac{2M}{M+m}V_{\mathrm{max}}=V_{\mathrm{max}}$$

the last equality holds due to our choice of $`V_{\mathrm{max}}`$. Therefore, the atoms after collisions remain faster than the piston, so the latter cannot “catch up” with them.

Hence, as long as $`|V(t)|<V_{\mathrm{max}}`$, the velocity of the piston $`V(t)`$ evolves as a Markov process with piecewise constant trajectories (a jump or step process). Moreover, this is a stationary (homogeneous in time) Markov process due to our requirements on the density $`p(x,v)`$. We will slightly change the function $`V(t)`$ so that it will evolve as a stationary Markov process unconditionally. Fix some $`\overline{V}(0,V_{\mathrm{max}})`$ and require that whenever the velocity $`V^{}`$ of the piston after a collision, in the notation of (1.1), exceeds $`\overline{V}`$ in absolute value, it gets “reflected” at $`\overline{V}`$, i.e. it instantaneously changes from $`V^{}`$ to $`V^{\prime \prime }`$ by the rules

$$\begin{array}{ccc}V^{}>+\overline{V}\hfill & & V^{\prime \prime }=+2\overline{V}V^{}\hfill \\ V^{}<\overline{V}\hfill & & V^{\prime \prime }=2\overline{V}V^{}\hfill \end{array}$$
(2.1)
These rules are in the spirit of random walks with reflecting boundary conditions.

We will denote the velocity of the piston in the so defined dynamics by $`W(t)`$. It is clear that $`W(t)=V(t)`$ for all $`t<T`$ such that $`sup_{0<t<T}|V(t)|<\overline{V}`$. We first study the Markov process $`W(t)`$ and later estimate the difference $`V(t)W(t)`$.

Denote by $`𝒫(u,dw;\mathrm{\Delta }t)`$ for every $`\mathrm{\Delta }t>0`$ the transition probability for the process $`W(t)`$, i.e.

$$P\left(W(t+\mathrm{\Delta }t)A/W(t)=u\right)=_A𝒫(u,dw;\mathrm{\Delta }t)$$

for every Borel set $`A\text{I}\text{R}`$. It is clear that the piston does not experience any collisions with the atoms during the interval $`(t,t+\mathrm{\Delta }t)`$ if and only if the trapezoidal domain

$$D=D(t,\mathrm{\Delta }t)=\{(x,v):\frac{vu}{xX(t)}<\frac{1}{\mathrm{\Delta }t},v_{\mathrm{min}}<|v|<v_{\mathrm{max}}\}$$
(2.2)
does not contain gas atoms. Therefore, the probability that $`W(t+\mathrm{\Delta }t)=W(t)=u`$ is

$$𝒫(u,\{u\};\mathrm{\Delta }t)=\mathrm{exp}\left(_Dp(x,v)𝑑x𝑑v\right)$$
(2.3)
To evaluate this integral, we partition the domain $`D=D(t,\mathrm{\Delta }t)`$ into two trapezoids, $`D=D^{}D^+`$ as shown on Fig. 2. Under our assumptions (1.4)

$$_Dp(x,v)𝑑x=\mathrm{\Delta }t\left[_{v_{\mathrm{min}}}^{v_{\mathrm{max}}}(vu)p_{}(v)𝑑v+_{v_{\mathrm{max}}}^{v_{\mathrm{min}}}(uv)p_+(v)𝑑v\right]$$

We introduce the following notation: for each $`k0`$ let

$$_{v_{\mathrm{min}}}^{v_{\mathrm{max}}}v^kp_{}(v)𝑑v=F_k^{}\mathrm{and}_{v_{\mathrm{max}}}^{v_{\mathrm{min}}}v^kp_+(v)𝑑v=F_k^+$$

and

$$Q_k=F_k^{}F_k^+$$

Then we obtain

$`{\displaystyle _D}p(x,v)𝑑x`$ $`=`$ $`\mathrm{\Delta }t\left[F_1^{}F_1^+(F_0^{}F_0^+)u\right]`$ (2.4)
$`=`$ $`\mathrm{\Delta }t\left[Q_1Q_0u\right]`$
where $`u=W(t)`$.

It is clear that for every $`u`$ and $`\mathrm{\Delta }t>0`$ the probability measure $`𝒫(u,dw;\mathrm{\Delta }t)`$, besides having an atom at $`u`$ with probability given by (2.3), has an absolutely continuous component with a positive density on the interval $`(\overline{V},\overline{V})`$ (in fact, that density is bounded away from zero for every fixed $`u`$ and $`\mathrm{\Delta }t>0`$).

###### Proposition 2.1

The stationary Markov process $`W(t)`$ has a unique stationary measure $`\mu _0`$, which is absolutely continuous on the interval $`(\overline{V},\overline{V})`$. Any other initial distribution converges to $`\mu _0`$ uniformly exponentially fast in time.

Proof. The process $`W(t)`$ satisfies the so called Doeblin condition \[D\]: there exist $`\epsilon >0`$ and $`\mathrm{\Delta }t>0`$ such that for every Borel set $`A(\overline{V},\overline{V})`$ and every $`u(\overline{V},\overline{V})`$ we have

$$m(A)<\epsilon _A𝒫(u,dw;\mathrm{\Delta }t)<1\epsilon $$

where $`m(A)`$ is the Lebesgue measure of $`A`$. Then our result follows from general theorems \[D\]. But we also outline a simple direct argument.

Since the transition probability $`𝒫(u,dw;\mathrm{\Delta }t)`$ defines, for each $`\mathrm{\Delta }t>0`$, a continuous map on the convex space of probability distributions on the interval $`(\overline{V},\overline{V})`$, the existence follows from a general Schauder-Tychonoff theorem \[DS\], p. 456. Since $`𝒫(u,dw;\mathrm{\Delta }t)`$ has an absolutely continuous component whose density is bounded away from zero, the stationary distribution $`\mu _0`$ is absolutely continuous with density $`f_0(w)c_0>0`$.

To prove the uniqueness of $`\mu _0`$, suppose $`\mu _0^{}\mu _0`$ is another stationary measure with density $`f_0^{}(w)c_0^{}>0`$. Then $`(1+c)f_0cf_0^{}`$ for small $`c>0`$ is also a stationary density. Let

$$\overline{c}=sup\{c>0:\underset{|w|<\overline{V}}{inf}[(1+c)f_0(w)cf_0^{}(w)]0\}$$

Then $`(1+\overline{c})f_0\overline{c}f_0^{}`$ is a stationary density, hence it must also be bounded away from zero on $`(\overline{V},\overline{V})`$. But this contradicts our choice of $`\overline{c}`$.

To prove the convergence, we take an arbitrary initial distribution $`\nu _0`$ and consider its image $`\nu _t`$ at time $`t`$. Fix a $`t>0`$. The measure $`\nu _t`$ has an absolutely continuous component whose density is bounded away from zero. Hence $`\nu _t=(1\alpha )\nu ^{(1)}+\alpha \mu _0`$ with some small $`\alpha >0`$ and some measure $`\nu ^{(1)}`$. The value of $`\alpha `$ depends on $`t`$ but not on $`\nu _0`$. This implies, by induction, that $`\nu _{kt}=(1\alpha )^k\nu ^{(k)}+[1(1\alpha )^k]\mu _0`$ for all $`k1`$, which proves the exponential convergence of $`\nu _t`$ to $`\mu _0`$. $`\mathrm{}`$

We now derive one useful equation. Fix a $`\mathrm{\Delta }t>0`$ and for every $`u(\overline{V},\overline{V})`$ denote by

$$E(u,\mathrm{\Delta }t)=w𝒫(u,dw;\mathrm{\Delta }t)$$

the conditional expectation of the process $`W(t+\mathrm{\Delta }t)`$ given that $`W(t)=u`$. The invariance of the measure $`\mu _0`$ implies

$$E(u,\mathrm{\Delta }t)𝑑\mu _0(u)=w𝑑\mu _0(w)$$
(2.5)
Now suppose that

$$E(u,\mathrm{\Delta }t)=u+g(u)\mathrm{\Delta }t+R(u,\mathrm{\Delta }t)$$
(2.6)
so that

$$|R(u,\mathrm{\Delta }t)|\mathrm{const}(\mathrm{\Delta }t)^2$$

Substituting this into (2.5) and taking the limit as $`\mathrm{\Delta }t0`$ gives

$$g(u)𝑑\mu _0(u)=0$$
(2.7)
Our next step is to derive (2.6) for small $`\mathrm{\Delta }t>0`$ and compute $`g(u)`$ explicitly. Let $`X(t)=X`$ and $`W(t)=u`$. Due to our assumption (1.4) the probability that the domain

$$D_1=\{(x,v):0<(Xx)\mathrm{sgn}v<\overline{V}\mathrm{\Delta }t,v_{\mathrm{min}}<|v|<v_{\mathrm{max}}\}$$

contains more than one atom is $`𝒪((\mathrm{\Delta }t)^2)`$, hence we can ignore it and assume that $`D_1`$ contains at most one atom. Note that $`D_1`$ is much larger than $`D`$ defined by (2.2). Now it is easy to see that the piston experiences at most one collision during the interval $`(t,t+\mathrm{\Delta }t)`$. Precisely, a collision occurs if and only if the domain $`D=D^+D^{}`$ contains an atom. When it does, and the atom has velocity $`v`$ at time $`t`$, then the piston’s velocity at time $`t+\mathrm{\Delta }t`$ is computed by (1.1):

$$W(t+\mathrm{\Delta }t)=u+\frac{2m}{M+m}(vu)$$
(2.8)
Therefore, the conditional expectation $`E(u,\mathrm{\Delta }t)`$ is

$$E(u,\mathrm{\Delta }t)=u+\frac{2m}{M+m}\left[_Dvp𝑑x𝑑vu_Dp𝑑x𝑑v\right]+𝒪((\mathrm{\Delta }t)^2)$$
(2.9)
(the last term comes from our assumption that $`D`$ contains at most one atom). While the second integral in (2.9) is already computed in (2.4), the calculation of the first one is similar. For every $`k0`$ we have

$`{\displaystyle _D}v^kp(x,v)𝑑x`$ $`=`$ $`\mathrm{\Delta }t\left[{\displaystyle _{v_{\mathrm{min}}}^{v_{\mathrm{max}}}}v^k(vu)p_{}(v)𝑑v+{\displaystyle _{v_{\mathrm{max}}}^{v_{\mathrm{min}}}}v^k(uv)p_+(v)𝑑v\right]`$ (2.10)
$`=`$ $`\mathrm{\Delta }t\left[F_{k+1}^{}F_{k+1}^+(F_k^{}F_k^+)u\right]`$
$`=`$ $`\mathrm{\Delta }t\left[Q_{k+1}Q_ku\right]`$
Thus,

$$E(u,\mathrm{\Delta }t)=u+\frac{2m}{M+m}\left[Q_22Q_1u+Q_0u^2\right]\mathrm{\Delta }t+𝒪((\mathrm{\Delta }t)^2)$$

This gives us

$$g(u)=\frac{2m}{M+m}\left[Q_22Q_1u+Q_0u^2\right]$$

We now apply the key identity (2.7) and arrive at

$$Q_22Q_1u+Q_0u^2=0$$

where $``$ means averaging with respect to the stationary measure $`\mu _0`$ (we note that $`Q_i`$ do not depend on $`u`$). Thus, the average velocity of the piston is given by

$$W=\frac{Q_0W^2+Q_2}{2Q_1}$$
(2.11)
The mechanical equilibrium is defined by the equality of pressures on both sides of the piston, i.e. by

$$Q_2=F_2^{}F_2^+=0$$

and we will restrict ourselves to this case in the rest of the paper. Then (2.11) reduces to

$$W=\frac{Q_0W^2}{2Q_1}$$
(2.12)
which is still different from zero, albeit very small.

It is well known that at mechanical equilibrium the average value of $`W^2`$ is $`𝒪(m/M)`$, see, e.g. \[H, DGL\]. But we will estimate it more accurately below.

Similarly to $`E(u,\mathrm{\Delta }t)`$, we define

$$E_2(u,\mathrm{\Delta }t)=w^2𝒫(u,dw;\mathrm{\Delta }t)$$

the conditional expectation of the square velocity $`W^2(t+\mathrm{\Delta }t)`$ given that $`W(t)=u`$. The stationarity of the distribution $`\mu _0`$ implies

$$E_2(u,\mathrm{\Delta }t)𝑑\mu _0(u)=u^2𝑑\mu _0(u)$$

Assume that

$$E_2(u,\mathrm{\Delta }t)=u^2+g_2(u)\mathrm{\Delta }t+R_2(u,\mathrm{\Delta }t)$$

with $`|R_2(u,\mathrm{\Delta }t)|\mathrm{const}(\mathrm{\Delta }t)^2`$. Then, as in (2.7)

$$g_2(u)𝑑\mu _0(u)=0$$
(2.13)
Now, squaring the equation (2.8) we obtain

$$W^2(t+\mathrm{\Delta }t)=u^2+\frac{4m}{M+m}(vu)u+\frac{4m^2}{(M+m)^2}(vu)^2$$

Then, as in (2.9),

$`E_2(u,\mathrm{\Delta }t)=u^2`$ $`+`$ $`{\displaystyle \frac{4m}{M+m}}{\displaystyle _D}(vu)up𝑑x𝑑v`$
$`+`$ $`{\displaystyle \frac{4m^2}{(M+m)^2}}{\displaystyle _D}(vu)^2p𝑑x𝑑v+𝒪((\mathrm{\Delta }t)^2)`$
Computing these integrals as before gives

$`g_2(u)`$ $`=`$ $`{\displaystyle \frac{4m}{M+m}}(Q_2u2Q_1u^2+Q_0u^3)`$
$`+{\displaystyle \frac{4m^2}{(M+m)^2}}(Q_33Q_2u+3Q_1u^2Q_0u^3)`$
We included the terms with $`Q_2`$ for the sake of clarity, even though we had assumed that $`Q_2=0`$. In the subsequent expressions we remove all $`Q_2`$’s. The equation (2.13) gives

$$2Q_1W^2Q_0W^3=\frac{m}{M+m}(Q_3+3Q_1W^2Q_0W^3)$$

For brevity, we denote $`\epsilon =m/M`$. By making one step further and analyzing $`E_3(u,\mathrm{\Delta }t)=w^3𝒫(u,dw;\mathrm{\Delta }t)`$ in the same manner one can see that $`W^3=𝒪(\epsilon ^2)`$, we omit details. Hence we arrive at

$$W^2=\frac{Q_3}{2Q_1}\epsilon +𝒪(\epsilon ^2)$$
(2.14)
Substituting this into (2.12) yields

$$W=\frac{Q_0Q_3}{4Q_1^2}\epsilon +𝒪(\epsilon ^2)$$
(2.15)
## 3 Thermodynamical analysis

In this section we express our rigorous results in statistical mechanical terms, such as temperature, density, pressure, and heat transfer.

First, we consider the heat flow across the piston, from right to left. Even though both gases are infinite, the heat transfer from one to the other is well defined and its rate can be computed. The thermodynamical definition of the heat transfer \[GF\] is

$$R_+=d_{}/dtd_{}/dtW$$
(3.1)
where $`d_{}/dt`$ is the rate of change of kinetic energy and $`d_{}/dt`$ the rate of change of momentum on the left hand side. By the conservation of energy and momentum during collisions, it is enough to compute the average change of energy and momentum of the piston when it collides with an atom on the left hand side.

According to (1.1), the momentum of the piston at one collision changes by

$$\delta _M=\frac{2mM}{M+m}(vW)$$

The conditional average change of momentum during the interval $`(t,t+\mathrm{\Delta }t)`$, given that $`W(t)=u`$, is

$`E_u(\delta _M)`$ $`=`$ $`{\displaystyle \frac{2mM}{M+m}}{\displaystyle _D^{}}(vu)𝑑x𝑑v+𝒪((\mathrm{\Delta }t)^2)`$
$`=`$ $`{\displaystyle \frac{2mM}{M+m}}(F_2^{}2F_1^{}u+F_0^{}u^2)\mathrm{\Delta }t+𝒪((\mathrm{\Delta }t)^2)`$
where only collisions with atoms on the left are taken into account. Negating and averaging with respect to the stationary measure $`\mu _0`$ gives the average momentum transfer to the left gas:

$$d_{}/dt=2m[F_2^{}2F_1^{}W+F_0^{}W^2]+𝒪(m^2/M)$$
(3.2)
By squaring (1.1) we find the change of the kinetic energy of the piston at one collision:

$$\delta _M=2m\epsilon v^2+2mVv2mV^2+\mathrm{}$$

where only essential terms are shown (the rest will eventually end up in $`𝒪(m^3/M^2)`$ in the expression (3.3) below, so we omit them). The conditional average change of kinetic energy during the interval $`(t,t+\mathrm{\Delta }t)`$, given that $`W(t)=u`$, is (again only collisions with atoms on the left are taken into account)

$$E_u(\delta _M)=2m(\epsilon F_3^{}\epsilon F_2^{}u+F_2^{}u2F_1^{}u^2+F_0^{}u^3)\mathrm{\Delta }t+\mathrm{}+𝒪((\mathrm{\Delta }t)^2)$$

Negating and averaging with respect to the stationary measure $`\mu _0`$ gives the average kinetic energy transfer to the left gas:

$$d_{}/dt=2m[\epsilon F_3^{}+F_2^{}W2F_1^{}W^2]+𝒪(m^3/M^2)$$
(3.3)
Combining (3.2) and (3.3) we obtain the expression for the heat transfer

$$R_+=2m[\epsilon F_3^{}2F_1^{}W^2]+𝒪(\epsilon ^2)$$
(3.4)
Substituting (2.14) gives

$$R_+=\frac{2m\epsilon }{Q_1}(F_1^+F_3^{}F_1^{}F_3^+)+𝒪(\epsilon ^2)$$
(3.5)
Now, even though the thermodynamical temperatures of the gases are not defined, unless their velocity distributions are Maxwellian, we can define the “effective temperatures” via the second moment of the velocity distributions:

$$T_\pm =k_B^1m\frac{v^2p_\pm 𝑑v}{p_\pm 𝑑v}=k_B^1m\frac{F_2^\pm }{F_0^\pm }$$

where $`k_B`$ is Boltzmann’s constant. We recall that the pressures of the gases are given by $`P_\pm =2mF_2^\pm `$ and their spatial densities by $`n_\pm =2F_0^\pm `$, hence the classical law $`P_\pm =n_\pm k_BT_\pm `$ holds.

The piston’s velocity distribution $`\mu _0`$ is, generally, not Maxwellian, but its effective temperature can be computed similarly:

$$T_M=k_B^1M(W^2W^2)k_B^1m\frac{Q_3}{2Q_1}$$

where we used (2.14) and (2.15). Hence $`T_M`$ is of the same order of magnitude as $`T_\pm `$, but for generic distributions $`p_\pm `$ there is no simple relation between the exact values of these temperatures.

However, we can derive some interesting relations assuming that the distributions $`p_\pm `$ are Maxwellian. Since this would violate our assumptions (1.4), our further analysis is heuristic. Assume that the densities $`p_\pm `$ satisfy

$$p_\pm (x,v)=n_\pm f_\pm (v)$$

where $`n_\pm `$ are (constant) spatial densities and

$$f_\pm (v)=\frac{1}{\sqrt{2\pi \sigma _\pm ^2}}\mathrm{exp}\left(\frac{v^2}{2\sigma _\pm ^2}\right)$$

the Maxwellian distributions with temperatures $`T_\pm =k_B^1m\sigma _\pm ^2`$. As before, we assume that the pressures are equal, i.e. $`n_{}T_{}=n_+T_+`$. It is then easy to compute

$$Q_0=\frac{1}{2}(n_{}n_+)$$

$$Q_1=\frac{1}{\sqrt{2\pi }}(n_{}\sigma _{}+n_+\sigma _+)=\frac{\sqrt{k_B}}{\sqrt{2\pi m}}\left(n_{}T_{}^{1/2}+n_+T_+^{1/2}\right)$$

$$Q_2=\frac{1}{2}(n_{}\sigma _{}^2n_+\sigma _+^2)=\frac{k_B}{2m}(n_{}T_{}n_+T_+)=0$$

$$Q_3=\sqrt{\frac{2}{\pi }}(n_{}\sigma _{}^3+n_+\sigma _+^3)=\frac{\sqrt{2k_B^3}}{\sqrt{\pi m^3}}\left(n_{}T_{}^{3/2}+n_+T_+^{3/2}\right)$$

Now the effective temperature of the piston is

$$T_M=\frac{mQ_3}{2k_BQ_1}=\frac{n_{}T_{}^{3/2}+n_+T_+^{3/2}}{n_{}T_{}^{1/2}+n_+T_+^{1/2}}=\sqrt{T_{}T_+}$$
(3.6)
where we used the equality of the pressures. The average velocity of the piston is

$`W`$ $`=`$ $`{\displaystyle \frac{\sqrt{2\pi k_Bm}}{4M}}{\displaystyle \frac{(n_{}n_+)(T_{}T_+)^{1/2}}{n_{}T_{}^{1/2}+n_+T_+^{1/2}}}+𝒪(\epsilon ^2)`$ (3.7)
$`=`$ $`{\displaystyle \frac{\sqrt{2\pi k_Bm}}{4M}}\left(\sqrt{T_+}\sqrt{T_{}}\right)+𝒪(\epsilon ^2)`$
Lastly, the heat transfer across the piston (from right to left) given by (3.5) is now

$`R_+`$ $`=`$ $`{\displaystyle \frac{2m\epsilon \sqrt{2k_B^3}}{\sqrt{\pi m^3}}}\times {\displaystyle \frac{n_{}n_+(T_+^{3/2}T_{}^{1/2}T_{}^{3/2}T_+^{1/2})}{n_{}T_{}^{1/2}+n_+T_+^{1/2}}}+𝒪(\epsilon ^2)`$ (3.8)
$`=`$ $`{\displaystyle \frac{k_B}{M}}\sqrt{{\displaystyle \frac{8k_Bm}{\pi }}}{\displaystyle \frac{n_{}n_+(T_{}T_+)^{1/2}}{n_{}T_{}^{1/2}+n_+T_+^{1/2}}}(T_+T_{})+𝒪(\epsilon ^2)`$
Hence, the heat flow is proportional to the temperature gradient, $`T_+T_{}`$, and the heat conductivity is thus

$$\kappa =\frac{k_B}{M}\sqrt{\frac{8k_Bm}{\pi }}\frac{n_{}n_+(T_{}T_+)^{1/2}}{n_{}T_{}^{1/2}+n_+T_+^{1/2}}$$
(3.9)
The formulas (3.6)–(3.9) were obtained earlier in \[GP, GF\] by means of kinetic theory.

Let us examine closely what happens when the pressures are equal, $`P_{}=P_+`$, but the temperatures are different. Without loss of generality, assume that $`T_{}<T_+`$, i.e. the gas on the left is cooler but denser and the gas on the right is hotter but sparser. Then $`W>0`$, i.e. the piston slowly drifts in the direction of the hotter gas. At the same time $`R_+>0`$, i.e. the heat slowly flows from the hotter gas to the cooler gas. This demonstrates that the mechanical equilibrium is not a stable state. If our gases were finite, there would be a continuing evolution toward thermal equilibrium, in which the temperatures become equal, too. Of course, this conclusion makes no sense in the idealized model of infinite gases.

We note that the evolution of the system is somewhat counterintuitive. When the piston moves to the right, i.e. $`V(t)>0`$, then by (1.2) the atoms on the right bounce off it with a higher speed and so gain energy, while the atoms on the left collide with the piston and slow down, hence lose some energy. When the piston moves to the left, i.e. $`V(t)<0`$, it is vice versa. Since, on the average, the heat flows from right to left, one may conclude that the piston’s movements to the left dominate. On the other hand, the piston slowly drifts to the right, so its displacements in that direction actually dominate.

To explain this “paradox”, we first recall that the average speed of the piston is of order $`W^2^{1/2}\epsilon ^{1/2}`$, which is much larger than the speed of the drift $`W\epsilon `$. Hence the piston jiggles back and forth much faster than it drifts in one direction. The heat flow is due to “jiggling” rather than “drifting”. Indeed, when the piston jiggles to the right, the slow atoms on the left have less chance to collide with the piston, since the relative velocity is small. Most of the collisions between the slow atoms on the left and the piston occur when the latter jiggles to the left. This is why the cooler atoms speed up and gain energy, on the average. On the other hand, the hotter atoms on the right are less sensitive to the variations in the piston’s velocity, and the reason why they cool down is different, as we explain next.

The piston’s vibrations are not spatially symmetric. During excursions to the right, relatively few collisions with atoms on the right occur, since they are more energetic and able to quickly reverse the velocity of the piston. On the contrary, when the piston drifts to the left, the atoms on the left are weak and it takes them longer to turn the piston back to the right. During these intervals, the fast atoms on the right continue hitting the piston and lose energy. This explains why the hotter atoms mostly slow down. We must admit though that the precise mechanism of the heat flow across the piston in our model remains unclear. Some physicists call it a “conspiracy” between the microscopic vibrations of the piston and the incoming atoms of the gases \[GF, GP\].

## 4 Justification of Markov approximation

Here we estimate the difference between the true velocity of the piston $`V(t)`$ and our Markov approximation $`W(t)`$.

Recall that $`W^2\epsilon =m/M`$, i.e. the typical speed of the piston (in the Markov approximation) is of order $`\sqrt{\epsilon }`$. The following theorem, which is essentially proved in \[CLS\] (see remarks below), estimates large deviations of the true velocity $`V(t)`$:

###### Theorem 4.1 (\[CLS\])

Let $`|V(0)|<\epsilon ^{1/2}`$. There are constants $`C,d>0`$ such that for every $`T>0`$ we have

$$P\left(\underset{0<t<T}{sup}|V(t)|>C\epsilon ^{1/2}\mathrm{ln}\epsilon ^1\right)<T^2\epsilon ^{d\mathrm{ln}\mathrm{ln}\epsilon }$$

This gives a good bound on large deviations during time intervals of length $`T\epsilon ^A`$ for any large but fixed constant $`A>0`$: it shows that the true velocity $`V(t)`$ of the piston remains $`𝒪(\epsilon ^{1/2}\mathrm{ln}\epsilon ^1)`$ with “overwhelming” probability.

In our construction of the Markov approximation $`W(t)`$ we can choose any small but positive cutoff value $`\overline{V}>0`$, and then assume that $`\epsilon `$ is so small that $`C\epsilon ^{1/2}\mathrm{ln}\epsilon ^1<\overline{V}`$. Then we arrive at

###### Corollary 4.2

For all sufficiently small $`\epsilon `$ and all $`T>0`$ we have

$$P\left(V(t)=W(t)t(0,T)\right)1T^2\epsilon ^{d\mathrm{ln}\mathrm{ln}\epsilon }$$

Therefore, during time intervals of length $`T\epsilon ^A`$ for any fixed constant $`A>0`$, the true velocity $`V(t)`$ of the piston coincides with its Markov approximation $`W(t)`$ with “overwhelming” probability. If $`V(t)`$ and $`W(t)`$ differ, though, then we have an obvious bound: $`|V(t)|v_{\mathrm{max}}`$. So we obtain analogues of our main formulas (2.15) and (2.14):

###### Corollary 4.3

Let $`\epsilon >0`$ be small and $`|V(0)|<\epsilon ^{1/2}`$. Then for any fixed $`A>0`$ and all $`0<t<T=\epsilon ^A`$ we have

$$E(V(t))=\frac{Q_0Q_3}{4Q_1^2}\epsilon +𝒪(\epsilon ^2)$$

and

$$E(V^2(t))=\frac{Q_3}{2Q_1}\epsilon +𝒪(\epsilon ^2)$$

where $`E()`$ is the mean value.

It follows from the results of \[CLS\] that whenever the true velocity $`V(t)`$ exceeds $`C\epsilon ^{1/2}\mathrm{ln}\epsilon ^1`$, then with “overwhelming” probability it will be driven back toward zero by further collisions with the gas atoms. Hence, the estimates of the last corollary should remain true for all $`t>0`$, but we do not pursue such a goal, because the drift of the piston becomes evident on the time scale of order $`\epsilon ^A`$ already for $`A>2`$, as we show next, so the estimates in Corollary 4.3 are sufficient for practical purposes.

Indeed, let $`X(0)=0`$ and consider the random position of the piston $`X(t)`$ at time $`t`$. Also, let $`Y(t)=_0^tW(s)𝑑s`$ be the corresponding function in the Markov approximation model. We obviously have

$$Y(t)=tW=\frac{Q_0Q_3}{4Q_1^2}\epsilon t+𝒪(\epsilon ^2t)$$
(4.1)
We also need to estimate the variance of $`Y(t)`$, which is a little harder. We claim that

$$\mathrm{Var}[Y(t)]Dt$$
(4.2)
with some constant $`D>0`$ independent of $`\epsilon `$. Below we outline the argument, suppressing some technical details that are based on the estimates obtained in \[CLS\].

It is standard that

$$\mathrm{Var}[Y(t)]=_0^t_0^t\mathrm{Cov}(W(s),W(u))𝑑s𝑑u$$

By the stationarity of the process $`W(t)`$, the covariance here equals

$$\mathrm{Cov}(W(s),W(u))=\rho _{su}\mathrm{Var}(W)$$

where $`\rho _{su}`$ is the correlation between $`W(s)`$ and $`W(u)`$, which only depends on $`|su|`$. Also recall that

$$\mathrm{Var}(W)=W^2W^2=\frac{Q_3}{2Q_1}\epsilon +𝒪(\epsilon ^2)$$

Therefore,

$$\mathrm{Var}[Y(t)]\mathrm{const}t\epsilon _0^t\rho _s𝑑s$$
(4.3)
It remains to estimate the decay of the velocity autocorrelation function $`\rho _s`$.

Suppose during the interval $`(u,u+s)`$ the piston experiences $`k`$ collisions with gas atoms with velocities $`v_1,\mathrm{},v_k`$ numbered in the order in which the collisions occur. Then (1.1) implies

$$W(u+s)=(1\alpha )^kW(u)+\alpha \underset{j=1}{\overset{k}{}}(1\alpha )^{kj}v_j$$
(4.4)
where $`\alpha =2m/(M+m)2\epsilon `$. This equation can be easily verified by induction on $`k`$.

Now, the number of collisions $`k`$ grows as const$`s`$, on the average, in fact

$$P(|kc_1s|c_2\sqrt{s}\mathrm{ln}s)s^2\epsilon ^{d\mathrm{ln}\mathrm{ln}\epsilon }$$
(4.5)
with some positive constants $`c_1,c_2`$, see \[CLS\]. For $`s<\epsilon ^A`$ this gives a good bound on large deviations, so we can safely assume that $`kc_1\epsilon /2`$. Next, the velocities $`v_i`$ are almost independent of $`W(u)`$. In fact, there is a dependence, since the time of collision between a particular atom and the piston depends on the velocity of the latter, but the correlation between $`W(u)`$ and velocities $`v_j`$ for every fixed $`j`$ is negative. This follows from a simple observation: if $`W(u)>0`$, the piston is more likely to collide with atoms on the right, for which $`v_j<0`$, and vice versa. Therefore, the correlation $`\rho _s`$ is determined by the first term in (4.4), hence

$$|\rho _s|(1\alpha )^{c_1s/2}$$

This implies

$$_0^t\rho _s𝑑s_0^{\mathrm{}}(1\alpha )^{c_1s/2}𝑑s=\frac{2}{c_1\mathrm{ln}(1\alpha )}\frac{1}{c_1\epsilon }$$

and combining with (4.3) proves (4.2).

The estimates (4.1) and (4.2) show that the piston’s displacement due to the drift exceeds its typical random fluctuations at times $`t=\epsilon ^A`$ for all $`A>2`$. On such time scales the drift is “physically evident” and it is possible to verify it experimentally. Numerical experiments of this sort were done by \[GF\].

We now comment on the proofs of Theorem 4.1 and the estimate (4.5). The corresponding results were proven in \[CLS\] in a slightly different context. That paper dealt with a heavy piston in a cubical container of side $`L`$ filled by an ideal gas. That model, too, reduced by a trivial projection onto the axis perpendicular to the piston surface to a one-dimensional gas on an interval $`[L/2,L/2]`$. The velocity cutoff assumptions in \[CLS\] were identical to our (1.4). The piston had mass $`M=aL^2`$ ($`a>0`$ was a constant) and the spatial density of the one-dimensional gases was proportional to $`M=aL^2`$, but a simple rescaling of time and space by $`M`$ reduces the spatial density to values $`𝒪(1)`$ and makes it independent of $`M`$. Hence, we arrive exactly at the same initial conditions as in the present paper, except the gases in \[CLS\] were finite – they evolved on an interval $`[aL^3/2,aL^3/2]`$ with the piston initially placed at $`X(0)=0`$.

Due to the velocity cutoff (1.4) the gas atoms in \[CLS\] cannot interact with the piston more than once unless they travel across the half interval $`[aL^3/2,aL^3/2]`$, reflect at an endpoint $`x=\pm aL^3/2`$ and then travel back to the piston in the middle. This takes time $`T=𝒪(L^3)`$, as it was proven in \[CLS\]. Thus, during the initial interval $`(0,T)`$ with $`T=𝒪(M^{3/2})`$ the evolution of the piston studied in \[CLS\] is identical to ours. Hence, all the results and estimates derived in Section 3 of \[CLS\] directly apply to our present context. Furthermore, with some minor obvious changes those results can be extended to arbitrary time intervals as we did in Theorem 4.1 and (4.5).

Lastly, we have assumed the velocity cutoff (1.4) here and in the paper \[CLS\] mostly for convenience. It might be true that the upper bound on velocities is quite essential in \[CLS\], since arbitrarily fast atoms could bounce back and forth between the piston and the wall many times during short time intervals thus making significant impact on the piston, as it was remarked in Section 5 of \[CLS\]. But in the present work, if we had arbitrarily fast particles, they would hit the piston once and get away never to come back, so their impact would be quite limited.

As for the lower bound on velocities, it seems to be less essential in both models. It was conjectured in Section 5 of \[CLS\] that the results of that paper could be extended to velocity distributions without lower bounds (i.e. with $`v_{\mathrm{min}}=0`$). If so, then the estimate in Theorem 4.1 implies that the typical velocity of the piston remains of order $`\epsilon ^{1/2}\mathrm{ln}\epsilon ^1`$. As long as this is true, only atoms with smaller velocities can experience two or more collisions with the piston, thus violating the Markovian character of the piston’s velocity process $`V(t)`$. It is customary to redefine the dynamics so that the collisions with slow atoms are ignored altogether, and then the resulting process $`W(t)`$ will be Markovian. The difference $`|V(t)W(t)|`$ should be estimated separately and presumably is negligible. Such a strategy was successfully implemented in earlier works \[H, DGL\] in somewhat different context, and it is likely to go through in our case as well. This is the subject of our future work.

Acknowledgement. The author thanks J. Lebowitz, J. Piasecki and Ya. Sinai for many useful discussions and encouragement. The author was partially supported by NSF grant DMS-0098788. This work was completed during the author’s stay at the Institute for Advanced Study with partial support by NSF grant DMS-9729992.


format: text

text:
# Non–equilibrium CO chemistry in the solar atmosphere

## 1 Introduction

Thirty years ago Noyes & Hall (1972) inferred very low brightness temperatures (T$`\mathrm{\hspace{0.17em}3700}`$ K) from their discovery of strong rovibrational CO lines at 4.7 $`\mu \mathrm{m}`$ observed close to the edge of the solar disk. Several years later, after such surprising observational results had been confirmed by Ayres & Testerman (1981), it was suggested that the low chromosphere might not be hot at all but could instead be permeated by CO-cooled ‘clouds’ at altitudes between 500 and 1000 kilometers above continuum optical depth unity (Ayres 1981). This led to controversy because other (UV and submillimeter) diagnostics had suggested the existence of a uniformly hot chromosphere with a minimum temperature of about 4400 K near 500 km and a temperature rise above this temperature-minimum region. The controversy over the existence of cool gas in the low chromosphere continues today (see Kalkofen 2001, Ayres 2002, Avrett 2003), after the publication of an abundance of literature on the subject with theoretical investigations concluding that the CO lines have LTE source functions (Ayres & Wiedemann 1989; Uitenbroek 2000a) and with new observations showing off-limb CO emission protruding hundreds of kilometers into the chromosphere (Solanki, Livingston & Ayres 1994; Clark et al. 1995), but also with the discovery that far-UV chromospheric lines observed on the solar disk always remain in emission at all positions and times (Carlsson, Judge & Wilhem 1997).

Over the last few years, it has become increasingly evident that the next crucial step towards a better understanding of the enigmatic thermal structure of the solar chromosphere is to carefully investigate the reliability of the assumption of instantaneous chemical equilibrium (ICE), which is currently used for calculating the molecular number densities in stellar atmospheres (see, in particular, Uitenbroek 2000a,b; see also Ayres & Rabin 1996, Avrett et al. 1996, Asensio Ramos & Trujillo Bueno 2003, Avrett 2003). Actually, both observations and simulations indicate that the solar chromosphere is a highly inhomogeneous and dynamic region of low density plasma whose thermal, dynamic and magnetic properies we need to decipher for unlocking new discoveries. If the ICE approximation turns out to be adequate for modeling the strongest CO lines close to the edge of the solar disk, then the available CO observations would really be indicating the existence of cool gas in the solar chromosphere. Otherwise, a natural resolution of the current chromospheric temperature discrepancy could perhaps emerge if the CO concentration turns out to be sufficiently lower than would be expected on the basis of chemical equilibrium. This Letter addresses this challenging issue by showing the first results of a detailed investigation on the non-equilibrium CO chemistry in the solar atmosphere.

## 2 Formulation of the problem

Our strategy consists in performing chemical evolution calculations in the radiation hydrodynamical simulations of solar chromospheric dynamics described by Carlsson & Stein (1997, 2002), that do not include CO cooling in the energy equation. Therefore, at each time step of the hydrodynamical simulation, we have fixed the ensuing thermodynamic conditions and calculated the corresponding CO number density by following the chemical evolution starting from the molecular concentrations of the previous time step.

The chemical evolution (CE) problem consists in calculating the molecular number densities at each time step of the hydrodynamical simulation by taking into account all relevant chemical reactions with their ensuing rates. In order to solve this problem, one has to consider the following evolution equation for each chemical species ‘$`i`$’ included in the model:

$`{\displaystyle \frac{dn_i}{dt}}={\displaystyle \underset{A}{}}{\displaystyle \underset{B}{}}{\displaystyle \underset{C}{}}k_{ABC}n_An_Bn_C+{\displaystyle \underset{A}{}}{\displaystyle \underset{B}{}}k_{AB}n_An_B+`$
$`+{\displaystyle \underset{A}{}}k_An_A{\displaystyle \underset{A}{}}{\displaystyle \underset{B}{}}k_{ABi}n_An_Bn_i`$
$`{\displaystyle \underset{A}{}}k_{Ai}n_An_ik_in_i,`$ (1)
where *three-body* reactions (first and fourth terms), *two-body* reactions (second and fifth terms) and *one-body* reactions (third and sixth terms) have been taken into account. When all the reaction rates for three-body ($`k_{ABC}`$), two-body ($`k_{AB}`$) and one-body reactions ($`k_A`$) are given, this set of equations for all the species included in the model form a set of first order ordinary differential equations which has to be solved numerically. Due to the *stiffness* of the system of equations (produced by the enormous differences in the abundances and their temporal variability), an implicit scheme has to be used. Two key ingredients have to be taken into account: the number of chemical species and the reaction rates for all possible reactions. After a careful investigation, we found that at least the following set of 13 species, which includes the most abundant diatomic molecules, is required to achieve a suitable description of the relevant chemical processes: H, C, O, N, He, CH, CO, H2, OH, NH, N2, NO and CN. We have verified that the inclusion of ionic species does not significantly affect the CO concentration for the typical atmospheric conditions encountered at heights $`h1000`$ km, although they should be ideally taken into account for very strong shocks capable of producing sizable changes in the degree of ionization. Concerning the reactions rates, we have used a reaction database created for the study of combustion mechanisms (Konnov 2000) which seems to be appropriate for the physical conditions in the solar atmospheric plasma. Other databases like UMIST (Bennet 1988) have been used for the study of the interstellar chemistry, but they suffer from the lack of certain important reactions which can efficiently take place in the Sun’s atmosphere (e.g. the catalytic three-body formation of molecular hydrogen: H+H+He$``$H2+He). We have also investigated the possible influence of CO photodissociations, which are one-body reactions. They are mainly produced by discrete photon absorptions at wavelengths between the Lyman cutoff (912 Å) and the dissociation threshold (1120 Å). By using the photodissociation rates of Mamon, Glassgold & Huggins (1988), we have verified that the contribution of photodissociation processes to the total CO concentration is negligible for the radiation field in the solar atmosphere.

## 3 Results and discussion

We used two time-series of snapshots from the above-mentioned radiation hydrodynamical simulations, each one lasting about 3600 seconds and showing the upward propagation of acoustic wave trains growing in amplitude with height until eventually producing shocks. The first one corresponds to a relatively strong photospheric disturbance showing well-developed cool phases and pronounced hot zones at chromospheric heights (see Carlsson & Stein 1997). The strongest dynamic cycle of this simulation produces a peak-to-peak line-core brightness temperature of 1000 K, concerning the strong 3-2 R14 CO-line synthesized at disk center with the ICE approximation. However, the brightness temperature variations in most of the cycles of this simulation are 400 K approximately, which is similar to the observed values found by Uitenbroek, Rabin & Noyes (1994) under excellent seeing conditions, but larger than those inferred from the temperature histograms in Ayres & Rabin (1996). We shall refer to this as the strongly dynamic case.111This is similar to that used by Uitenbroek (2000a) and Ayres (2002) in their ICE modeling, although they considered only a segment of 190 s corresponding to the most dynamic cycle of the full simulation. The second simulation corresponds to a much less intense photospheric disturbance (see Carlsson & Stein 2002), which for strong CO lines synthesized at disk center produces a peak-to-peak line-core brightness temperature fluctuation which is always smaller than 400 K. We shall refer to this as the weakly dynamic case.

Starting from the molecular concentration given by the ICE approximation, we have followed the chemical evolution in order to obtain the temporal variation of the CO number density ($`N_{\mathrm{CO}}`$) at each height in the simulated solar atmosphere. As expected, the CO concentration is anticorrelated with the local temperature variations, yielding relatively low $`N_{\mathrm{CO}}`$ values in the hot phases and relatively high $`N_{\mathrm{CO}}`$ values in the cool phases. However, the amplitude of the local $`N_{\mathrm{CO}}`$ fluctuation is smaller than that given by the ICE approximation. In relatively low density regions characteristic of the outer atmospheric layers (e.g. around $`h=1000`$ km) the ICE approximation does a fairly good job during the hot phases, but it overestimates the CO number density during the cool phases. In contrast, in relatively high density regions characteristic of photospheric layers (e.g. around $`h=400`$ km) the ICE approximation underestimates the $`N_{\mathrm{CO}}`$ values during the hot phases. Thus, the CO abundance does not react instantaneously to the changes in the temperature because of the finite reaction rates. If one carries out a linear analysis by introducing a small temperature perturbation in a medium of given density and temperature and calculates the relaxation time needed to recover the original CO concentration one finds that, for a given density, the relaxation time is the larger the cooler the medium where the temperature perturbation is introduced. Similarly, for a given temperature of the unperturbed medium, the relaxation time increases rapidly with decreasing density. Short relaxation times are typical of high-temperature and high-density media (e.g. $`t_{\mathrm{relax}}\mathrm{\hspace{0.17em}0.006}`$ s for $`n_\mathrm{H}=10^{16}\mathrm{cm}^3`$ and $`T=6000`$ K), while long relaxation times are characteristic of low-temperature and low-density situations (e.g. $`t_{\mathrm{relax}}\mathrm{\hspace{0.17em}400}`$ s for $`n_\mathrm{H}=10^{14}\mathrm{cm}^3`$ and $`T=4000`$ K). One would find a broad range of relaxation times at a fixed height (e.g. $`h=700`$ km) in the dynamic atmosphere that results from the above-mentioned hydrodynamical simulations, simply because the existing rarefactions, compressions and temperature fluctuations are continually changing the atmospheric conditions. Obviously, the situation is highly non-linear and the relaxation time concept, although useful, loses its meaning in the real Sun. Therefore, any firm conclusion needs to be achieved via detailed numerical simulations.

Fig. 1 refers to the strongly dynamic case. The left panel shows the temporal variation in the brightness temperatures of the line-core emergent intensities at $`\mu =1`$ and $`\mu =0.1`$ in the strong 3-2 R14 CO-line. The right panel gives the ensuing fluctuations of the atmospheric height where the line-core optical depth is unity, which we use as an indicator of the ‘representative height’ where the CO line-core radiation originates. As seen in the figure, the ICE approximation does a sufficiently good job concerning the synthesis of the emergent CO spectrum at the solar disk center ($`\mu =1`$), but it largely underestimates the line-core emergent intensities at $`\mu =0.1`$ during the cool phases of the hydrodynamical simulations, producing brightness temperatures that are typically 500 K lower than those computed with the non-equilibrium CO concentration. Clearly, this is because during the cool phases the ICE approximation overestimates the ‘heights of line formation’ by about 300 km, concerning the synthesis of strong CO lines at $`\mu =0.1`$ in the strongly dynamic case. This happens because during such cool episodes the CO number density in the outer atmospheric regions is smaller than what is stipulated by the ICE approximation, resulting in decreased CO opacity in the solar chromosphere. Interestingly, in the weakly dynamic case which has smaller kinetic temperature fluctuations (but still larger than the fluctuations of observed brightness temperatures!) the ICE approximation does a much better job even at $`\mu =0.1`$.

The reader may find surprising our conclusion that the ICE approximation is suitable for modeling the CO spectrum at the solar disk center, given that Uitenbroek (2000b) found that the spatially averaged line cores of weak CO lines synthesized in the three-dimensional (3D) hydrodynamical model of Stein & Nordlund (1989) are overly dark compared to the observed ATMOS spectrum described by Farmer & Norton (1989). As we shall show in a forthcoming publication (Asensio Ramos & Trujillo Bueno; in preparation), this is because such a 3D hydrodynamical model of the solar photosphere is too cold in the CO line forming region. In fact, our ICE synthesis of the 7-6 R68 CO line in the improved 3D hydrodynamical model of Asplund et al. (2000) shows a notable agreement with the observed ATMOS spectrum, which constitutes an additional indication of the realism of the most recent 3D hydrodynamical simulations of solar surface convection (see Shchukina & Trujillo Bueno 2001 concerning the iron spectrum in such a 3D hydrodynamical model).

Fig. 2 contrasts the time-averaged CO concentration obtained assuming ICE at each time step of the strongly dynamic simulation case with that resulting from the chemical evolution. Note that the ICE approximation leads to a significant overabundance of CO in the outer layers of the atmosphere (i.e. above 700 km). Thus, the CO number density can be relatively low in such outer layers, in spite of the fact that the temporally averaged temperature of the Carlsson & Stein (1997) simulations decreases outwards and does not show any chromospheric temperature rise. As expected from the previously reported results, the ICE approximation does a sufficiently good job concerning the modeling of the temporally-averaged CO spectrum at the solar disk center ($`\mu =1`$). In contrast, the emergent CO spectrum computed close to the edge of the solar disk (i.e. at $`\mu =0.1`$) shows significantly stronger CO-lines when the ICE approximation is used, especially concerning strong CO lines like the 3-2 R14 one, for which the line-core brightness temperature is about 100 K lower than that obtained using the non-equilibrium CO number densities.

Finally, we turn our attention to the modeling issue of the off-limb CO emissions, which we have carried out by solving the radiative transfer equation in spherical coordinates at each time step of the hydrodynamical simulation222This is a suitable strategy for contrasting ICE and non-ICE results, as done also in Fig. 1. It is however clear that a truly realistic modeling of the available observations should take into account the geometrical foreshortening effect resulting from the fact that the solar chromosphere is a highly inhomogeneous medium in which the raypath goes through many different structures.. The off-limb distances where (time-averaged) CO emission has been observed depend on the line (see Ayres 2002, for a summary of the available off-limb observational results): the off-limb emission extension of strong lines like the 3-2 R14 line lies between 0.55” and 0.7” above the 4.7 $`\mu `$m continuum limb, while weaker lines like the 7-6 R68 line extend to $``$0.4”. Fig. 3 refers to the strongly dynamic case. It shows that the atmospheric region where we can have a significant off-limb emission is extremely large when the ICE approximation is used, i.e. much larger than when the spectral synthesis is carried out using the non-equilibrium CO number densities. The dashed and solid lines show the corresponding height variation of the temporally-averaged off-limb emission at the core of the strong 3-2 R14 CO-line. They show that the non-equilibrium CO chemistry improves the agreement with the available off-limb observations. First, as seen in Fig. 3, the ICE approximation predicts that appreciable emission in the 3-2 R14 line should continue to relatively large off-limb distances, while partial eclipse measurements by Clark et al (1995) show a rapid disappearance of CO emission at heights above 700 km, which is more in line with our chemical evolution calculations. Second, the representative off-limb emission extension where the normalized intensity falls to 50% of the on-disk value is larger when the non-equilibrium CO concentration is used, which helps to improve the agreement with the observations of translimb emission extensions, although they still seem to be below the most recently observed values (see Ayres 2002).

## 4 Concluding comment

Our results indicate that the CO line radiation observed close to the edge of the solar disk comes from atmospheric heights not greater than 700 km, approximately, and that the gas in these regions of the low chromosphere must be much cooler than indicated by the cool phases of the Carlsson & Stein (1997) hydrodynamical simulations. Lower temperatures will probably increase the relaxation times needed to reach the molecular equilibrium concentrations. In a forthcoming publication we will show what happens when the hydrodynamical simulations themselves are carried out taking the CO cooling into account in a way consistent with the non-equilibrium evolution of the molecular number densities.

## 5 Acknowledgements

We thank Gene Avrett, Tom Ayres and Han Uitenbroek for useful remarks. This work has been supported by the Spanish Ministerio de Ciencia y Tecnología through project AYA2001-1649, by the Norwegian Research Council through grant 146467/420 and by the European Commission via the Solar Magnetism Network.