According to thermodynamics every adiabatic system and (with no external energy added) will reach thermodynamic equilibrium or an ergodic state (2 law of thermodynamics entropy can not decrease in a closed system). State at which Temperature (or mean kinetic energy of the particles) will be same everywhere. So the answer to the question should be simple and straightforward $dT/dh = 0$.

But is it the case?

Let's open https://en.wikipedia.org/wiki/Lapse_rate page on Wikipedia and will see that adiabatic system will always create and support gradient $$dT/dh < 0$ and equal to some constant.

One can argue that this is because this system do not reach equilibrium and do not relax much. But look at calculation, if some additional relaxation took place then the gradient must decay and reach 0 at some time, but it is not the case.

One can argue and say that ok gradient is there but there is no contradiction with second law. But a thought experiment conducted by Maxwell with two columns of different gases. He showed that if the calculations are correct, then gradients of different magnitude will be generated in these columns. And therefore, in a system where these two gas columns are isolated everywhere except at the top, there will certainly be a flow of heat from a colder body to a hotter one.

Also, because gradient is exists Boltzmann distribution are mistaken.

It is simply striking how two such contradictory statements, namely the universality of the second law of thermodynamics and the gradient in the gas column in the gravitational field, are able to coexist, this is pure schizophrenia.

I also made a simple computational model which clearly shows that adibatic lapse rate you can find it on https://github.com/MaratZakirov/playground/blob/master/ideal_gas.py or in the answer of this question.

