Équation de Boltzmann-Peierls

En physique du solide l'équation de Boltzmann-Peierls décrit l'évolution de la fonction de distribution des phonons dans un solide cristallin. Elle a été établie par Rudolf Peierls en 1929. Diverses approximations de type BGK ont été proposées, la plus précise par Joseph Callaway en 1959. Dans les dernières décennies des efforts ont porté sur le développement de méthodes analogues à celles de la mécanique des fluides, permettant une grande efficacité de la résolution du problème moyennant quelques simplifications acceptables pour les applications pratiques.

Cette équation décrit le transfert d'énergie thermique au niveau microscopique, dans un domaine où l'approximation diffusive (loi de Fourier) n'est pas valide, le libre parcours moyen des phonons étant du même ordre de grandeur ou plus grand que la longueur caractéristique associée à la température $\textstyle {\frac {T}{||\nabla T||}}$. Ce type de problème est présent en microélectronique et la méthode s'applique également aux électrons dans les solides.

L'équation de Boltzmann-Peierls

Les phonons sont des bosons décrits par leur densité numérique $\textstyle \,f_{\mathbf {k} ,s}(\mathbf {x} ,t,\mathbf {k} )\,$ où $\textstyle \mathbf {x}$ est la variable d'espace, $\textstyle t$ le temps, $\textstyle \mathbf {k}$ le vecteur d'onde supposé continu et $\textstyle \,s=1,2,3\,$ désigne la polarisation. Cette grandeur est liée à la densité d'états $\textstyle \rho _{\mathbf {k} ,s}(k)$ caractéristique du matériau et à la distribution statistique de ces états $\textstyle n_{\mathbf {k} ,s}(\mathbf {x} ,t)$ par $\textstyle f_{\mathbf {k} ,s}=\rho _{\mathbf {k} ,s}n_{\mathbf {k} ,s}$.

Dans ce qui suit on s'intéressera à une polarisation particulière, les calculs étant les mêmes pour toutes : l'indice s disparaît donc.

L'équation de conservation de $\textstyle f_{\mathbf {k} }$ s'écrit :

${\frac {\mathrm {d} f_{\mathbf {k} }}{\mathrm {d} t}}={\frac {\partial f_{\mathbf {k} }}{\partial t}}+\mathbf {v} _{\mathbf {k} }\cdot \nabla f_{\mathbf {k} }={\frac {\mathrm {d} f_{\mathbf {k} }}{\mathrm {d} T}}\left({\frac {\partial T}{\partial t}}+\mathbf {v} _{\mathbf {k} }\cdot \nabla T\right)=S_{\mathbf {k} }$

où $\textstyle T$ est la température thermodynamique et $\textstyle \mathbf {v} _{\mathbf {k} }$ est la vitesse de propagation ou vitesse de groupe liée à la pulsation de l'onde $\textstyle \omega _{\mathbf {k} }$ par $\textstyle \mathbf {v} _{\mathbf {k} }={\frac {\mathrm {d} \omega _{\mathbf {k} }}{\mathrm {d} k}}{\frac {\mathbf {k} }{k}}$. La relation $\textstyle \omega _{\mathbf {k} }(k)$ constitue la relation de dispersion.

La connaissance du terme source $\textstyle S_{\mathbf {k} }$ est au centre du problème. Prenons l'exemple simple de la diffusion d'un phonon sur un centre diffuseur comme un défaut du réseau cristallin. La probabilité de transition $\textstyle \mathbf {k} \to \mathbf {k} '$ est proportionnelle à $\textstyle \Lambda _{\mathbf {k} \to \mathbf {k} '}n_{\mathbf {k} }(1+n_{\mathbf {k} }')$ où $\textstyle \Lambda _{\mathbf {k} \to \mathbf {k} '}$ est la probabilité de transition spontanée.

En prenant en compte la réversibilité $\textstyle \Lambda _{\mathbf {k} \to \mathbf {k} '}\,=\,\Lambda _{\mathbf {k} '\to \mathbf {k} }$ on obtient le terme source pour la diffusion :

