Distribution de Boltzmann

1. Introduction

Ce document est une introduction aux méthodes de Monte-Carlo appliquées à la physique statistique. Un exemple particulier est considéré : un système formé de N particules indépendantes à énergie quantifiée, dont l'énergie totale est constante. Après avoir introduit la distribution de Boltzmann pour les états des particules, on fera une simulation de Monte-Carlo pour mettre en évidence cette distribution.

2. Définition du système

On considère un système formé de N particules indépendantes et discernables. Pour chaque particule, il y a M états possibles, numérotés de 0 à M-1. À chaque état m est associée une énergie ε_m. Ces énergies sont régulièrement espacées. On note ε l'espacement entre deux niveaux voisins. Le nombre M d'états accessibles peut éventuellement être infini.

Un exemple physique de système de ce type est un cristal, dont chaque atome est modélisé par un oscillateur harmonique. On démontre en mécanique quantique que l'énergie de l'état m d'un oscillateur harmonique est :

On suppose ici que tous ces oscillateurs sont indépendants, ce qui n'est pas exact dans un cristal réel.

La figure suivante montre les niveaux d'énergie de quelques particules du système. Il est important de bien noter qu'à chaque niveau d'énergie correspond un et un seul état de la particule. On dit que ces niveaux d'énergie sont non dégénérés.

Les particules sont indépendantes dans le sens où l'énergie d'interaction entre les particules est négligeable devant l'énergie des particules. Il y a tout de même des petites interactions entre particules qui permettent au système d'atteindre l'équilibre thermodynamique.

Ce système à énergie discrète a été étudié par Boltzmann bien avant la naissance de la mécanique quantique. Pour Boltzmann, la discrétisation de l'énergie était un moyen de rendre les états dénombrables. Cette démarche a été reprise par Planck pour l'étude du rayonnement du corps noir.

3. Micro-états

Un micro-état du système est l'état de tous les degrés de liberté du système. Dans le cas présent, un micro-état est défini par l'état de chacune de ses particules. Il est donc défini par N nombres entiers compris entre 0 et M-1 :

On suppose que le système est isolé. Son énergie totale est alors une constante. Celle-ci s'écrit :

où e_i est l'énergie de la particule i.

L'énergie totale est nécessairement multiple de ε. On pose donc :

ce qui donne :

On recherche le nombre Ω de micro-états pour cette énergie totale. Il est égal au nombre de manières de répartir l'énergie totale (le nombre Q) sur les N particules. La figure suivante permet de comprendre le dénombrement à faire. Les cercles représentent les Q éléments d'énergie ε. Les croix sont les frontières qui permettent de répartir l'énergie sur les différentes particules; il y en a N-1.

Il s'agit de répartir Q+N-1 symboles (croix et cercles) sur autant d'emplacement. Pour le premier symbole, il y a Q+N-1 possibilités, pour le deuxième Q+N-2, pour le troisième Q+N-3, etc. Le nombre total de possibilités est donc le produit factoriel (Q+N-1)!. L'échange de deux cercles donne deux micro-états identiques. Il faut donc diviser ce nombre par Q!. L'échange de deux croix donne aussi deux micro-états identiques. Il faut donc diviser par (N-1)!. On obtient finalement le nombre de micro-états :

On cherche à évaluer ce nombre pour un système comportant un très grand nombre de particules, avec un très grand nombre d'éléments d'énergie à répartir. Pour cela, on calcule son logarithme néperien et on utilise l'approximation de Stirling :

En négligeant 1 devant N, on obtient :

En introduisant l'énergie moyenne (divisée par ε) par particule q=Q/N, on obtient :

import numpy
from matplotlib.pyplot import *

def omega(q):
    return (1+q)*numpy.log(1+q)-q*numpy.log(q)
    
q = numpy.linspace(1e-2,100,1000)
figure()
w = omega(q)
plot(q,w)
xlabel("q")
ylabel(r"$\ln(\Omega)/N$")
grid()
                  

Si par exemple l'énergie moyenne par particule est 20 (l'énergie totale est alors 20Nε), alors le nombre de micro-états est 4^N. Lorsque N est grand, ce nombre est extrêmement grand. Même pour N=100, il dépasse de très loin les capacités de calcul des ordinateurs. Il est donc impossible d'explorer systématiquement tous ces états.

