Les travaux de recherches que j'ai réalisés depuis mon master (2014), portent sur la physique statistique hors-équilibre dans des contextes biologiques, hydrodynamiques et électromagnétiques. Mon approche consiste principalement à résoudre analytiquement des problèmes physiques en développant des expressions asymptotiques, tout en réalisant des résolutions numériques d'équations et des simulations stochastiques, et en confrontant ces résultats à des expériences. Mes contributions ont porté sur sept thèmes majeurs:
1. Réduction du temps de recherche des lymphocytes NK en présence de spectateurs [P1]. Les lymphocytes NK cherchent à éliminer les cellules infectées, et rencontrent des cellules spectatrices favorisant leur migration pendant leur recherche. J'ai développé un modèle numérique pour comprendre les expériences in vitro, montrant que le temps de recherche est réduit en présence de ces cellules spectatrices.
2. Schémas d'approximation quantitatifs pour la transition vitreuse [P2]. La transition vitreuse est une transition du premier ordre dépendant de la cinétique du refroidissement, ne se réalisant pas à une température donnée (à pression fixée), au contraire de la transition liquide-solide. J'ai réalisé un développement systématique de l'équation d'état autour de sa solution exacte en dimension infinie, permettant d'obtenir les propriétés thermodynamiques de la transition vitreuse (pour $d\ge3$) et de construire le diagramme de phase.
3. Dispersion dans les milieux complexes [P3,P4,P5,P6,P7]. La dispersion permet de comprendre l’étalement d’un nuage de particules initialement proches dans les milieux complexes. Dans ces milieux, la distribution de Boltzmann n’est pas atteinte dans le régime des temps longs et cet étalement est caractérisé par la diffusivité effective. J'ai étudié la dispersion dans des microcanaux en m'intéressant uniquement à son lien avec la géométrie de confinement, permettant de caractériser trois régimes de dispersion analytiquement. Puis, j'ai analysé la dispersion dans un réseau d'obstacles, identifiable à celle dans les microcanaux. Enfin, j'ai regardé cette dispersion en ajoutant un potentiel attractif sur la surface des obstacles, qui est étrangement accélérée par rapport à des obstacles réfléchissants malgré l'ajout d'un piège énergétique.
4. Vortex browniens créés par les pièges optiques [P8,P9,P10]. Les pièges optiques sont souvent supposés harmoniques, or ils sont généralement anharmoniques et la pression de radiation du laser crée une force non conservative. Des courants stationnaires sont alors présents formant des vortex browniens. Je me suis intéressé au régime sous-amorti atteint dans des expériences menées en parallèle pour des basses pressions. J'ai réalisé une étude analytique et numérique du courant stationnaire et la densité spectrale, vérifiées expérimentalement.
5. Temps d'échappement d'une cellule et organisation du cytosquelette [P11,P12]. Le transport intracellulaire de diverses particules-cargos vers la synapse immunologique est crucial pour le bon fonctionnent des cellules. Le temps de ce transport peut être minimisé en fonction de l'organisation du cytosquelette, et en particulier de la taille du cortex d'actine. J'ai étudié le temps de premier passage dans la limite d'une petite synapse et d'un fin cortex pour une géométrie simplifiée du cytosquelette pour des particules passives, et montré qu'il pouvait être minimisé en présence d'un mécanisme forçant les particules à aller vers la membrane cellulaire.
6. Mouvement collectif de spins avec un alignement local [P13,P14,P16,P15,P17,P18]. Les mouvements collectifs de grands amas de particules sont largement observés dans la nature. Il s'agit d'un phénomène hors-équilibre décrit par une transition de phase liquide-gaz sans région supercritique. J'ai étudié les modèles de Potts et d'horloge actifs pour lesquels les états internes des particules, correspondant à leurs directions de mouvement, s’alignent localement. J'ai analysé la transition de phase, faisant apparaître une réorientation de la phase liquide et un état bloqué en présence d'interactions. J'ai également étudié des modèles à plusieurs espèces pouvant faire intervenir des interactions d'anti-alignement et des non-réciprocités.
7. Sédimentation de la matière active: étude des courants stationnaires [P19]. La matière active présente comportements hors équilibre complexes tels que l'accumulation sur les parois,et accompagnés de courants stationnaires émergeant en l'absence d'alignement. J'ai étudié ces courants stationnaires sous l’effet de la gravité en présence d’une paroi, en considérant un modèle simplifié de particules browniennes actives. J'ai montré que l'accumulation de particules sur les parois est associée à la présence de plusieurs vortex.
Ce travail a été effectué au cours de mon stage de 5 mois du master 1 sous la direction de H. Rieger et K. Schwarz à Sarrebruck (Allemagne), pendant le premier semestre 2014. Il a donné lieu à une publication [P1].
Pendant ce stage, j'ai modélisé le processus de recherche des cellules tueuses naturelles (CTN) ou lymphocytes NK en présence de cellules spectatrices, en collaboration avec l'équipe de physique expérimentale de B. Qu à Hombourg (Allemagne). Les CTN cherchent à éliminer les cellules infectées par un pathogène (virus) ou tuméfiées, appelées ci-après cellules cibles, pour éviter la propagation de l'infection. Pendant leur recherche, elles rencontrent d'autres types de cellules non-cibles (saines) qui sont spectatrices mais qui augmentent l'efficacité et la migration de ces CTN. Cela est dû à la production de peroxyde d'hydrogène (H2O2) par ces cellules spectatrices. Les expériences in-vitro ont montré que la vitesse et la persistance des CTN étaient accrues en leur présence, sans destruction de leur capacité létale.
Nous avons alors développé des modèles diffusifs bidimensionnels, dans des domaines périodiques, permettant de comprendre les observations expérimentales. Les CTN ont été représentées comme des cellules browniennes discoïdes cherchant les cellules cibles immobiles instantanément détruites lorsqu'elles sont trouvées. Les cellules spectatrices, immobiles également, augmentent la diffusivité des CTN localement dans un rayon fini. Nous avons utilisé la méthode de Monte Carlo cinétique [1.1,1.2] d'une part, ainsi qu'un modèle sur réseau d'autre part pour reproduire cette dynamique stochastique numériquement. Ces simulations numériques montrent que le temps moyen pour localiser les cellules cibles, défini comme le temps de demi-vie des cibles, est réduit en présence des cellules spectatrices, autant par leur nombre que leur rayon d'action. Cet impact est supérieur à la présence de simples obstacles, même si ceux-ci favorisent déjà la recherche en diminuant la surface accessible. Notre étude numérique montre donc que la recherche des CTN est plus efficace en présence des cellules spectatrices, reproduisant parfaitement les observations des expériences in-vitro.
Ce travail a été effectué au cours de mon stage de 2 mois du master 2 sous la direction de F. Zamponi au LPTENS de Paris, pendant l'hiver 2015. Il a donné lieu à ma première publication [P2].
Pendant ce stage, j'ai étudié les propriétés thermodynamiques de la transition vitreuse en développant un schéma d'approximation en partant de leurs solutions exactes en dimension infinie. La transition vitreuse est une transition du premier ordre aléatoire, appartenant à la classe d'universalité des verres de spin [2.1,2.2]. La phase vitreuse n'est donc pas unique, au contraire du cristal par exemple. Le diagramme de phase est alors extrêmement complexe, exhibant plusieurs états vitreux caractérisés par différentes propriétés thermodynamiques (par exemple l'équation d'état) qui dépendent généralement de la cinétique de la transition vitreuse.
Comme plusieurs modèles statistiques, les propriétés thermodynamiques de la transition vitreuse sont solvables en dimension infinie avec une structure de champ moyen [2.3]. Nous avons alors réalisé un développement systématique autour de cette solution en dimension infinie, similaire à un développement haute température ou basse densité. En utilisant la méthode des répliques [2.4], pour m copies du système original, nous avons alors dérivé des expressions analytiques pour l'entropie et l'équation d'état de la phase vitreuse en toutes dimensions d. Ces équations ne dépendent que de la thermodynamique et de la structure du liquide, et en particulier de la corrélation de paire gliq et de l'entropie Sliq, et ont une structure comparable à celles obtenues avec la théorie du couplage de modes (TCM) [2.5] valable uniquement pour d=3. Seules des solutions approximatives de gliq et Sliq existent dont les plus connues sont nommées Hyper-Netted-Chain et Percus-Yevick et dérivées à partir de l'équation du viriel [2.6].
À partir de l'équation d'état et de l'entropie (aussi appelée complexité Σeq) du verre, ainsi que des valeurs numériques de la corrélation de paire du liquide, nous avons calculé numériquement, en dimension d≥3, le facteur de non-ergodicité et les densités: (i) de la transition dynamique φd pour laquelle le liquide devient infiniment visqueux et plusieurs état vitreux apparaissent, et à partir de laquelle l'équation d'état pour m=1 admet des solutions; (ii) de la transition de Kauzmann φK pour laquelle le nombre d'état vitreux devient subexponentiel avec une crise entropique, Σeq=0 pour m=1, correspondant à la transition vitreuse idéale du second ordre; (iii) de la transition de jamming (m=0) allant de la densité de seuil φth à la densité de conditionnement compact du verre φGCP, correspondant respectivement aux transitions dynamique et de Kauzmann. Pour la dimension d=3, nous avons obtenu des valeurs numériques (φd≈0.53, φK≈0.62, φth≈0.45 et φGCP≈0.68) quasi-similaires aux précédentes études [2.7,2.8], basées notamment sur la TCM. Nous avons également obtenus des valeurs pour toutes dimensions d≥3. Ce travail a donc permis de relier les propriétés thermodynamiques de la transition vitreuse en dimension d=3 et en dimension infinie.
Ce travail a été effectué au cours de ma thèse [PhD] réalisée sous la direction de D. S. Dean et T. Guérin au LOMA de Bordeaux, entre Octobre 2015 et Septembre 2018. Il a donné lieu à cinq publications [P3,P4,P5,P6,P7].
Pendant la première partie de ma thèse, j'ai analysé la dispersion de particules browniennes au sein de certains milieux complexes: les microcanaux et les milieux hétérogènes (réseaux d'obstacles). La diffusion dans les milieux simples, sans confinement entropique, est parfaitement caractérisée par les travaux réalisés au cours du XIXe siècle posant les bases du mouvement brownien, des phénomènes diffusifs et de la mécanique statistique; et du début du XXe siècle montrant que le mouvement brownien est un processus diffusif et posant les bases de la dynamique stochastique.
Or, la diffusion dans les milieux complexes est moins comprise, bien que largement étudiée. L'étude de la dispersion permet de comprendre l'étalement d'un nuage de particules browniennes initialement proches au cours du temps dans les milieux complexes. C'est un phénomène qui apparaît dans de nombreux contextes, comme le mélange ou le tri de particules, l'étalement de tâches de polluants, ou la cinétique des réactions chimiques. Nous nous sommes intéressés principalement à la dispersion dans les milieux confinés (microcanaux et réseaux d'obstacles) importante dans les systèmes biologiques [3.1], la microfluidique [3.2,3.3] ou les milieux poreux [3.4]. Pour ces systèmes, la distribution de Boltzmann n’est pas atteinte dans le régime des temps longs. Cependant, la taille typique de l'étalement des particules browniennes, donnée par le déplacement quadratique moyen, augmente au cours du temps et est caractérisée par la diffusivité effective $D_e$ inférieure à la diffusivité microscopique $D_0$, principal effet du confinement. Au cours de cette étude nous avons développé plusieurs expressions analytiques, valables dans différents régimes asymptotiques et vérifiées par des simulations numériques.
Dans un premier temps, nous avons étudié la dispersion de particules browniennes dans les microcanaux périodiques, en nous intéressant uniquement au lien entre la géométrie du confinement et la dispersion. En effet, même pour des particules browniennes ponctuelles et sans interactions, dans un fluide au repos, il est bien connu que le confinement peut induire une réduction significative du ratio $D_e/D_0$. Les canaux considérés sont symétriques et bidimensionnels avec des parois situées en $y=\pm h(x)$ de période $L$, dans le système de coordonnées cartésiennes $(x,y)$. La diffusivité effective est alors reliée au déplacement quadratique moyen dans l'axe du canal $\langle x^2 \rangle \simeq 2 D_e t$ dans la limite des temps longs. La première théorie quantitative a été introduite par Jacobs [3.6], basée sur une réduction de dimensionnalité de l’espace et considérée aujourd'hui comme standard sous le nom de Fick-Jacobs (FJ). Elle permet de comprendre l'origine de cette diminution de la diffusivité. Une entropie effective $s(x) = k_B \ln h(x)$ peut être définie en comptant le nombre de configurations à la position $x$, en supposant un équilibre rapide dans la direction transverse. Les régions étroites et larges agissent alors respectivement comme des barrières et des pièges entropiques, ralentissant la dispersion des particules. La diffusivité effective est alors calculée en considérant les solutions exactes des processus diffusifs unidimensionnels [3.7]. De nombreuses études ont proposé des améliorations à cette approximation au cours des dernières décennies [3.8,3.9,3.10,3.11], en se basant pour la plupart sur une réduction de dimensionnalité introduisant une diffusivité unidimensionnelle $D(x)$, bien qu'elle mène généralement à des processus non-markoviens.
Tout d'abord, nous avons utilisé les formules de Kubo de T. Guérin et D. S. Dean [3.12,3.13], dérivées à partir de l'équation de Fokker-Planck et donnant une expression exacte de la diffusivité effective à partir d'une fonction auxiliaire $f(x,y)$, solution de l'équation de Laplace dans le canal. Ces formules permettent d'obtenir des résultats numériques exacts en déterminant la fonction auxiliaire par la méthode des éléments finis avec le logiciel FreeFem++ [3.5]; ainsi que des résultats analytiques asymptotiques. En considérant la fonction analytique $w(x+iy)= {\rm Re} f(x,y)$, nous avons montré que la diffusivité effective satisfait alors $${\rm Im} w\left(x+ih(x)\right) = h(x)-C, \qquad \frac{D_e}{D_0}= \frac{C}{\langle h \rangle}, \tag{3.1}$$ où la constante $C$ peut être déterminée grâce à la condition de périodicité. La notation $\langle h \rangle$ représente la moyenne spatiale de la hauteur $h(x)$ sur une période. Il s'agit d'un formalisme très compact, simplifiant les calculs.
Nous avons alors réalisé une classification des différents régimes de dispersion en fonction de la géométrie du canal. En définissant la hauteur minimale $a$ et la variation de hauteur $H$, la géométrie du canal est caractérisée par deux nombres sans dimensions : $\varepsilon = a/L$ et $\xi = H/a$ avec un profil donné par la hauteur $h(x) = a + H g(x/L)$ où $g(u)$ est une fonction $1$-périodique à valeurs dans l'intervalle $[0,1]$. Un exemple est donné dans l'encart de la Fig. 1a, pour lequel une trajectoire brownienne a été tracée. Dans la limite des canaux étroits ($\varepsilon \ll 1$), en définissant $h(x) = \varepsilon \zeta(x)$, nous avons dérivé à l'aide d'une méthode perturbative l'expression: $$\frac{D_e}{D_0} = \frac{1}{\langle h\rangle\langle h^{-1} \rangle } \left[1 - \frac{\varepsilon^2}{3} \frac{\langle \zeta'^2/\zeta\rangle}{\langle \zeta^{-1} \rangle} + {\cal O}(\varepsilon^4)\right]. \tag{3.2}$$ L'ordre dominant est donné par la diffusivité de FJ qui nous permet de distinguer deux types de canaux en fonction du comportement près de la hauteur minimale du canal, dans la limite de fortes ondulations $\xi \gg 1$. En supposant que la hauteur s'écrit $h(x) = a + \Lambda |x|^\nu$ près de son minimum, la diffusivité de FJ se réécrit en fonction du temps de premier passage à travers un entonnoir [3.14] pour un exposant $\nu>1$, et ne dépend pas de la géométrie du reste du canal. Pour un exposant $\nu < 1$, la diffusivité de FJ tends vers une constante dans la limite $\xi \rightarrow \infty$. Ces deux régimes de dispersion, contrôlés par la géométrie près de la hauteur minimale, définissent respectivement les canaux lisses et rugueux. De plus, nous avons montré que les résultats dérivés en considérant une réduction de dimensionnalité [3.8,3.9,3.10] sont bien valides, au moins jusqu'à l'ordre $\varepsilon^6$.
Dans la limite opposée des canaux larges ($\varepsilon \gg1$), qui n'avait jamais été abordée dans la littérature, la dispersion est dominée par les particules se situant dans la région $|y| < a$ où les particules contribuent au déplacement longitudinal alors que celles situées dans la région $|y| > a$ sont majoritairement piégées par les parois du canal. L'ordre dominant de la diffusivité effective est donc obtenu en considérant que les excursions dans les régions lentes n'y contribuent pas. Il s'agit du résultat observé pour la dispersion dans des peignes [3.15]. À partir de notre formalisme complexe et plusieurs transformations conformes du domaine périodique, nous avons dérivé l'ordre sous-dominant de la diffusivité effective qui ne dépend pas du profil du canal et faisant intervenir une constante universelle $\ln 2/\pi$: $$\frac{D_e}{D_0} = \frac{a}{\langle h\rangle} \left[ 1 + \frac{\ln 2}{\pi \varepsilon} + {\cal O}\left(\varepsilon^{-2}\right) \right]. \tag{3.3}$$
Dans les canaux de taille intermédiaire (pour des valeurs de $\varepsilon$ intermédiaires), nous avons défini un approximant de Padé permettant de relier les résultats asymptotiques pour $\varepsilon \ll1$ et $\varepsilon \gg 1$. Celui-ci approxime assez bien les valeurs numériques dans le régime intermédiaire pour les canaux lisses, mais présente quelques défauts pour les fortes ondulations ($\xi \gg1$). La dispersion est alors contrôlée par la présence de petites ouvertures. En considérant perturbativement les petites ouvertures [3.16,3.17], nous avons montré que la diffusivité effective s'écrit $$\frac{D_e}{D_0} \simeq \frac{\pi L}{2 \langle h \rangle} \frac{1}{\ln(2\kappa/\varepsilon)}, \tag{3.4}$$ où $\kappa$ dépend de la pseudo-fonction de Green [3.18] du domaine sans ouverture. Dans la limite $H \gg L$, cette constante a pour valeur $\kappa= 2/\pi$. À partir de ce résultat, la diffusivité effective peut s'écrire en fonction du temps de premier passage pour aller d'une petite ouverture, considérée réfléchissante, à la petite ouverture suivante [3.17]. Ce temps diffère ici du temps moyen de premier passage global pour atteindre une petite ouverture, par la correction du terme logarithmique donnée par $\kappa$. Nous avons alors obtenu trois régimes de dispersion contrôlés par la géométrie des canaux : canaux étroits, canaux larges et canaux présentant des petites ouvertures. Sur la Fig. 1a, nous avons représenté la diffusivité effective en fonction des paramètres $\varepsilon$ et $\xi$ pour le canal tracé en encart, en comparant les solutions numériques des formules de Kubo et les expressions analytiques obtenues pour ces trois régimes. Le diagramme de validité de ces expressions est montré sur la Fig. 1b, avec des lignes de transitions en $\varepsilon =1$ et $\varepsilon \xi^{1/\nu}=1$ pour les canaux lisses ($\nu ≥ 1$) ou $\varepsilon \xi=1$ pour les canaux rugueux ($\nu < 1$).
Enfin, nous avons abordé la dispersion dans les canaux étroits présentant des discontinuités de profil. Le résultat asymptotique développé pour les canaux étroits [Eq. (3.2)] n'est valable qu'à l'ordre dominant, dans la mesure où la dérivée de $\zeta(x)$ diverge. En utilisant notre formalisme complexe et en réalisant des transformations conformes du domaine, la diffusivité effective s'écrit pour $N$ discontinuités situées en $x=x_i$: $$\frac{D_e}{D_0} = \frac{1}{\langle h\rangle\langle h^{-1}\rangle} \left[ 1- \frac{\varepsilon}{L \langle \zeta^{-1} \rangle}\sum_{i=1}^N \gamma(\lambda_i) + {\cal O}(\varepsilon^2) \right], \tag{3.5}$$ où $\lambda_i$ est un paramètre relié à la $i^e$ discontinuité. Le piégeage des particules est donc plus important proche des discontinuités (terme d'ordre $\varepsilon$) que dans les sections continues (terme d'ordre $\varepsilon^2$), faisant diminuer tous deux la valeur de la diffusivité effective. La fonction $\gamma(\lambda)$ dépend de la géométrie de la discontinuité, et nous avons dérivé deux expressions exactes: $$\gamma_d(\lambda)=\frac{1+\lambda^2}{\pi\lambda} \ln \left\vert \frac{1+\lambda}{1-\lambda} \right\vert - \frac{2}{\pi} \ln \left\vert \frac{4\lambda}{1-\lambda^2}\right\vert; \qquad \gamma_c(\lambda)=-\frac{4}{\pi}\ln \left(\sin \frac{\pi \lambda}{2}\right), \tag{3.6}$$ où $\gamma_d(h_+/h_-)$ est valable pour une discontinuité simple avec une hauteur passant de $h_-$ à $h_+$ (canal 1 de la Fig. 1c), et $\gamma_c(h_0/h_-)$ pour des cloisons réduisant la hauteur au seul point de discontinuité $h_0 < h_-=h_+$ (canal 2 de la Fig. 1c). Sur la Fig. 1c, nous avons représenté la valeur de $\gamma$ obtenue à partir de la diffusivité effective, en comparant la solution numérique des formules de Kubo et des deux expressions analytiques dérivées. Nos résultats sont généralisables à d'autres formes de discontinuité. Il est à noter que ces fonctions $\gamma(\lambda)$ sont exactes et peuvent être reliées aux taux de piégeages introduits au niveau des discontinuités en considérant une réduction de dimensionnalité, calculés approximativement dans [3.19,3.20].
Nous avons donc décrit les régimes de dispersion pour toutes les géométries possibles d'un canal, et proposé une expression asymptotique de la diffusivité effective pour chacun d'eux en les reliant à celles utilisées dans la littérature. Ces expressions analytiques peuvent cependant être améliorées en calculant les ordres asymptotiques suivants.
Dans un second temps, nous avons analysé la dispersion de particules browniennes ponctuelles dans un réseau périodique d'obstacles sphériques attractifs. La dispersion dans les milieux hétérogènes apparaît dans de nombreux contextes physiques, tels que les milieux poreux, les suspensions colloïdales ou les cellules vivantes, et permet de décrire les phénomènes de diffusion anormale [3.22]. La seule présence des obstacles sphériques (de rayon $R$) implique un ralentissement de la dispersion due à un piégeage entropique des particules. Il s’agit d’un mécanisme de piégeage entropique identique à celui mentionné dans la partie A. En dimension $d=2$, la dispersion dans un réseau d'obstacles est alors analogue à celle dans un canal, de même périodicité $L$ et de hauteur $h(x)=L/2 - \sqrt{R^2-x^2}$ [3.21], comme illustré sur la Fig. 2a. Cela nous permet de déterminer deux régimes spécifiques de dispersion dépendant uniquement de la fraction volumique des obstacles $\varphi = \pi R^2/L^2$. Dans la limite diluée ($\varphi \ll 1$), en utilisant le formalisme complexe développé pour les canaux [Eq. (3.1)], la diffusivité effective s'écrit $D_e/D_0 \simeq 1/(1+\varphi)$. Ce résultat est équivalent à celui de Maxwell, démontré pour la conductivité effective d’un milieu électrique [3.23,3.24]. Il est exact jusqu'à l'ordre $\varphi^4$ et est valide pour des fractions volumiques $\varphi < 0.4$. Dans la limite encombrée ($\varphi \simeq \pi/4$), l'approximation de FJ peut être appliquée. En remarquant que les formules de Kubo ne font intervenir que la normale à la sphère $n_x = x/R$, nous avons proposé une diffusivité unidimensionnelle $D(x) \simeq D_0 (1 - n_x^2/3)$ en remplaçant $h'(x)$ par $n_x$, équivalentes lors de la réduction de dimensionnalité [3.8,3.9,3.10]. Cette expression donne une meilleure approximation de la diffusivité effective que celles développées précédemment [3.21] et est valable pour des fractions volumiques à $\varphi > 0.67$. Sur la Fig. 2b, nous avons tracé la diffusivité effective en fonction de la fraction volumique $\varphi$, en comparant les solutions numériques des formules de Kubo (exactes) et les expressions analytiques mentionnées ci-dessus.
Ensuite, nous avons étudié l'effet d'un potentiel attractif $V(r)$ localisé sur la surface des sphères, sur la dispersion dans les réseaux d'obstacles. Les particules browniennes sont piégées dans les minima locaux de ce potentiel, entraînant un ralentissement bien connu et caractérisé par le temps moyen d'échappement du minimum de potentiel $t_{\rm esc} \sim \exp(\beta\Delta E)$, suivant la loi d'Arrhénius avec $\Delta E$ la barrière d'énergie à franchir et $\beta$ l'inverse de la température. Cependant, l'effet du couplage des deux barrières entropique (obstacles) et énergétique (potentiel attractif) sur la dispersion n'est pas si simple. Récemment, il a été mis en évidence par des simulations numériques [3.25] que la dispersion pouvait être accélérée pour un potentiel attractif, tout en respectant l'inégalité $D_e ≤D_0$. Ce résultat est contre-intuitif par le fait que la présence d'un minimum de potentiel puisse dépiéger une particule brownienne du confinement entropique des sphères. Cet effet pourrait cependant être compris si les surfaces des obstacles étaient suffisamment proches pour que les particules puissent passer de l'une à l'autre sans explorer tout l'espace. Mais ce phénomène apparaît également dans la limite diluée! Le potentiel que nous avons considéré est attractif $V(r)=-V_a <0$ sur $R<r<R_a$, présentant une barrière $V(r)=V_r>0$ sur $R_a<r<R_r$ et nul au delà (voir Fig. 2c). À l'aide des formules de Kubo [3.12,3.13], nous avons dérivé une expression de la diffusivité effective dans la limite diluée, cohérente avec celle obtenue pour un puits de potentiel carré ($V_r=0$) [3.26,3.27,3.28]. Nous avons alors remarqué que la dispersion était optimisée dans la limite de faible portée du potentiel $\alpha_a = (R_a-R)/R\ll 1$ et $\alpha_r = (R_r-R_a)/R\ll 1$, pour laquelle la diffusivité effective s'écrit en dimension $d=2$: $$\frac{D_e}{D_0} = 1- \frac{2(\tau+1)\xi^2 - (1-\tau)\xi +1}{(\tau+1)\xi+1} \varphi + {\cal O}(\varphi^2), \tag{3.7}$$ pour $\xi = \alpha_a \exp(\beta V_a)$ et $\tau = \alpha_a \alpha_r \exp[\beta (V_a+V_r)]$ fixés. Cette expression est tracée sur la Fig. 2d par rapport à la diffusivité effective maximale $D_e^*/D_0 \simeq 1 - 0.657 \varphi$. La dispersion est alors optimisée pour ce potentiel attractif si $D_e/D_0 > 1 - \varphi$, région définie par $\xi < 1/(1+\tau)$. De plus, les paramètres $\xi$ et $\tau$ peuvent être exprimés en fonction des temps de premier passage pour atteindre la surface $t_\Sigma$ en partant de l'extérieur (renormalisé par le volume) et pour s'échapper de celle-ci $t_{\rm esc}$ via les relations $\xi = t_{\rm esc}/t_{\Sigma}$ et $\tau = D_0 t_{\Sigma} / R^2$.
Enfin, nous avons étudié la dispersion des particules browniennes en considérant que la diffusion est médiée par la surface des sphères [3.29]. En considérant $D_\Sigma$ la diffusivité sur la surface, $D_{\rm B}$ la diffusivité dans le volume, $k$ le taux d'absorption sur la surface et $\lambda$ le taux de désorption, nous avons dérivé l'expression de la diffusivité effective à l'aide de formules de Kubo modifiées pour ce problème, par rapport aux équations de Fokker-Planck utilisées dans [3.29]. Cette expression revient à l'Eq. (3.7) dans la limite diluée pour $D_{\rm B}=D_\Sigma$, $\xi=k D_{\rm B}/\lambda R$ et $\tau=1/kR$. La dispersion dans un potentiel attractif localisé sur la surface des sphère est donc analogue à un processus de diffusion médiée par la surface des sphères et les taux d'absorption $k$ et de désorption $\lambda$ peuvent être reliés aux temps de premier passage définis précédemment: $t_{\rm esc} = \lambda^{-1}$ et ${t_\Sigma} = R/k D_{\rm B}$.
Nous avons donc décrit la dispersion dans un réseau d'obstacles en regardant l'effet entropique de ceux-ci d'une part, et le double effet entropique-énergique avec un potentiel attractif sur leur surface d'autre part. Le premier est parfaitement caractérisé tandis que le second a été analysé seulement dans la limite diluée montrant que la dispersion pouvait être accélérée par rapport à des obstacles réfléchissants.
Ce travail a été effectué au cours de ma thèse [PhD] réalisée sous la direction de D. S. Dean et T. Guérin au LOMA de Bordeaux, entre Octobre 2015 et Septembre 2018. Il a donné lieu à trois publications [P8,P9,P10].
Pendant la seconde partie de ma thèse, en collaboration avec Y. Louyer et Y. Amarouchene pour la partie expérimentale, j'ai analysé la dynamique stochastique de particules browniennes dans un piège optique en présence d'une force non conservative. Ces particules sont soumises à un piège purement énergétique, contrairement à celles de la partie 3 où elles étaient confinées par un piège entropique. Les pièges (ou pinces) optiques développés par Ashkin dans les années 1970 [4.2] permettent de manipuler des objets de petites dimensions (cellules vivantes, colloïdes, nanoparticules) à l'aide d'un faisceau laser et de mesurer l'intensité de petites forces. Ils ont permis de caractériser l'élasticité des brins d'ADN et vérifier des résultats de la physique statistique hors-équilibre et de la mécanique quantique. Le piège optique y est essentiellement supposé harmonique, or celui-ci est généralement anharmonique et la pression de radiation du laser crée une force non conservative dans sa direction longitudinale [4.3]. Le système est donc hors-équilibre et l'état stationnaire ne vérifie pas la distribution de Boltzmann, conduisant à la présence de courants. Les particules suivent alors une circulation toroïdale dans ce flux stationnaire, appelée vortex brownien [4.4] et observée à températures et pressions ambiantes, caractéristiques du régime suramorti [4.5].
Nous nous sommes intéressés au régime sous-amorti qui n'avait jamais été étudié, avec un coefficient de friction $\gamma$ faible pour les pressions atteintes dans les expériences menées en parallèle (i.e. $P<10 {\rm mbar}$). L'équation de Langevin sous-amortie pour la trajectoire ${\bf x}(t)$ de la particule brownienne en dimension $d=3$ est $$m \ddot {\bf x} = -\gamma \dot {\bf x} +{\bf F}({\bf x}) +\sqrt{2 k_B T\gamma}{\bf \xi}, \tag{4.1}$$ avec $m$ la masse de la particule brownienne et $T$ la température. ${\bf \xi}$ représente le bruit thermique, considéré comme un bruit blanc gaussien. La force ${\bf F}({\bf x})$ créée par le piège optique [4.3] est composée d'une partie conservative dérivant du potentiel harmonique $V_{\rm harm}({\bf x}) = \kappa \left(x^2+ y^2 + \eta z^2\right)/2$ avec une raideur $\kappa$ dans les directions transverses et $\eta \kappa$ dans la direction longitudinale ($\eta \simeq 0.2$); et d'une partie non-conservative ${\bf F}_{\rm scat}({\bf x}) = \varepsilon \kappa a [ 1-(x^2+y^2)/a^2 ] {\bf e_z}$, appelée force de scattering, pour une largeur de faisceau $a$ et d'amplitude $\varepsilon$. Nous avons réalisé des simulations numériques [4.1] en intégrant l'équation de Langevin sous-amortie [Eq. (4.1)], permettant d'obtenir l'état stationnaire des vortex browniens pour toute amplitude $\varepsilon$ et tout facteur de qualité $Q= \sqrt{\kappa m}/\gamma$.
Tout d'abord, nous avons dérivé la densité de probabilité stationnaire dans le régime sous-amorti en considérant la force de scattering comme une perturbation ($\varepsilon \ll 1$), et donc une déviation à la distribution de Maxwell-Boltzmann $P_0({\bf x},{\bf v})$ dans l'espace des phases $({\bf x},{\bf v})$. En intégrant sur les vitesses ${\bf v}$, nous avons alors obtenu une expression pour la densité de probabilité marginale $P_0({\bf x})$ et le courant ${\bf J} ({\bf x})$ formant le vortex brownien (voir Fig. 3a). La géométrie de ce vortex est indépendante du taux d'amortissement caractérisé par $Q$, tandis que son amplitude, proportionnelle à la circulation du courant, présente un maximum en fonction du facteur de qualité. Le courant dans l'espace des vitesses se présente également sous la forme d'un vortex dont seule l'amplitude dépend de $Q$, avec une évolution strictement croissante.
Ensuite, nous avons dérivé l'expression exacte de la densité spectrale $S_{zz}(\omega)$ du processus $z(t)$ correspondant à la transformée de Fourier de la corrélation $\langle z(t) z(0) \rangle$. Le processus $z$ peut être décomposé en deux processus indépendants : $z_{\rm B}$ soumis uniquement au bruit thermique et $z_{\rm nc}$ soumis à la force non conservative faisant intervenir les carrés des processus $x$ et $y$, dont les densités spectrales ont une forme lorentzienne avec un maximum en $\omega=\omega_0 = \sqrt{\kappa/m}$ et une largeur $\gamma/m$ à demi-hauteur. La densité spectrale $S_{zz}(\omega)$ présente deux pics pour $\omega= \sqrt \eta \omega_0$ et $\omega = 2 \omega_0$, correspondant à la signature respective des processus conservatif et non-conservatif. La force non-conservative influe également sur l'état stationnaire via le plateau situé en $\omega=0$, se comportant comme $1/\gamma$ dans la limite sous-amortie.
Pour finir, nous avons dérivé une expression exacte de la densité de probabilité stationnaire $P_s({\bf x})$ dans le régime suramorti, qui n'avait été calculé que dans la limite $\varepsilon \ll 1$ [4.5]. Les processus d'Ornstein-Uhlenbeck $x$, $y$ et $z_{\rm B}$ sont décrits par la distribution de Boltzmann dans l'état stationnaire. Seul le processus $z_{\rm nc}$ est hors-équilibre. Sa densité de probabilité $P_{\rm nc}$ peut être écrite sous la forme d'une intégrale de chemin sur les trajectoires gaussiennes $x(t)$ et $y(t)$, se réduisant dans l'espace de Laplace au calcul du déterminant d'un opérateur [4.6,4.7] et admettant une forme gaussienne en $x$ et $y$. Dans l'espace réel, la densité de probabilité stationnaire $P_s({\bf x})$ s'écrit comme la convolution de $\widehat P_{\rm nc}$ avec $\widehat P_{\rm B}$, la distribution du processus $z_{\rm B}$ dans l'espace de Laplace. Cette expression exacte nous permet de retrouver la théorie de perturbation d'ordre $\varepsilon$ et d'obtenir l'ordre suivant, autant pour la densité de probabilité que pour le courant (voir Fig. 3b).
Nous avons donc caractérisé les vortex browniens dans le régime sous-amorti pour des faibles amplitudes et dans le régime suramorti pour tout valeur d'amplitude, en dérivant une expression analytique pour la densité de probabilité, le courant et la densité spectrale en accord avec les observations expérimentales. Seul le régime sous-amorti pour des grandes amplitudes n'a pas été complètement analysé.
Ce travail a été effectué au cours de mon post-doctorat dans le groupe de H. Rieger à Sarrebruck (Allemagne), commencé en Novembre 2018. Il a donné lieu à deux publications [P11,P12].
Au cours de mon post-doctorat, j'ai analysé le temps mis par une particule brownienne pour s'échapper d'une cellule vivante, en considérant un modèle simplifié de l'organisation de son cytosquelette. Le transport intracellulaire de diverses particules-cargos (protéines, vésicules, mitochondries par exemple) vers la synapse immunologique située sur la membrane cellulaire, est crucial pour le bon fonctionnent des cellules et des organismes. Le cytosquelette est auto-organisé par un réseau de filaments, majoritairement constitué par des microtubules polarisés et des filaments d'actine, transportant rapidement les particules entre les différentes régions de la cellule. Ces deux types de filaments séparent le cytosquelette en deux régions : les microtubules émergeant radialement du centrosome, et les filaments d'actine distribués aléatoirement formant un fin cortex près de la membrane cytoplasmique [5.1]. Les particules-cargos réalisent une recherche intermittente [5.2,5.3], en alternant aléatoirement entre un transport diffusif dans le cytoplasme et un transport balistique sur les filaments [5.4]. Ce transport balistique est assisté par des moteurs moléculaires liant les particules-cargos aux filaments: la dynéine (du centrosome au cortex) et la kinésine (du cortex au centrosome) sur les microtubules, et la myosine sur les filaments d'actine. Il a été montré récemment que le temps pour atteindre la synapse via cette recherche intermittente pouvait être minimisé en jouant sur la taille du cortex, grâce à une diffusion effective plus importante dans celui-ci [5.5,5.6].
Nous nous sommes alors intéressés au calcul du temps moyen de premier passage (MFPT, pour mean first passage time) pour atteindre la synapse immunologique en considérant une géométrie simplifiée du cytosquelette (voir Fig. 4a), pour des particules passives. Ces particules réalisent un mouvement diffusif dans un disque unitaire composé de deux régions concentriques, avec une diffusivité $D_2$ dans la région centrale et $D_1 \ge D_2$ dans le cortex de taille $\Delta$, cohérent avec les études récentes [5.5,5.6]. De plus, une différence de potentiel $\Delta V = V_2 - V_1$ est présente entre ces deux régions, et la synapse est modélisée par un arc de cercle absorbant de taille $\varepsilon$. Le MFPT satisfait une équation de Poisson $D_i \nabla^2 t_i = -1$ dans chacune des deux régions avec une condition aux limites mixte : $t=0$ sur la synapse absorbante et $\partial_n t = 0$ sur le reste de la membrane réfléchissante [5.7]. Dans le cas d'un disque homogène ($D_1=D_2$ et $\Delta V=0$), une expression du MFPT a été dérivée dans la limite de petite ouverture $\varepsilon \ll 1$ [5.8] ainsi qu'une expression exacte [5.9] démontrée à partir de transformations conformes.
Nous avons principalement étudié le MFPT global (GMFPT), pour lequel la position de départ est moyennée par rapport à la densité de probabilité stationnaire, satisfaisant la distribution de Boltzmann. Les équations aux dérivées partielles du MFPT peuvent être résolues numériquement grâce à la méthode des éléments finis, avec FreeFem++ [5.10], donnant accès à une valeur numérique du GMFPT, équivalente aux simulations numériques réalisées avec la méthode de Monte Carlo cinétique [5.11,5.12]. Nous avons dérivé une expression analytique du MFPT et du GMFPT dans la limite $\varepsilon \ll 1$ en supposant la petite ouverture comme une perturbation [5.13,5.14,5.15], ainsi que dans la limite d'un fin cortex ($\Delta \ll 1$), similaire au problème de diffusion médiée par la surface [5.16].
Lorsque les particules ne peuvent pas entrer dans la région centrale, pour une interaction de cœur dur ($\Delta V= +\infty$), le GMFPT peut être minimisé pour un cortex de taille $\Delta < 1$ et ce minimum est parfaitement décrit dans la limite $\varepsilon \ll 1$. Sans différence de potentiel ($\Delta V=0$), le GMFPT est une fonction strictement décroissante de l'angle $\varepsilon$ et de la taille du cortex $\Delta$ et strictement croissante du ratio des diffusivités $D_1/D_2$. Aucune optimisation ne peut alors être réalisée. Cela peut s'expliquer par l'absence de biais forçant les particules à aller vers le cortex, comme en présence du cœur répulsif, d'une diffusion médiée par la surface [5.16] ou d'un transport balistique intermittent [5.5,5.6]. Malgré cela, nous avons pu montrer que le point de départ pour lequel le MFPT est maximal évolue discontinûment en raison d'une diffusion hétérogène (ralentissant les particules dans la région centrale), allant du point opposé à l'ouverture (situé dans le cortex) pour $D_1 = D_2$ au centre dans la limite $D_1 \gg D_2$ (voir Fig. 4b). Dans le cas général, nous avons montré que GMPFT présente un minimum pour $\Delta < 1$ si et seulement si la différence de potentiel $\Delta V$ est suffisamment grande, avec une valeur limite dépendant du ratio des diffusivité $D_1/D_2$ et de l'ouverture $\varepsilon$ (voir Fig. 4c). Le GMFPT tend vers une constante indépendante de la différence de potentiel dans la limite $\Delta \ll 1$, tandis que la limite $\varepsilon \ll 1$ capture parfaitement le minimum quand il existe.
Nous avons donc dérivé des expressions analytiques pour le MFPT pour une diffusivité hétérogène dans la limite des petites ouvertures, et montré qu'il était nécessaire d'avoir un mécanisme forçant les particules à aller vers la membrane cellulaire pour minimiser le MFPT avec la taille du cortex d'actine.
Ce travail a été effectué au cours de mon post-doctorat dans le groupe de H. Rieger à Sarrebruck (Allemagne), commencé en Novembre 2018. Il a donné lieu à quatre publications [P13,P14,P15,P16,P17,P18].
Pendant mon post-doctorat, en collaboration avec R. Paul, S. Chatterjee, et M. Karmakar de l'IACS de Calcutta (Inde), et J. D. Noh et C.-U. Woo de l'Université de Séoul (Corée du Sud), j'ai également étudié le mouvement collectif de particules actives s'alignant localement. Les mouvements collectifs sont largement observés dans la nature (nuées d'oiseaux, troupeaux de mammifères, bancs de poissons, essaims de bactéries) et également étudiés dans divers systèmes artificiels (gouttes de liquides, colloïdes roulants, disques polarisés vibrants) [6.1]. Ils sont constitués d'un grand nombre de particules actives qui consomment de l'énergie pour s'autopropulser et interagissent en s'alignant entre elles. Un mouvement synchrone de grands amas de particules, appelés nuées, émerge spontanément pour des grandes densités ou des faibles fluctuations thermiques. Cette transition vers les nuées (flocking transition en anglais) est un phénomène hors-équilibre décrit par de nombreux modèles théoriques. Le premier d'entre eux est le modèle de Vicsek (VM, pour Vicsek model) [6.2,6.3] pour lequel les particules actives réalisent un mouvement balistique avec une vitesse constante dans une direction s'alignant avec celles de leurs voisines. Cet alignement est similaire à une interaction ferromagnétique dans un bain thermique, et un mouvement collectif émerge grâce à la brisure spontanée de la symétrie de rotation. Une théorie hydrodynamique [6.4] a relié le VM à la classe d'universalité du modèle XY. Récemment, une version discrétisée sur un réseau bidimensionnel a été étudiée, le modèle d'Ising actif (AIM, pour active Ising model) [6.5,6.6], pour lequel les particules actives sautent vers la gauche ou la droite avec un biais, et s'alignent localement comme des spins d'Ising. La transition vers un mouvement collectif est très similaire au VM, par contre ce modèle permet d'obtenir de nombreux résultats analytiques.
Dans un premier temps, nous avons considéré une généralisation pour le modèle de Potts actif (APM, pour active Potts model) à $q$ états sur réseau, pour lequel les particules ont $q$ états internes décrits par le spin $\sigma_k$ correspondant à leurs directions de mouvement, et s'alignent localement avec une interaction de Potts défini avec l'hamiltonien $${\cal H}_i = -\frac{1}{2\rho_i} \sum_{k \ne l} (q\delta_{\sigma_k,\sigma_l}-1), \tag{6.1}$$ où $(k,l)$ désignent les particules sur le site $i$ et $\rho_i$ leur nombre. Dans ce résumé, je présenterai principalement les résultats pour $q=4$ sur un réseau carré, reproductibles pour $q=6$ sur un réseau triangulaire, en considérant $N$ particules sur un réseau de $L^2$ sites, pour une densité moyenne $\rho_0=N/L^2$. Pendant un intervalle de temps $\Delta t$, chaque particule dans un état de spin $\sigma$ peut soit changer de spin vers l'état $\sigma'$ avec une probabilité $\exp(- \beta \Delta {\cal H}_i) \Delta t$, pour une température $T=\beta^{-1}$, soit sauter vers un site voisin avec une probabilité $D(1-\epsilon/3)\Delta t$ dans une direction aléatoire ou $(4D\epsilon/3)\Delta t$ dans la direction $\sigma$. La probabilité totale de saut $4D\Delta t$ est indépendante du biais $\epsilon$. Des simulations numériques sont alors réalisées en utilisant l'algorithme de Metropolis [6.7], optimisé pour $\Delta t = [4D + \exp(4\beta)]^{-1}$.
Dans la limite hydrodynamique ($L\gg1$), nous avons dérivé à partir du modèle microscopique les équations pour la densité de particules dans l'état $\sigma$, notée $\rho_\sigma({\bf x},t)$, s'écrivant sous la forme : $$\partial_t \rho_\sigma = D_\parallel \partial_\parallel^2 \rho_\sigma + D_\perp \partial_\perp^2 \rho_\sigma - v \partial_\parallel \rho_\sigma + \sum_{\sigma' \ne \sigma } K_{\sigma \sigma'}, \tag{6.2}$$ avec une diffusivité longitudinale $D_\parallel = D(1+\epsilon/3)$, une diffusivité transverse $D_\perp = D(1-\epsilon/3)$ et une vitesse d'autopropulsion $v=4D\epsilon/3$. Cette équation peut être résolue par la méthode des éléments finis avec FreeFem++ [6.8] en discrétisant le temps $t$.
Quand le biais est nul ($\epsilon=0$), une transition de phase du premier ordre a lieu pour une densité $\rho_*(T)$ causée par la brisure spontanée de la symétrie. Elle sépare une phase désordonnée à basse densité, pour laquelle tous les états sont identiquement peuplés avec une magnétisation $m_0 = 0$, d'une phase ordonnée à haute densité, pour laquelle un état est surpeuplé par rapport aux trois autres avec une magnétisation $m_0 > 0$. Pour une température supérieure à $T_c \simeq 2.4$, seule la phase désordonnée existe.
Pour un biais non nul ($\epsilon>0$), nous observons une séparation de phase sous forme de bandes et d'allées, respectivement pour des petits et grands biais (voir Fig. 5a). La phase ordonnée (liquide) se déplace sur la phase désordonnée (gaz) immobile avec un mouvement transverse pour les bandes et longitudinal pour les allées. La réorientation d'une bande en une allée doit son origine à la diminution de la diffusion transverse $D_\perp$ et à l'augmentation de la vitesse $v$ quand le biais devient grand, stabilisant les allées. Bien que cette réorientation est absente du VM et de l'AIM, la transition vers le mouvement collectif reste similaire à une transition de phase liquide-gaz sans région supercritique, définissant les binodaux $\rho_g$ et $\rho_l$ des phases gazeuse et liquide, permettant de construire les diagrammes de phase biais-densité (Fig. 5a) et température-densité (Fig. 5b). L'étude de la stabilité des solutions homogènes de l'Eq. (6.2) ($K_{\sigma \sigma'}=0$), obtenue en y ajoutant une petite perturbation, permet de dériver une expression analytique du biais pour lequel la réorientation a lieu.
Nous avons ensuite ajouté une restriction sur le nombre de particules $\rho_i$ pouvant occupé un site via deux formalismes : (i) en fixant un nombre maximal de particules MPS par site en empêchant le saut d'une nouvelle particule si il est atteint, ou (ii) en ajoutant un potentiel d'interaction $U \rho_i(\rho_i-1)$, revenant à multiplier la probabilité de saut par $\kappa(\rho_i) = \exp(-2 \beta U \rho_i)$. L'équation hydrodynamique [Eq. (6.2)] est légèrement modifiée avec un courant maintenant égal à $J_\parallel = D_\parallel[\kappa(\rho)\partial_\parallel \rho_\sigma - \kappa'(\rho)\rho_\sigma \partial_\parallel \rho] + v \kappa(\rho)\rho_\sigma$ dans la direction longitudinale et $J_\perp = D_\perp[\kappa(\rho)\partial_\perp \rho_\sigma - \kappa'(\rho)\rho_\sigma \partial_\perp \rho]$ dans la direction transverse, avec $\kappa(\rho_i) = 1 - \rho_i/{\rm MPS}$ dans le cas (i).
Avec une faible restrictions, les résultats discutés précédemment restent valables, cependant lorsque les restrictions deviennent importantes (MPS petit ou $U$ grand), le mouvement collectif en bandes et en allées est arrêté, formant un nouvel état nommé bloqué et présent pour des basses températures $T$, des grands biais $\epsilon$, ou des densités $\rho_0$ élevées. Ce phénomène avait déjà été observé pour des particules actives avec interactions, pour lesquels une séparation de phase induite par la motilité où des zones densément peuplées coexistent avec des zones moins denses [6.9,6.10,6.11]. Nous avons alors construit les diagrammes de phases modifiés en présence de ces restrictions (Fig. 5c). Nous avons également réalisé une étude de la stabilité des solutions homogènes de l'équation hydrodynamique modifiées, montrant que le biais pour lequel la réorientation a lieu augmente avec la restriction (quand MPS diminue, ou quand $U$ augmente). Cependant, l'existence de l'état bloqué ne peut pas être comprise avec cette étude de stabilité.
Dans un second temps, nous avons considéré le modèle d'horloge à $q$ états hors réseau comme une discrétisation naturelle du VM, pour lequel les particules ont $q$ états internes décrits par l'angle $\theta_k$, où l'interaction de Potts donnée par l'Eq. (6.1) est remplacée par $${\cal H}_i = -\frac{1}{2\rho_i} \sum_{k \ne l} \cos(\theta_k-\theta_l), \tag{6.3}$$ où $(k,l)$ représentent les particules dans le voisinage de la particule $i$. La limite $q\to \infty$, correspondant au modèle XY actif, appartient à la même classe d'universalité quel le VM. Pendant un intervalle de temps $\Delta t$, chaque particule dans un état $\theta$ peut soit changer d'état vers $\theta'$ avec une probabilité $\exp(- \beta \Delta {\cal H}_i) \Delta t$, soit se déplacer dans une direction aléatoire avec une probabilité $D(1-\epsilon)\Delta t$ ou dans la direction $\theta$ avec une probabilité $D\epsilon \Delta t$. Des simulations numériques sont alors réalisées en utilisant l'algorithme de Metropolis [6.7], optimisé pour $\Delta t = [D + \exp(2\beta)]^{-1}$. Un modèle légèrement différent a été étudié dans [6.12], et une autre discrétisation du VM a été proposée dans [6.13] avec des résultats qualitativement similaires.
Pour un petit nombre de directions ($q \lesssim 5$), nous avons constaté que la transition vers le mouvement collectif de l'ACM a le même comportement que celui de l'APM présenté précédemment, avec une réorientation de bandes en allées pour des grands biais $\epsilon$. La séparation de phase est macroscopique (une seule bande est observée dans l'état stationnaire) et la phase liquide conserve une orientation fixe au cours du temps pour une direction $\theta$ donnée. Ce comportement peut être expliqué par des fluctuations du nombre de particules $\Delta n^2 = \langle n^2 \rangle - \langle n \rangle^2$ normales, i.e. linéaire avec le nombre de particules: $\Delta n^2 \sim \langle n \rangle$. En revanche, pour un plus grand nombre de directions ($q \gtrsim 8$), la transition vers le mouvement collectif devient équivalent à celui du VM, sans réorientation pour des grands biais $\epsilon$. La séparation de phase est microscopique (plusieurs bandes peuvent être observées dans l'état stationnaire avec une taille maximale pour une bande) et la phase liquide ne conserve pas une orientation fixe au cours du temps. Cela peut être expliqué par des fluctuations du nombre de particules gigantesques où $\Delta n^2 \sim \langle n \rangle^{1.67}$.
Nous avons donc décrit le mouvement collectif de spins avec une interaction ferromagnétique locale. La transition vers celui-ci est similaire à une transition de phase de type liquide-gaz sans région supercritique, comme pour l'AIM et le VM. De plus, nous avons observé une réorientation des bandes en allées, absente des modèles précédents, et expliqué son origine. En ajoutant une restriction sur le nombre de particules sur un site, nous avons obtenu un état bloqué, qui avait déjà été observé dans le contexte d'autres systèmes actifs. Finalement, nous avons réussi à relié ces résultats avec le VM pour un modèle d'horloge actif à $q$ états.
Les systèmes complexes sont typiquement hétérogènes, car les individus varient dans leurs propriétés, leur réponse à l'environnement externe et leurs interactions mutuelles. En particulier, de nombreux systèmes biologiques possédant un mouvement collectif peuvent être modélisés par des particules autopropulsées avec des interactions hétérogènes, ce qui motive l'étude de populations avec plusieurs espèces.
Dans un premier temps, nous avons considéré la variante à deux espèces du VM (TSVM, pour two-species Vicsek model), constituée de deux types de particules autopropulsées qui tendent à s'aligner avec les particules de leur propre espèce (interaction ferromagnétique) et à s'anti-aligner avec celles de l'autre espèce (interaction anti-ferromagnétique). Ce système peut être modélisé par $N$ particules dans un domaine de taille $L_x \times L_y$, dotées chacun d'une position ${\bf r}_i^t$, d'une orientation $\theta_i^t$ et d'un spin $s_i^t = \pm 1$ caractérisant l'espèce ($s_i^t = +1$ pour l'espèce A et $s_i^t = -1$ pour l'espèce B), avec autant de particules dans chacune des deux espèces. L'évolution temporelle suit des équations de Langevin similaire au VM : $\theta_i^{t+1} = \langle \theta_i^t \rangle + \eta \xi_i^t$ et ${\bf r}_i^{t+1} = {\bf r}_i^t + {\bf e}(\theta_i^{t+1})$, où $\xi_i^t$ est un bruit blanc uniforme dans $[-\pi,\pi]$ et ${\bf e}(\theta)=(\cos \theta, \sin \theta)$. Le calcul de l'angle moyen ressenti par la particule $i$ dans son voisinage ${\cal N}_i$ est modifié pour correspondre à la dynamique des deux espèces tel que $$\langle {\bf e}(\theta_i^t) \rangle = \sum_{j \in {\cal N}_i} s_i^t s_j^t {\bf e}(\theta_j^t). \tag{6.4}$$ Ces équations de Langevin couplées peuvent être simulées numériquement avec un algorithme de dynamique moléculaire. Le mouvement collectif est analogue au modèle original à une espèce avec une transition de phase liquide-gaz sans région supercritique, et présente une séparation de phase microscopique où plusieurs bandes de chaque espèce se déplacent sur une phase gazeuse immobile.
Cependant, deux états dynamiques apparaissent dans la région de coexistence : l'état PF (parallel flocking), où toutes les bandes des deux espèces se déplacent dans la même direction, et l'état APF (anti-parallel flocking), où les bandes des deux espèces se déplacent dans des directions opposées (voir Fig. 6a). L'état PF n'existe que pour une densité inférieure à $\rho_{\bf PF}$ où les bandes de chaque espèce ne se recouvrent pas, tandis que la phase liquide est dans l'état APF. Nous avons alors construit le diagramme de phase bruit-densité correspondant au TSVM (voir Fig. 6a). Les états PF et APF effectuent des transitions stochastiques de l'un à l'autre, dans la région PF+APF du diagramme de phase (partie faible densité, grand bruit de la région de coexistence). La fréquence de ces transition et les temps de persistance de chaque état dépend drastiquement du rapport entre la largeur des bandes et la taille longitudinale du système : les transitions sont beaucoup plus fréquentes quand $L_x \lesssim 200$.
Dans un second temps, nous avons considéré la variante à deux espèces de l'AIM (TSAIM, pour two-species active Ising model). Ce système peut être modélisé par $N$ particules dans un réseau périodique de taille $L_x \times L_y$, dotées chacun d'un spin $\sigma_i = \pm 1$ caractérisant l'orientation ($\sigma_i = +1$ vers la droite et $\sigma_i = -1$ vers la gauche) et d'un spin $s_i = \pm 1$ caractérisant l'espèce ($s_i = +1$ pour l'espèce A et $s_i = -1$ pour l'espèce B). Pendant un intervalle de temps $\Delta t$, chaque particule dans un état $(\sigma,s)$ peut soit changer d'état vers $(-\sigma,s)$ avec une probabilité $p_1$ ou vers $(\sigma,-s)$ avec une probabilité $p_2$, définies par $$p_1 = \gamma_1 \exp(- 2\beta_1 J_{ss'} m_{s'}) \Delta t, \qquad p_2 = \gamma_2 \exp(-2\beta_2 \mu) \Delta t, \tag{6.5}$$ où $J_{ss'}$ est le couplage entre les espèces, $m_s = \langle \sigma_k \rangle_{k \in s}$ est la magnétisation de l'espèce $s$, et $\mu = \langle s_k \rangle_k$ est la magnétisation entre les espèces, soit se déplacer dans la direction ${\bf p}$ avec une probabilité $D(1+\epsilon\sigma {\bf p} \cdot {\bf e_x})\Delta t$. Deux températures effectives peuvent être considérées : $T_1 = \beta^{-1}$ et $T_2 = \beta_2^{-1}$. Des simulations numériques sont alors réalisées en utilisant l'algorithme de Metropolis [6.7], optimisé pour $\Delta t = [4D + \gamma_1 \exp(2\beta_1) + \gamma_2 \exp(2\beta_2)]^{-1}$. Nous avons considéré trois sous-modèles : (i) sans changement d'espèce ($\gamma_2 = 0$) et avec une interaction réciproque entre les espèces pour $J_{AA} = J_{BB} = 1$ et $J_{AB} = J_{BA} = -1$, revenant à la version discrète du TSVM discuté précédemment; (ii) toujours avec interaction réciproque mais avec un changement d'espèce ($\gamma_2 = 0.5$), correspondant à une extension active du modèle d'Ashkin-Teller [6.14]; et (iii) sans changement d'espèce ($\gamma_2 = 0$) mais avec une interaction non-réciproque entre les espèces : $0 \le J_{AB} = - J_{BA} = J_{\rm NR} \le 1$, signifiant que les particules de l'espèce A tendent à s'aligner avec celles de B (interaction ferromagnétique) tandis que celles de B tendent à s'anti-aligner avec celles de A (interaction anti-ferromagnétique). Pour chacun de ces sous-modèles, nous avons dérivés une théorie hydrodynamique reproduisant parfaitement tous les résultats mentionnés ci-après.
Tout d'abord, pour le sous-modèle (i), le mouvement collectif est similaire au modèle original à une espèce avec une transition de phase liquide-gaz sans région supercritique, et présente une séparation de phase macroscopique où une bande de chaque espèce se déplace sur la phase gazeuse immobile, sous la forme d'états PF et APF, comme pour le TSVM (voir Fig. 6b). Cependant, l'état PF existe pour toute densité et est plus stable que l'état APF, en partant d'une configuration désordonnée dans la limite hydrodynamique ($L_x,L_y \gg 1$), autant dans la région de coexistence que pour la phase liquide. L'état PF à haute densité, observé pour la phase liquide, se compose d'une bande de chaque espèce se déplaçant dans la même direction et remplissant tout l'espace. De plus, contrairement au TSVM, aucune transition stochastique ne s'effectue entre les états PF et APF, suggérant que les fluctuations du nombre de particules gigantesques du VM en sont responsables. Nous avons alors construit les diagrammes de phase biais-densité et température-densité (voir Fig. 6b).
Pour le sous-modèle (ii), la nature de la transition vers le mouvement collectif est plus riche, en raison d'une transition de phase supplémentaire pour le changement d'espèce, pour une densité dépendant uniquement de $T_2$. Lorsque $T_2 \lesssim T_1$, cette transition de phase se produit pour la phase gazeuse séparant un gaz paramagnétique, où l'orientation et l'espèce sont désordonnées, d'un gaz ferromagnétique, où l'orientation est désordonnée et l'espèce ordonnée. Le mouvement collectif présente une séparation de phase macroscopique sous deux formes stables : (a) une seule espèce survit avec une bande se déplaçant sur la phase gazeuse, et (b) les deux espèces survivent dans un état PF constitué d'une bande de chaque espèce suivie par une phase gazeuse de la même espèce. De même, l'état liquide est soit composé d'une seule espèce, soit de deux espèces sous la forme d'un état PF à haute densité, déjà observé pour le sous-modèle (i). Lorsque $T_2 \gg T_1$, la transition de phase du changement d'espèce se produit pour une densité plus grande que la transition vers le mouvement collectif. Le mouvement collectif présente une séparation de phase microscopique où plusieurs bandes de chaque espèce se déplacent dans un état PF sur une phase gazeux paramagnétique, causé par la frustration d'un état avec une orientation ordonnée mais une espèce désordonnée. La phase liquide est, quand à elle, sous la forme d'un état PF à haute densité dont le nombre de bandes dépend de la condition initiale. À partir de ces différentes phases, nous avons alors construit les diagrammes de phase température $T_1$-densité, température $T_2$-densité et biais-densité.
Pour le sous-modèle (iii), la nature de la transition vers le mouvement collectif peut être étudié en augmentant le paramètre $J_{\rm NR}$. Quand $J_{\rm NR}=0$, les deux espèces n'ont aucune interactions entre elles, et suivent chacune d'entre elles les propriétés de l'AIM. Quand $J_{\rm NR}$ est augmenté, l'interaction entre les deux espèces est principalement défavorable à l'espèce B, provocant une disparition de son mouvement collectif. Pour $J_{\rm NR} \sim 0.5$, il est possible que plus aucun mouvement collectif ne soit constaté (état gazeux), dû à l'interaction non réciproque. Cependant, quand $J_{\rm NR}$ est suffisamment grand, le mouvement collectif apparaît à nouveau pour les deux espèces où l'espèce B suit l'espèce A en formant un état course-et-poursuite (ou run-and-chase), non observé dans les études récentes de systèmes non-réciproques [6.15,6.16] (voir Fig. 6c). Cet état se forme lorsque les particules de l'espèce A viennent s'accumuler à l'avant de la bande d'espèce B, provocant un ralentissement global de la dynamique. Nous avons alors construit les diagrammes de phase couplage-température et couplage-densité (voir Fig. 6c). Quand le biais est nul ($\epsilon=0$), trois états peuvent être observés : un état désordonné dans la limite des basses densités et hautes températures, un état oscillant où la magnétisation de l'espèce A oscille en quadrature de phase avec celle de l'espèce B, et un état ordonné dans la limite des hautes densités et basses températures. Ce résultat vient d'être également observé pour le modèle d'Ising non-réciproque [6.17].
Récemment, des études ont montré qu'un petit obstacle ou une gouttelette artificiellement excitée peuvent déstabiliser la phase liquide du VM et de l'AIM [6.18,6.19]. Pour ce dernier, une fluctuation émergeant spontanément peut déstabiliser la phase liquide [6.20]. Cette instabilité est également présente pour la phase liquide du TSAIM quand la diffusion est faible par rapport à la vitesse de propulsion des particules. Pour le sous-modèle (i), cette instabilité conduit à un nouvel état stable, également sous la forme d'un état PF mais avec plusieurs bandes, où l'orientation a un ordre à longue portée tandis que l'espèce a un ordre à courte portée. Pour le sous-modèle (ii), cette instabilité conduit à un nouvel état stable similaire à l'AIM [6.20], où l'orientation a un ordre à courte portée tandis que l'espèce a un ordre à longue portée. Pour le sous-modèle (iii), cette instabilité conduit à un état désordonné, en raison de l'interaction non-réciproque. De plus, nous avons observé un état bloqué (MIIP pour motility induced interface pinning) dans la limite de faible températures, similairement à l'AIM [6.20], quel que soit le sous-modèle étudié. Sa structure est similaire à celui discuté pour l'APM.
Nous avons donc décrit le mouvement collectif pour des modèles actifs à deux espèces en présence d'interactions ferromagnétiques et anti-ferromagnétiques. La transition vers celui-ci est similaire à une transition de phase de type liquide-gaz sans région supercritique, comme pour l'AIM et le VM, en présence d'interactions réciproques. De plus, nous avons observé la présence de deux états dynamiques : PF et APF, ayant des comportements totalement différents suivant le modèle étudié (état liquide APF et transition entre états pour le TSVM, état liquide PF et aucune transition pour le TSAIM). Finalement, nous avons également étudié l'impact du changement d'espèce et de la présence d'interactions non-réciproques, ainsi que la stabilité de la phase liquide pour le TSAIM.
Ce travail a été effectué au cours de mon post-doctorat dans le groupe de H. Rieger à Sarrebruck (Allemagne), commencé en Novembre 2018. Il a donné lieu à une publication [P19].
Au cours de mon post-doctorat, j'ai étudié les courants stationnaires produits par la matière active sous l'effet de la gravité en présence d'une paroi, en considérant un modèle simplifié de particules browniennes actives. La matière active est constituée de particules auto-propulsées, comme les microorganismes, bactéries ou colloïdes actifs, qui consomment de l'énergie à petite échelle et génèrent un mouvement persistant, maintenant le système hors équilibre [7.1]. Les états stationnaires sont non triviaux en présence de parois ou d'obstacles, avec des zones d'accumulation de particules [7.2]. Un modèle minimal fréquemment étudié est celui des particules browniennes actives (ABP, pour active Brownian particles) [7.3], montrant une séparation de phase induite par la motilité (MIPS, pour motility induced phase separation), séparant une phase dense et lente d'une phase diluée et rapide [7.4,7.5]. Ce mécanisme est également responsable de l'accumulation de particules actives sur les parois répulsives. Bien que ce phénomène soit bien compris dans le cas des ABP, l'effet des bords et des forces d'interaction stériques sous un champ externe (e.g. la gravité) reste mal exploré. Des expériences récentes ont montré que les ABP peuvent sédimenter sous gravité [7.6,7.7,7.8] ou même remonter les parois grâce à leur activité, illustrant des comportements analogues à l'action capillaire, mais sans forces attractives [7.9,7.10]. Enfin, des courants stationnaires émergent souvent, même en l'absence de mécanismes explicites d'alignement, révélant une organisation complexe et coopérative des particules actives [7.11,7.12].
Nous avons étudié un modèle d'ABP avec interactions et soumises à la gravité, en considérant $N$ particules circulaires de rayon $R$ dans un domaine de taille $L_x \times L_y$ avec des parois réfléchissantes. Les particules sont autopropulsées avec une vitesse $v_s$ constante dans la direction $\theta_i$, performant une diffusion rotationnelle, avec une diffusivité $D_r$. Les coordonnées $\{{\bf r}_i(t), \theta_i(t)\}$ des particules obéissent aux équations de Langevin : $$\dot{{\bf r}_i}=v_s{\bf \hat{e}}(\theta_i)-v_g{\bf \hat{y}}+\frac{{\bf F}_i}{\gamma}, \qquad \dot{\theta_i}= \sqrt{2D_r}\eta, \tag{7.1}$$ où $\gamma$ est le coefficient de friction, $\eta$ est un bruit thermique, considéré comme un bruit blanc gaussien, et ${\bf \hat{e}}(\theta_i)=(\cos \theta_i, \sin \theta_i)$. Le mouvement d'une particule est gouverné par la vitesse d'autopropulsion $v_s$, la vitesse de sédimentation $v_g$, et l'interaction ${\bf F}_i$ avec les particules voisines, modélisé par une force élastique de raideur $k$ dans un rayon inférieur à $2R$, et les parois du domaine. Le problème peut être étudié avec le nombre de Péclet ${\rm Pe} = v_s/aD_r$, le ratio entre les vitesses de sédimentation et d'autopropulsion $\alpha=v_g/v_s$, et l'intensité de la force de répulsion $F_0 = k a / \gamma v_s$. Les équations de Langevin couplées (7.1) peuvent être simulées numériquement avec un algorithme de dynamique moléculaire, donnant accès aux densités de particules $\rho({\bf r})$ et de courant ${\bf J}({\bf r})$ stationnaires. L'amplitude du courant est généralement donné par le rotationnel $A({\bf r}) = \partial_x J_y - \partial_y J_x$.
En présence d'interaction ($F_0 = 100$), les particules forment un ménisque sur les parois verticales (voir Fig. 7a) avec une élévation capillaire $\Delta h \sim \lambda_{\rm sed}^{1.8}$, où $\lambda_{\rm sed} = {\rm Pe}/\alpha$ est la longueur caractéristique de la sédimentation loin des parois. La présence de ce ménisque peut être expliquée par la présence d'une MIPS (entre une phase dense, presque cristalline, et une phase diluée) quand l'interaction entre les particules est forte. La formation de ce ménisque est déterminée par un courant stationnaire, formé par deux vortex, centrés à la base de celui-ci (voir Fig. 7b), avec une taille et une intensité augmentant avec l'activité des particules. L'origine de ces vortex peut être retracée à partir du confinement des particules : déjà, l'état stationnaire des ABP sans interaction et sans gravitation montre des courant qui s'organisent symétriquement dans les coins du domaine en huit vortex. La gravitation déforme cette configuration de vortex vers le bas, laissant deux vortex principaux sur les parois verticales. Les interactions répulsives entre les ABP ne modifient cette situation que lorsque la MIPS apparaît et forme une région dense au fond, ce qui pousse le centre des vortex vers le haut, en direction du ménisque.
En absence d'interaction ($F_0 = 0$), les équations de Langevin (7.1) ne sont plus couplées, permettant d'écrire une équation de Fokker-Planck pour la densité de probabilité $p({\bf r}, \theta; t)$ : $$\partial_t p = \nabla \cdot [ D_t \nabla p - (v_s {\bf e}_\theta - v_g {\bf \hat{y}} ) ] + D_r \partial_\theta^2 p, \tag{7.2}$$ où $D_t$ est la diffusivité translationnelle ajoutée pour stabiliser la dynamique. Cette équation peut être résolue numériquement grâce à la méthode des éléments finis, avec FreeFem++ [7.13], donnant accès à une valeur numérique pour les densités de particules $\rho({\bf r}) = \int d\theta p({\bf r}, \theta)$ et de courant ${\bf J}({\bf r}) = -D_t \nabla \rho({\bf r}) + v_s {\bf P}({\bf r}) - v_g \rho({\bf r}) {\bf \hat{y}}$ stationnaires, où $P({\bf r}) = \int d\theta {\bf e}_\theta p({\bf r}, \theta)$ est le vecteur polarisation. Même si la MIPS est absente, les particules s'accumulent proche des parois verticales (voir Fig. 7c), avec une élévation capillaire cependant moins forte : $\Delta h \sim \lambda_{\rm sed}^{1.1}$. Cette accumulation peut être encore expliquée par la présence d'un courant stationnaire, formé par trois vortex, situé dans les coins du domaine (voir Fig. 7d), avec une taille diminuant et une intensité augmentant avec l'activité des particules.
Nous avons donc décrit la sédimentation de particules actives faisant apparaître un courant stationnaire. Les particules s'accumulent proche des parois verticales avec une élévation capillaire dépendant de la longueur de sédimentation. Cette accumulation peut être expliquée par la présence d'un ou plusieurs vortex dont la géométrie et l'intensité dépend de l'activité des particules.