$S_{\mathbf {k} }={\frac {1}{2}}\sum _{\mathbf {k} }\sum _{\mathbf {k'} }\Lambda _{\mathbf {k} \to \mathbf {k} '}\left(n_{\mathbf {k} }-n_{\mathbf {k'} }\right)$

le facteur $\textstyle {\frac {1}{2}}$ étant destiné à corriger les doubles comptes.

On peut faire de même pour les processus à 3 phonons dont les probabilités de transition sont plus complexes car faisant intervenir une approximation à l'ordre 3 au moins du potentiel interatomique cristallin (il n'y a pas d'interaction phonon-phonon à l'ordre 2).

Le calcul de ces termes est coûteux et la résolution de l'équation de Boltzmann par une méthode stochastique ou déterministe est elle-même coûteuse (calcul d'un grand nombre de fonctions de distribution en tout point et à tout instant, un terme source complexe). Des approximations ont donc été recherchées.

Première approximation

On peut écrire en première approche le terme source sous une forme utilisant la méthode de Bhatnagar-Gross-Krook comme un simple terme de rappel vers la valeur d'équilibre $\textstyle f_{\mathbf {k} }^{(0)}=\rho _{k}n_{\mathbf {k} }^{(0)}$ où $\textstyle n_{\mathbf {k} }^{(0)}$ est donné par la statistique de Bose-Einstein pour un potentiel chimique nul :

$\textstyle n_{\mathbf {k} }^{(0)}=\left(e^{\frac {\hbar \omega _{\mathbf {k} }}{k_{B}T}}-1\right)^{-1}$

Ce rappel s'effectue avec le temps caractéristique $\textstyle \tau _{\mathbf {k} }$ :

$S_{\mathbf {k} }=-{\frac {f_{\mathbf {k} }-f_{\mathbf {k} }^{(0)}}{\tau _{\mathbf {k} }}}$

Démonstration

Pour un milieu homogène où la distribution d'équilibre est $f_{\mathbf {k} }^{(0)}$ l'équation de Boltzmann-Peierls pour $\textstyle f_{\mathbf {k} }(t)$ dans l'approximation de relaxation se réduit à :

${\frac {\mathrm {d} f_{\mathbf {k} }}{\mathrm {d} t}}=\tau _{\mathbf {k} }^{-1}(f_{\mathbf {k} }^{(0)}-f_{\mathbf {k} })$

qui décrit un phénomène de relaxation vers l'équilibre avec la constante de temps $\textstyle \tau _{\mathbf {k} }$ :

$f_{\mathbf {k} }(t)=f_{\mathbf {k} }(0)+(f_{\mathbf {k} }^{(0)}-f_{\mathbf {k} }(0))e^{-{\frac {t}{\tau _{\mathbf {k} }}}}$

Si à présent dans l'équation de Boltzmann on suppose un faible écart à l'équilibre on peut confondre la dérivée de $\textstyle n_{\mathbf {k} }(\mathbf {x} ,t)$ avec celle de $\textstyle n_{\mathbf {k} }^{(0)}(T(\mathbf {x} ,t))$ :

${\frac {\mathrm {d} f_{\mathbf {k} }}{\mathrm {d} t}}\approx {\frac {\mathrm {d} f_{\mathbf {k} }^{(0)}}{\mathrm {d} t}}$

D'où la solution du problème stationnaire :

$f_{\mathbf {k} }-f_{\mathbf {k} }^{(0)}=-\mathbf {v} _{\mathbf {k} }\cdot \nabla T\,\tau _{\mathbf {k} }\,{\frac {\hbar \omega _{\mathbf {k} }}{k_{B}T^{2}}}\,{\frac {e^{\frac {\hbar \omega _{\mathbf {k} }}{k_{B}T}}}{\left(e^{\frac {\hbar \omega _{\mathbf {k} }}{k_{B}T}}-1\right)^{2}}}$

Cette approximation est bonne lorsque l'on traite des processus dits « résistifs » (car participant à la résistivité thermique) qui ne conservent pas la quantité de mouvement comme la diffusion dite processus umklapp ou pour les phénomènes de diffusion d'un phonon sur un centre diffuseur (un défaut du cristal).

Lorsque plusieurs phénomènes se superposent ce temps caractéristique est obtenu en appliquant la règle de Matthiessen qui exprime l'indépendance des phénomènes :

$\tau ^{-1}=\sum _{i}\tau _{i}^{-1}$

Deuxième approximation (Callaway)

Si, en plus des processus résistifs, il existe des processus dits « normaux » qui respectent la conservation de la quantité de mouvement une expression différente est proposée par Callaway :

${\frac {\partial f_{\mathbf {k} }}{\partial t}}+\mathbf {v} _{\mathbf {k} }\cdot \nabla f_{\mathbf {k} }=-{\frac {f_{\mathbf {k} }-f_{\mathbf {k} }^{(0)}}{\tau _{\mathbf {k} }}}-{\frac {f_{\mathbf {k} }-f_{\mathbf {k} }^{(1)}}{\tau _{\mathbf {k} }^{*}}}$

où $\textstyle f_{\mathbf {k} }^{(1)}$ est la distribution de Bose-Einstein en repère mobile à la vitesse de dérive $\textstyle \mathbf {u} (\mathbf {x} ,t)$ :

$f_{\mathbf {k} }^{(1)}={\frac {1}{e^{{\frac {\hbar }{k_{B}T}}\left(\omega _{\mathbf {k} }+\mathbf {u} \cdot \mathbf {k} \right)}-1}}\approx n_{\mathbf {k} }^{(0)}+{\frac {\hbar \mathbf {u} \cdot \mathbf {k} }{k_{B}T}}{\frac {e^{\frac {\hbar \omega _{\mathbf {k} }}{k_{B}T}}}{\left(e^{\frac {\hbar \omega _{\mathbf {k} }}{k_{B}T}}-1\right)^{2}}}$

Au contraire de la distribution de Bose-Einstein cette distribution n'est pas isotrope. $\textstyle f_{\mathbf {k} }^{(1)}$ correspond à une distribution de Bose-Einstein observée dans un repère lorentzien se déplaçant à la vitesse $\textstyle \mathbf {u}$ par rapport au solide. Ce terme participe à la propagation et donc indirectement à la résistivité thermique.

Démonstration

La thermodynamique des milieux hors équilibre montre comment rechercher la solution du problème en introduisant le minimum d'information pour sa résolution. La solution correspond au maximum d'entropie à énergie fixée et, s'il y a lieu, à quantité de mouvement fixée. Cette entropie s'écrit :

${\begin{array}{rcl}S&=&k_{B}\rho \int _{\mathbf {k} }\left[\left(n_{\mathbf {k} }+1\right)\log {(n_{\mathbf {k} }+1)}-n_{\mathbf {k} }\log {n_{\mathbf {k} }}\right]\mathrm {d} \mathbf {k} \\[0.5em]&=&k_{B}\int _{\mathbf {k} }\left[\left(\rho +f_{\mathbf {k} }\right)\log {\left(1+{\frac {f_{\mathbf {k} }}{\rho }}\right)}-f_{\mathbf {k} }\log {\frac {f_{\mathbf {k} }}{\rho }}\right]\mathrm {d} \mathbf {k} \end{array}}$

La ou les contraintes s'écriront à l'aide de multiplicateurs de Lagrange notés $\beta$ et $\mathbf {u}$. Dans le cas des processus résistifs on maximise la quantité suivante,

$\textstyle S+\beta \left(\int _{\mathbf {k} }\hbar \omega _{\mathbf {k} }f_{\mathbf {k} }\mathrm {d} \mathbf {k} -E\right)$

l'énergie volumique $E$ étant fixée, ce qui conduit à la valeur habituelle :

$\textstyle \beta ={\frac {1}{k_{B}T}}$

et donc à la distribution de Bose-Einstein.

Dans le cas des processus normaux la quantité à maximiser est :

$S+\beta \left(\int _{\mathbf {k} }\hbar \omega _{\mathbf {k} }f_{\mathbf {k} }\mathrm {d} \mathbf {k} -E\right)+\mathbf {u} \cdot \left(\int _{\mathbf {k} }\hbar \mathbf {k} f_{\mathbf {k} }\mathrm {d} \mathbf {k} -\mathbf {p} \right)$

la quantité de mouvement $\mathbf {p}$ étant également fixée. Ceci conduit à la distribution $\textstyle f_{\mathbf {k} }^{(1)}$ donnée ci-dessus. La vitesse de dérive $\mathbf {u}$ n'est pas calculable analytiquement sauf cas particulier.

Méthodes aux moments

À une époque plus récente les physiciens se sont intéressés aux méthodes utilisées en mécanique des fluides et en transfert radiatif, particulièrement aux méthodes basées sur l'utilisation des moments de la fonction de distribution comme la méthode de Grad ou la méthode MN. On peut remarquer que, dans l'hypothèse du continu, $\textstyle {\frac {\mathbf {k} }{k}}$ est une distribution angulaire et donc que la quantité $I_{\mathbf {k} ,s}=\textstyle \hbar \omega _{\mathbf {k} }v_{\mathbf {k} }f_{\mathbf {k} ,s}{\frac {\mathbf {k} }{k}}$ est l'exacte contrepartie de la luminance spectrale du rayonnement de photons.

Moments de l'équation

Prenons l'exemple simple d'un milieu isotrope dans lequel on ne retient que les ondes acoustiques dont la vitesse de propagation $\textstyle v$ est indépendante du nombre d'onde (approximation $\textstyle k\to 0$) et connue. La relation de dispersion se réduit à $\textstyle \omega _{\mathbf {k} }=vk$. Dans ce qui suit on ignore les polarisations.

On définit les moments suivants ($\textstyle \otimes$ est le produit tensoriel) :

$E=\int _{\mathbf {k} }\hbar kvf_{\mathbf {k} }\mathrm {d} \mathbf {k} \quad$ énergie volumique ;

$\mathbf {F} =\int _{\mathbf {k} }\hbar kv^{2}f_{\mathbf {k} }{\frac {\mathbf {k} }{k}}\mathrm {d} \mathbf {k} \quad$ flux d'énergie, et sa valeur réduite $\textstyle \mathbf {f} ={\frac {\mathbf {F} }{vE}}$ ;

${\mathsf {P}}=\int _{\mathbf {k} }\hbar kvf_{\mathbf {k} }{\frac {\mathbf {k} }{k}}\otimes {\frac {\mathbf {k} }{k}}\mathrm {d} \mathbf {k} \quad$ tenseur de pression et sa valeur normalisée (tenseur d'Eddington) $\textstyle {\mathsf {D}}={\frac {\mathsf {P}}{E}}$.

Ces quantités sont les moments de Hausdorff d'une quantité non-négative sur la sphère unité. On montre que $\textstyle \mathbf {F}$ est un vecteur propre de $\textstyle {\mathsf {P}}$ d'où :

${\mathsf {D}}\mathbf {f} =\chi \mathbf {f}$

où $\textstyle \chi$ est le facteur d'Eddington.

Le tenseur d'Eddington s'exprime en fonction du facteur d'Eddington par :

${\mathsf {D}}={\frac {1-\chi }{2}}\,{\mathsf {I}}+{\frac {3\chi -1}{2}}{\frac {\mathbf {f} }{f}}\otimes {\frac {\mathbf {f} }{f}}\,,\quad {\frac {1}{3}}\leq \chi \leq 1$

où $\textstyle {\mathsf {I}}$ est le tenseur unité. $\textstyle \chi ={\frac {1}{3}}$ correspond à l'isotropie et $\textstyle \chi =1$ à un faisceau parallèle. La solution générale possède la symétrie de révolution autour de $\textstyle \mathbf {f}$.

On obtient les équations aux moments de la façon suivante :

• On multiplie l'équation de Bolzmann-Peierls par l'énergie du phonon $\textstyle \hbar kv$ et on intègre sur k :

$\int _{\mathbf {k} }{\frac {\partial f_{\mathbf {k} }}{\partial t}}\hbar kv\mathrm {d} \mathbf {k} +\int _{\mathbf {k} }\hbar kv^{2}{\frac {\mathbf {k} }{k}}\cdot \nabla f_{\mathbf {k} }\,\mathrm {d} \mathbf {k} =\int _{\mathbf {k} }S_{\mathbf {k} }\hbar kv\mathrm {d} \mathbf {k} =0$

Le second terme de cette équation est nul du fait de la conservation de l'énergie de tous les processus.

Par identification on obtient :

${\frac {\partial E}{\partial t}}+\nabla \cdot \mathbf {F} =0$

• On multiplie l'équation de Bolzmann-Peierls par $\textstyle \hbar kv{\frac {\mathbf {k} }{k}}$ et on intègre sur k :

$\int _{\mathbf {k} }{\frac {\partial f_{\mathbf {k} }}{\partial t}}\hbar kv{\frac {\mathbf {k} }{k}}\mathrm {d} \mathbf {k} +\int _{\mathbf {k} }\hbar kv^{2}{\frac {\mathbf {k} }{k}}\cdot \nabla f_{\mathbf {k} }{\frac {\mathbf {k} }{k}}\mathrm {d} \mathbf {k} =\int _{\mathbf {k} }S_{\mathbf {k} }\hbar kv{\frac {\mathbf {k} }{k}}\mathrm {d} \mathbf {k}$

Dans le second membre les processus élastiques conduisent à un bilan nul. On suppose de plus que tous les processus résistifs sont décrits par une approximation BGK à une seule constante de temps $\textstyle \tau$. Alors par identification on obtient :

${\frac {\partial \mathbf {F} }{\partial t}}+v^{2}\nabla \cdot {\mathsf {P}}=-{\frac {\mathbf {F} }{\tau }}$

Comme tout système aux moments il faut à présent faire une hypothèse pour « fermer » (compléter) le système par nature incomplet puisque l'on a deux équations pour les trois inconnues $\textstyle E,\,\mathbf {f} ,\,\chi$.

Fermeture entropique

On va comme ci-dessus maximiser l'entropie en respectant les moments définis ci-dessus. Cette méthode constitue la méthode M₁. La solution est analytique :

$\chi ={\frac {3+4f^{2}}{5+2{\sqrt {4-3f^{2}}}}}$

Le système à présent complet est hyperbolique et l'on peut utiliser pour le résoudre les méthodes très efficaces qui ont été développées pour la mécanique des fluides.

On peut reconstruire la fonction de distribution qui n'apparaît pas explicitement dans le problème :

$n_{\mathbf {k} }={\frac {1}{e^{{\frac {\hbar kv}{k_{B}T}}\left[g_{1}(f)+g_{2}(f)\mathbf {f} \cdot \mathbf {k} \right]}-1}}$

où g₁ et g₂ sont analytiques. Cette expression comporte une « vitesse de dérive » $\textstyle \,\mathbf {u} =g_{2}(f)\,\mathbf {f} \,$ qui s'interprète comme la manifestation locale du gradient de peuplement de phonons dans chaque direction, avec un maximum dans le sens du flux d'énergie.

L'analyse asymptotique du système aux dérivées partielles permet d'exhiber une limite diffusive :

$\mathbf {F} =-{\frac {4avT^{3}\tau }{3}}\nabla T$

où a est l'équivalent de la constante radiative modifiée pour prendre en compte le nombre de polarisations possibles (3 au lieu de 2 pour les photons) $\textstyle a={\frac {12\pi ^{5}k_{B}^{4}}{15h^{3}v^{3}}}$

En identifiant le coefficient avec la conductivité thermique $\textstyle \lambda$ on obtient le temps caractéristique :

$\tau ={\frac {3\lambda }{4avT^{3}}}$

Résolution

Transfert de phonons dans une lame de silicium . Comparaison calcul-mesure. « Fourier » désigne le calcul de conduction classique.

Au final on est ramené à un système aux dérivées partielles hyperbolique comportant une équation scalaire sur E et une équation vectorielles sur $\textstyle \mathbf {f}$ pour lequel il existe de nombreux solveurs numériques efficaces développés dans le cadre de la mécanique des fluides.

Il y a cependant une difficulté à écrire les conditions aux limites au bord du matériau. On remarquera d'abord que la notion de phonon repose sur l'hypothèse d'un milieu infini. On supposera toutefois que l'équation de Boltzmann-Peierls reste valide au voisinage de la surface. En l'absence de transmission (on écarte le problème de l'interface entre deux matériaux) on écrit les conditions aux limites comme on le fait dans le cas de la couche de Knudsen :

• la réflexion est soit spéculaire, soit diffuse. On peut également superposer les deux types de réflexion en les pondérant par un coefficient d'« accomodation » ;
• on peut également imposer une distribution d'équilibre à la paroi : cette condition, difficile à justifier physiquement, est utilisée pour étudier des problèmes spécifiques.

Cette approche est également valide pour l'interface entre deux matériaux (résistance de Kapitza).