On peut donner une estimation plus simple du nombre d'état lorsque la condition suivante est vérifiée :

On suppose de plus que Q est défini à plus ou moins une unité. La figure suivante montre une interprétation géométrique pour 2 particules :

Le quadrillage a un pas d'une unité. Les points vérifiant

sont dans le domaine colorié. Si les conditions ci-dessus sont vérifiées, on peut considérer que le nombre de points dans ce domaine est proportionnel à Q. En dimension 3, il est proportionnel à 2. En dimension N, il est proportionnel à N-1. Comme N est grand devant 1, on obtient la loi d'échelle suivante pour le nombre d'états :

4. Probabilité des micro-états

Pour une énergie E donnée, il y a un très grand nombre de micro-états correspondants, d'autant plus grand que N est grand. Pour un système thermodynamique, le nombre de micro-états correspondant à une énergie donnée est extrêmement grand.

La physique statistique repose sur deux hypothèses fondamentales. La première consiste à considérer qu'un système à l'équilibre thermodynamique, en l'occurence un système ayant une énergie E bien définie, parcourt tous les micro-états qui lui sont accessibles, en un temps plus ou moins long. Un système vérifiant cette propriété est dit ergodique. Lorsqu'un système est ergodique, on peut s'intéresser à ses micro-états d'un point de vue probabiliste seulement, sans chercher à savoir exactement comment le système passe d'un micro-état à un autre. La seconde hypothèse est la suivante :

Pour un système isolé d'énergie E, la probabilité d'un micro-état s'écrit donc :

où Ω(N,E) est le nombre de micro-états qui ont l'énergie E. Cette hypothèse fondamentale a été introduite par Boltzmann, qui a aussi défini l'entropie :

Nous avons calculé ci-dessus l'entropie pour le système de particules indépendantes.

Le choix de la fonction logarithme pour définir l'entropie permet d'assurer son extensivité. Considérons une division en deux parties égales du système. Le nombre de micro-états de l'ensemble est le produit des nombres de micro-états des deux parties :

On obtient ainsi :

La relation $\mathrm{(12)}$ nous donne l'estimation suivante :

Pour un système thermodynamique, il s'agit donc d'une fonction très rapidement croissante de l'énergie. L'entropie est alors :

La température absolue est définie par :

Cette estimation conduit donc à une énergie proportionnelle à la température.

Cette expression est une estimation. La forme exacte de l'énergie en fonction de la température dépend du système considéré. Néanmoins, cette forme de l'énergie proportionnelle à la température est classique en thermodynamique puisqu'elle correspond à une capacité thermique constante.

5. Distribution de Boltzmann

On s'intéresse à présent à une particule, dans le système dont l'énergie totale est constante. Il est bien évident que l'énergie d'une particule n'est pas du tout constante, car elle échange de l'énergie avec ses voisines. Toutes les énergies possibles inférieures à Mε sont possibles pour cette particule. Rappelons qu'à chaque énergie d'une particule i correspond un et un seul état de cette particule, donné par le nombre m_i. On s'intéresse à la probabilité de cet état de la particule, qui est certainement une fonction de e. Cette probabilité est proportionnelle au nombre de micro-états du système pour lesquels cette particule a l'énergie e, c'est-à-dire le nombre de micro-états d'un système de N-1 particules dont l'énergie totale est E-e :

où C est une constante. Comme l'énergie e est extrêmement petite devant l'énergie totale du système E, on peut faire un développement limité. On remarque tout d'abord que si N est très grand, on peut confondre N-1 et N. On commence par faire une estimation en utilisant la relation $\mathrm{(17)}$ :