Here I list some findings which I made while discussed this question and made my model:

  1. If you consider collisions of perfect gas particles it always leads to just velocity exchange (Newton's cradle as analogy), this statement can be easily proved mathematically, because masses are the same and collision are rigid and radius of the particle is negligible. This is the real reason why you should not take into account perfect gas collisions because it does not introduce into the model any new properties.

  2. Despite the fact that Boltzmann and others derived their distributions for an ideal gas, implying the property of ergodicity of the system, in reality there is no energy mixing for the ideal gas model and particle collisions will not help at all here (see the previous paragraph). In reality, a certain entity is needed that would mix the energies of the particles and I introduced such an entity, and just after that the gradient manifested itself in all its glory.


3 Answers 3


Apparently not clear, the key point is that a system at equilibrium (with or without an external field) must have the same temperature everywhere. Missing that, there would be a net energy flux between hotter and colder parts of the system, violating the assumption of thermal equilibrium.

The above statement is a basic thermodynamics fact and can be easily derived by the entropy's maximum principle. Therefore, it is a consequence of the second principle of thermodynamics.

The atmosphere's temperature profile cannot be used as a counterexample: the atmosphere is not a system at equilibrium.

What about numerical simulation?

It is not surprising that a perfect gas does not get thermal equilibrium. Perfect gas has no mechanism to equilibrate. It is a non-ergodic system, and it is useless fr numerical simulations of thermodynamic systems. Some interaction between particles must be present to have a real thermodynamic system. The perfect gas should be taken as a limiting behavior of really interacting systems.

To clarify the previous comments, let me summarize a few facts about the interplay between ergodicity and thermodynamic behavior. Notice that I am trying to convey the main physical ideas more than trying to get the same concepts' best mathematical formulation.

A key property of any thermodynamic system is its ability to relax towards equilibrium if isolated and not initially at equilibrium. Such behavior is ensured if the system's dynamics is enough disordered to ensure that all the relevant time correlation functions among observable quantities decay to zero within the experimental time of observation. Saying in another way, a thermodynamic system loses memory of its initial state. Formally, such property of the dynamics is called mixing. If a dynamic system is mixing is also ergodic. Ergodicity is a weaker condition than mixing. It can be stated as the property that for almost every initial condition, the system's trajectory in the positions/velocities space (the phase space) visits all parts of the phase space that the system moves in. An important result of the dynamic system theory is that mixing dynamics is also ergodic. Inversely, a non-ergodic system cannot be mixing.

That the ideal gas is not ergodic may be clear by thinking a simple initial condition: a cubic box, half of the particles are at rest, and half have the same speed. Part of the available phase space will never get visited by such a system. Moreover, the particles' subsystem at rest has zero temperature and the remaining a finite temperature. Clearly, this is neither an ergodic system nor a system at thermodynamic equilibrium.

To get a mixing system, it is enough to add an even small interaction among particles or with the walls to introduce a dynamics chaotic enough to recover the mixing property. In a mixing system, one can start with whatever velocity distribution, and waiting enough, it is possible to get a well-equilibrated system in an interacting system.

I also notice that neither Maxwell-Boltzmann nor the uniform distribution is the correct velocity distribution at equilibrium in an isolated system. Even if one starts with one distribution, the velocity distribution will evolve towards the correct equilibrium values after some relaxation time, depending on the thermodynamic state. Monitoring the time evolution of the velocity distribution should be enough to show the phenomenon, at least when starting with a uniform distribution. Since the microcanonical velocity distribution and Maxwell-Boltzmann are very close for a system of a few thousand particles, I do not think it would be easy to notice the difference. However, a careful measurement of temperature at different heights should be enough. Moreover, it is also important for this kind of study to estimate the statistical error on results before any quantitative conclusion could be drawn.

  • $\begingroup$ Warning! You are not allowed to use unproven (in this model) things to prove others also unproven. Please firstly prove that the gradient is zero and then lets talk about the entropy that will be correct. You have code on screen and may simulate us much time as you like I think 3700 ticks is enough. As I said I also tried uniform distribution but it showed $dT/dh < 0$ please check it by yourself. $\endgroup$ Nov 23, 2020 at 12:51
  • $\begingroup$ @MaralZakirov I had to re-edit my previous answer after realizing the big physical mistake you have in your code. No way to get thermodynamics by simulating a perfect gas. $\endgroup$ Nov 23, 2020 at 12:54
  • $\begingroup$ please tell me more. $\endgroup$ Nov 23, 2020 at 12:56
  • $\begingroup$ @MaratZakirov I can try to elaborate adding something more to my answer. Unfortunately, I cannot do it in the next two hours. I'll turn to my answer later, $\endgroup$ Nov 23, 2020 at 12:58
  • $\begingroup$ @GiorgioP is referringto the fact that you have no collision mechanism in your gas, and therefore no wayto relax to equilibrium. If you put in an initial velocity distribution and bounce elastically off the walls, the distributon will stay unchanged. $\endgroup$
    – mike stone
    Nov 23, 2020 at 13:36

the OP does not say what he means by "his equation" but I assuming that the the OP's question is about Boltzmann's law $$ \rho(h)\propto e^{-mgh/kT} $$ for the density profile of an isothermal atmosphere, and not one about thermal equilibrium. This simple atmospheric density law assumes that the atmosphere is isothermal

There is no reason for the distribution in a real atmospheric column to be isothermal. Indeed, in the lower part of the Earth's atmosphere where it is stirred by convection, the temperature falls off with height at roughly the adiabatic lapse rate. This is because if a parcel of air moves up into a lower pressure region it expands and so cools. Similarly a parcel moving downwards is compressed and gets hotter.

Of course a non-uniform temperature is not in thermal equilibrium, only in mechanical equilibrium. For thermal equilibrium one does not assume that the temperature is constant one can prove it in suitable statistical mechanical settings.

  • $\begingroup$ The article that you showed me on the wiki claims exactly what I got on the simulation, the gradient exists, and it does not go away, and therefore Boltzmann's calculations and the omnipotence of the second law of thermodynamics are questioned, to put it mildly $\endgroup$ Nov 25, 2020 at 12:43

Taken from Wikipedia page Lapse_rate


A formal definition from the Glossary of Meteorology [Todd S. Glickman (June 2000). Glossary of Meteorology (2nd ed.). American Meteorological Society, Boston. ISBN 978-1-878220-34-9. (Glossary of Meteorology)] is:

The decrease of an atmospheric variable with height, the variable being temperature unless otherwise specified.

Typically, the lapse rate is the negative of the rate of temperature change with altitude change:

$Г = -\frac{dT}{dz}$

where $Γ$ (sometimes $L$) is the lapse rate given in units of temperature divided by units of altitude, $T$ is temperature, and $z$ is altitude.

Mathematics of the adiabatic lapse rate

These calculations use a very simple model of an atmosphere, either dry or moist, within a still vertical column at equilibrium (my comment: equilibrium Karl!!!).

Dry adiabatic lapse rate

Thermodynamics defines an adiabatic process as:

$PdV = \frac{-VdP}{\gamma}$

the first law of thermodynamics can be written as

$mc_vdT - \frac{VdP}{\gamma} = 0$ (my comment: adiabatic process no external energy flow)

Also, since $α = V / m$ and $γ = c_p/c_v$, we can show that:

$c_pdT - \alpha dP=0$

where $c_p$ is the specific heat at constant pressure and $α$ is the specific volume.

Assuming an atmosphere in hydrostatic equilibrium:

$dP=-\rho gdz$

where $g$ is the standard gravity and $\rho$ is the density. Combining these two equations to eliminate the pressure, one arrives at the result for the dry adiabatic lapse rate (DALR)

lapse rate (DALR),[11]

$Γ_d = -\frac{dT}{dz} = \frac{g}{c_p} = 9.8 ∘C/km$

My comment:

So the answer for my question is where is no reason $dT/dh$ to be equal to 0 in the gravitational field. And a simulation I made on python shows exactly that higher layers of gas columns are cooler and stays cooler regardless of time ticks I spent to simulate. If we recall the thought experiment that Maxwell and later Tsiolkovsky did with two gas columns, then we get a direct violation of the principle of universal applicability of the second law of thermodynamics (because different gases have different lapse rate), recognized among many physicists. Let the one who read this now live with this insane contradiction.


Updated experiment with:

  1. Maxwell distribution of velocities (normal distribution of velocity projections).
  2. Random turns of velocity. It does not help at all for ergodicity of the system.
  3. Introducing energy mixer at the bottom which enables the system to become ergodic.


# Ideal GAS model
#%matplotlib inline
#from IPython.display import HTML
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation


# Macro parameters
x0 = 0; x1 = 16
y0 = 0; y1 = 16
# TODO set to 8000 for evaluation
N = 8400
E = 20000
L = 10
period = 30
ENANLEMIX = 'maxwell'


# State
X = np.random.uniform((x0, y0), (x0 + 1, y0 + 1), (N, 2))
# TODO understand velocity distribution
    V = 13 * np.random.randn(N, 2)
    V = 30 * (np.random.rand(N, 2) - 0.5)

g = np.array([0, -9.8])
dt = 0.01

lines = np.array([[x0, y0, x1, y0],
                  [x0, y1, x1, y1],
                  [x0, y0, x0, y1],
                  [x1, y0, x1, y1]])

# For statistics storing
x_data = []
v_data = []
Ek_data = []
Ep_data = []
Hkh_data = []
Hk2h_data = []
Hph_data = []
Hfh_data = []

# Every particle in ideal gas model has same amount of energy
Ef_const = X[:, 1] * (-g[1]) + (V * V / 2).sum(1)
dE = (Ef_const.mean() + 4 * Ef_const.std()) / L
dE2 = (np.sqrt(Ef_const).mean() + 4 * np.sqrt(Ef_const).std()) / L
dh = (y1 - y0) / L

for epoch in range(E):
    X_n = X + V * dt + 0.5 * g * dt ** 2
    V = V + g * dt

    # Fix velocities to return particles in volume
    Ix = (X_n[:, 0] < x0) + (X_n[:, 0] > x1)
    Iy = (X_n[:, 1] < y0) + (X_n[:, 1] > y1)
    Iy0 = X_n[:, 1] < y0
    Iy1 = X_n[:, 1] > y1

    V[Ix, 0] *= -1
    V[Iy0, 1] = np.abs(V[Iy0, 1])
    V[Iy1, 1] = -np.abs(V[Iy1, 1])

    if ENANLEMIX == 'uniform' and Iy0.sum() >= 2:
        alpha = 0.4
        n = Iy0.sum()

        # Storing kinetic energy of affected particles
        E = np.linalg.norm(V[Iy0], axis=1, keepdims=True) ** 2

        # Creating random mixing matrix
        M = np.random.uniform(low=0, high=1, size=(n, n))
        M = M / M.sum(axis=1, keepdims=True)
        M = np.eye(n) * (1 - alpha) + alpha * M

        E_n = np.dot(M, E)
        E_n = E_n * (E.sum() / E_n.sum())

        V[Iy0] = (V[Iy0] / np.sqrt(E)) * np.sqrt(E_n)

    if ENANLEMIX == 'maxwell' and Iy0.sum() >= 2:
        alpha = 0.4
        n = Iy0.sum()

        # Storing kinetic energy of affected particles
        E = np.linalg.norm(V[Iy0], axis=1, keepdims=True) ** 2

        E_n = np.clip(E.std() * np.random.randn(n, 1) + E.mean(), a_min=0, a_max=1000000)
        E_n = E_n * (E.sum() / E_n.sum())
        E_n = E * (1 - alpha) + E_n * alpha
        E_n = E_n * (E.sum() / E_n.sum())

        V[Iy0] = (V[Iy0] / np.sqrt(E)) * np.sqrt(E_n)

        eps = 1.0
        # Find situation when particles are just in inner volume
        I_n_pos_in = (X_n[:, 0] > (x0 + eps)) * (X_n[:, 0] < (x1 - eps)) * (X_n[:, 1] > (y0 + eps)) * (X_n[:, 1] < (y1 - eps))
        I_o_pos_not = ~((X[:, 0] > (x0 + eps)) * (X[:, 0] < (x1 - eps)) * (X[:, 1] > (y0 + eps)) * (X[:, 1] < (y1 - eps)))
        I_rr = I_n_pos_in * I_o_pos_not

        ralpha = np.random.uniform(low=-0.1, high=0.1, size=(I_rr).sum())
        Rrm = np.array([[np.cos(ralpha), -np.sin(ralpha)],
                        [np.sin(ralpha), np.cos(ralpha)]]).transpose([2, 0, 1])

        # Calculate potential new velocities
        V[I_rr] = np.stack([(V[I_rr] * Rrm[..., 0]).sum(1),
                            (V[I_rr] * Rrm[..., 1]).sum(1)], axis=1)

        if (V != V).any():
            assert 0

    X = X_n

    if epoch % period == 0:
        if (V != V).any():
            assert 0

        # Store statistics
        Ek_data.append((V * V / 2).sum(1))
        Ep_data.append(X[:, 1] * (-g[1]))

        # Calculating entropy
        Ek_lev = np.clip(Ek_data[-1] / dE, 0, L - 1).astype(int)
        Ek2_lev = np.clip(np.sqrt(Ek_data[-1]) / dE2, 0, L - 1).astype(int)
        Ep_lev = np.clip(Ep_data[-1] / dE, 0, L - 1).astype(int)
        Ef_lev = np.clip((Ep_data[-1] + Ek_data[-1]) / dE, 0, L - 1).astype(int)
        h_lev = np.clip((X[:, 1] - y0) / dh, 0, L - 1).astype(int)
        p_kh = np.bincount(h_lev * L + Ek_lev, minlength = L * L) / N
        p_k2h = np.bincount(h_lev * L + Ek2_lev, minlength = L * L) / N
        p_ph = np.bincount(h_lev * L + Ep_lev, minlength=L * L) / N
        p_pf = np.bincount(h_lev * L + Ef_lev, minlength=L * L) / N
        Hkh_data.append(-np.sum(p_kh * np.log(p_kh + 0.00001)))
        Hk2h_data.append(-np.sum(p_k2h * np.log(p_k2h + 0.00001)))
        Hph_data.append(-np.sum(p_ph * np.log(p_ph + 0.00001)))
        Hfh_data.append(-np.sum(p_pf * np.log(p_pf + 0.00001)))

x_data = np.stack(x_data)
v_data = np.stack(v_data)
Ek_data = np.array(Ek_data)
Ep_data = np.array(Ep_data)
Hkh_data = np.array(Hkh_data)
Hk2h_data = np.array(Hk2h_data)
Hph_data = np.array(Hph_data)
Hfh_data = np.array(Hfh_data)

print('Stalled particles:', (x_data[-234:, :, 1].max(axis=0) < 0).sum())

# Print dT/dh = const * dE(Ek)/dh of stationary state
Levs = 14
stEp = len(Ek_data) // 6
Ek0 = Ek_data[-stEp:]
h0 = x_data[-stEp:, :, 1]

T_whole = []
N_whole = []
for e in range(stEp):
    N_cur = np.zeros(Levs)
    T_cur = np.zeros(Levs)

    Ek0_cur = Ek0[e]
    h0_cur = (Levs * (h0[e] + 0.2) / (y1 - y0 + 1)).astype(int)

    for l in range(Levs):
        if (h0_cur == l).sum() > 0:
            T_cur[l] = Ek0_cur[h0_cur == l].mean()
            N_cur[l] = (h0_cur == l).sum()


N_whole = np.stack(N_whole).mean(0)
T_whole = np.stack(T_whole).mean(0)

print('T_whole', T_whole)
print('N_whole', N_whole, N_whole.sum())

fig, ax1 = plt.subplots()

ax1.plot(N_whole, color='tab:blue')
ax1.set_xlabel('h level')
ax1.set_ylabel('Num of particles')
ax2 = ax1.twinx()
ax2.plot(T_whole, color='tab:red')
ax2.set_ylabel('T or average kinetic energy of particles')


# Print correlation
Ek = Ek_data
Ep = Ep_data
Ef = Ek_data + Ep_data
h = x_data[..., 1]
rk = ((Ek - Ek.mean(1, keepdims=True)) * (h - h.mean(1, keepdims=True))).mean(1) / Ek.std(1) / h.std(1)
rp = ((Ep - Ep.mean(1, keepdims=True)) * (h - h.mean(1, keepdims=True))).mean(1) / Ep.std(1) / h.std(1)
rf = ((Ef - Ef.mean(1, keepdims=True)) * (h - h.mean(1, keepdims=True))).mean(1) / Ef.std(1) / h.std(1)
plt.legend(['Ek with h', 'Ep with h', 'Ef with h'])


plt.ylim(0, Ef_const.sum() * 2)
plt.plot(Ep_data.sum(1) + Ek_data.sum(1))
plt.legend(['Current full energy', 'Old full energy'])


plt.ylabel('H (Entropy)')
plt.legend(['H(Ek, h)', 'H(v, h)', 'H(Ep, h)', 'H(Ef, h)'])


fig, _ = plt.subplots()
plt.xlim(x0-1, x1+1)
plt.ylim(y0-1, y1+1)

if not MAKEANIM:
    if not SAVEFIG:

scatter = plt.scatter(x_data[0, :, 0], x_data[0, :, 1], s=4)
def animate(i):
    scatter.set_array(Ek_data[i] + Ep_data[i])
    return scatter,

a = FuncAnimation(fig, animate, frames=len(x_data), interval=20, blit=True, repeat=True)

Graphs for checking:

Temperature gradient red line which shows that gradient is constant which is totally agreeds with the lapse rate theory.

enter image description here

Full, Potential, Kinetic energy correlation with height. From this picture you may see clear correlation between kinetic energy (~temperature of gas particle) and height

enter image description here

Full energy of the system does not change

enter image description here

  • $\begingroup$ Your code is highly suspect to me. I don't see any indication of an appropriate collision simulation here, which is required to bring this system to equilibrium. Random 3D collisions between hard spheres are more complex than you seem to think, and involve a random transfer of energy from the faster particle to the slower depending on the collision's impact parameter. To test this, try to simulate an ideal gas in a box. If you start all of the particles with the same energy, then after a sufficient number of collisions the velocity distribution will be Maxwellian. If your simulation [...] $\endgroup$
    – J. Murray
    Dec 9, 2020 at 14:45
  • $\begingroup$ [...] does not have this feature, then your simulation is wrong. Numerically simulating collisions is non-trivial. More generally, you are now in the realm of computational non-equilibrium statistical mechanics, which is very deep water. Regarding your statement "Let the one who read this now live with this insane contradiction," I would suggest you approach this very complex study with a bit more humility. $\endgroup$
    – J. Murray
    Dec 9, 2020 at 14:49
  • $\begingroup$ I do not have to simulate collisions due to properties of perfect gas and I wrote about this. "If you consider collisions of perfect gas particles it always leads to just velocity exchange (Newton's cradle as analogy), this statement can be easily proved mathematically, because masses are the same and collision are rigid and radius of the particle is negligible. This is the real reason why you should not take into account perfect gas collisions because it does not introduce into the model any new properties." $\endgroup$ Dec 9, 2020 at 15:02
  • 1
    $\begingroup$ That is only true if the particles collide with zero impact parameter. Real collisions occur at glancing angles, and so the energy transfer could be extremely small (if they just brush past each other), or the "velocity exchange" that you refer to, or anything in between. $\endgroup$
    – J. Murray
    Dec 9, 2020 at 15:04
  • $\begingroup$ See e.g. this paper on an algorithm for simulating collisions. As you may notice, it is not trivial. Collisions are not easy to simulate. $\endgroup$
    – J. Murray
    Dec 9, 2020 at 15:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.