Pour effectuer le développement limité, il faut tenir compte du fait que l'énergie E est de l'ordre de NkT et que N est très grand (il tend vers l'infini à la limite thermodynamique). On écrit donc :

La limite de cette expression lorsque N tend vers l'infini est :

Cette estimation suggère de faire un développement limité du logarithme de la probabilité à l'ordre 1 :

En ne gardant que la partie dépendant de e et en introduisant l'entropie définie plus haut, on obtient :

On a donc :

La constante C est déterminée en écrivant que la somme des probabilités pour tous les états de la particule doit être égale à 1 (condition de normalisation) :

On obtient finalement la probabilité qu'une particule soit dans l'état m :

Cette distribution est la distribution de Boltzmann. La démonstration ci-dessus montre qu'elle se généralise à un système en contact avec un thermostat de température T, qui peut être constitué d'une ou plusieurs particules. La probabilité d'un micro-état μ de ce système (qui n'est pas un système isolé), s'écrit :

La somme écrite au dénominateur est une somme sur tous les micro-états du système.

Cette probabilité ne dépend que de l'énergie du micro-état et de la température du thermostat. Il faut bien noter qu'il s'agit de la probabilité d'un micro-état, et non pas de la probabilité que le système ait une énergie donnée. Pour une énergie donnée, le nombre de micro-états correspondant est Ω(N,E). Dans le cas d'une seule particule, il n'y a qu'un seul état. Lorsque N est grand, ce nombre est très grand, de l'ordre de E^N.

6. Simulation de Monte-Carlo

import numpy
import numpy.random
import math
from matplotlib.pyplot import *
                

def tirages(E,N,epsilon,emax,nt):
    nh = int(emax/epsilon)
    histogramme = numpy.zeros(nh,dtype=numpy.float64)
    energies = numpy.arange(nh)*epsilon
    for k in range(nt):
        partition = numpy.random.random_sample(N-1)*E
        partition = numpy.sort(partition)
        partition = numpy.append(partition,E)
        p = 0
        for i in range(N):
            e = partition[i]-p
            p = partition[i]
            m = math.floor(e/epsilon)
            if m<nh:
                histogramme[m] += 1
    return (energies,histogramme/(nt*N))

                

6.c. Simulation

Un centaine de particules est un nombre assez grand pour mettre en évidence la distribution de Boltzmann. On trace l'histogramme pour trois valeurs différentes de l'énergie.

N=100
nt=100000
emax = 100
epsilon = 0.5

figure()
E1=200
(e1,h1) = tirages(E1,N,epsilon,emax,nt)
plot(e1,h1,'o',label='E=200')
E2=400
(e2,h2) = tirages(E2,N,epsilon,emax,nt)
plot(e2,h2,'o',label='E=400')
E3=800
(e3,h3) = tirages(E3,N,epsilon,emax,nt)
plot(e3,h3,'o',label='E=800')
xlabel("energie")
ylabel("proba")
legend(loc='upper right')
axis([0,30,0,0.3])
grid()
                

figA.pdf

L'aspect des courbes suggère une forme exponentielle. Les états les plus probables pour une particule sont les états de faible énergie et l'état le plus probable est celui d'énergie nulle. L'énergie moyenne d'une particule est E/N. Elle augmente avec l'énergie totale, alors que l'état le plus probable reste toujours celui d'énergie nulle. Ce résultat n'a rien d'intuitif. Un raisonnement rapide aurait pu nous suggérer une énergie la plus probable proche de l'énergie moyenne; ce n'est pas du tout le cas.

Pour vérifier la forme exponentielle, on trace le logarithme de la probabilité (divisée par p(0)), multipliée par E/N.

figure()
plot(e1,numpy.log(h1/h1[0])*E1/N,label='E=200')
plot(e2,numpy.log(h2/h2[0])*E2/N,label='E=400')
plot(e3,numpy.log(h3/h3[0])*E3/N,label='E=800')
xlabel("energie")
ylabel("E/N*ln(p)")
legend(loc='upper right')
axis([0,50,-60,0])
grid()
                  

figB.pdf

Cela confirme la forme exponentielle et suggère l'expression suivante :

$$p\left(e\right)=p\left(0\right)e^{-\frac{Ne}{E}}\text{(32)}$$

Par identification avec la distribution de Boltzmann, on en déduit la température :

$$kT=\frac{E}{N}\text{(33)}$$

ce qui est précisément la relation $\mathrm{(20)}$. Pour cette simulation, la température est égale à l'énergie moyenne des particules.

Voir aussi Simulation de Monte-Carlo d'un système à énergie constante.