# Flots protostellaires poussés par un jet dans une enveloppe stratifiée: revisiter un modèle à l'ère d'ALMA

Mialy Rabetanety Rabenanahary

► To cite this version:

Mialy Rabetanety Rabenanahary. Flots protostellaires poussés par un jet dans une enveloppe stratifiée: revisiter un modèle à l'ère d'ALMA. Astrophysique stellaire et solaire [astro-ph.SR]. Observatoire de Paris, 2022. Français. NNT: tel-04122023

HAL Id: tel-04122023

https://hal.science/tel-04122023

Submitted on 8 Jun 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Flots protostellaires poussés par un jet dans une enveloppe stratifiée : revisiter un modèle à l'ère d'ALMA

Soutenue par

Mialy Rabetanety
RABENANAHARY

Le 06 décembre 2022

Ecole doctorale n° 127

AAIF - Astronomie et
Astrophysique d'Ile de
France

Spécialité

Astronomie et Astrophysique

## Composition du jury :

Martina Wiedner
Directeur de recherche
Observatoire de Paris

Présidente du jury

Bertrand Lefloch
Directeur de recherche
IPAG Grenoble

Rapporteur

Jean-Marc Huré
Professeur
LAB Université de Bordeaux

Rapporteur

Lars Egstrøm Kristensen
Associate professor
University of Copenhagen

Examinateur

Alejandro Raga
Professeur
UNAM Mexico

Examinateur

Sylvie Cabrit
Astronome
Observatoire de Paris

Directrice de thèse

Guillaume Pineau des Forêts
Professeur émérite
Université Paris-Saclay

Invité

Zakaria Meliani
Astronome adjoint
Observatoire de Paris

Invité
Remerciements

Une thèse est une aventure de tous les instants, avec tous ses moments compliqués et ses défis quotidiens. Je tiens tout d'abord à remercier mes encadrants et collaborateurs, sans qui cette thèse aurait été impossible. Merci beaucoup à Sylvie, Guillaume et Zakaria, pour leur bienveillance, leurs encouragements, leurs enseignements, et leur patience tout au long de ces trois années riches en moments forts. Merci beaucoup à Martina pour avoir été une super marraine au LERMA, pour ses nombreux conseils et son suivi, ainsi que pour sa gentillesse. Un grand merci à Antoine, qui m'a aidé à me lancer dans cette thèse,

Une thèse, c'est aussi trois ans dans un laboratoire (ou plusieurs, selon le jour de la semaine) avec des collègues. C'est pourquoi je tiens à remercier tous mes camarades doctorants pour leurs discussions passionnantes et les repas animés à la cantine. Merci tout particulièrement à Daniel, Aristide, Barbara, Lu, Pierre, Guillaume. Je remercie aussi Benoît, dont les nombreux passages au laboratoire ont toujours été l'occasion d'échanges intellectuels très riches sur les flots protostellaires et les disques d'accrétion, mais aussi de très bons conseils pour faire face aux aléas d'une telle aventure.

L'astrophysique offre la chance de rencontrer et d'interagir avec des chercheurs du monde entier. Je tiens ainsi à remercier tous les membres de la collaboration ALMA-FLOW pour cette incroyable opportunité. Même si notre premier proposal sur ALMA n'est pas passé, je remercie tout particulièrement Doug, Lars, Bertrand et Hsien pour le plaisir d'échanger, écouter, discuter, et expérimenter un tel travail.

Je remercie aussi le personnel administratif du LERMA et de l'Observatoire pour avoir permis le déroulement de cette thèse dans les meilleures conditions possibles.

Je remercie les rapporteurs, pour leur lecture attentive de ce manuscrit, et les membres du jury pour avoir pris le temps de venir à Paris assister à cette soutenance inoubliable. Je les remercie pour leurs félicitations, et pour ce magnifique dîner convivial. À cette occasion, je tiens à remercier Alex, Lars, Bertrand, Jean-Marc pour les nombreux conseils, l'écoute sur mes projets d'avenir, et pour m'avoir convaincu dans la poursuite de mes objectifs et de mes ambitions.

Merci à l'équipe informatique du LERMA et de l'Observatoire, notamment à Tarik, pour m'avoir aidé à réaliser des simulations numériques avec le meilleur matériel et des ressources informatiques aussi performantes. Merci à l'équipe d'encadrement de l'Ecole thématique AstroSim 2019 ayant eu lieu au CINES, pour ces deux semaines inoubliables au milieu de tous ces serveurs dernier cri,

Pour cette thèse, je suis reconnaissant du soutien financier apporté par le Programme National Physique et Chimie du Milieu Interstellaire (PCMI) du CNRS/INSU avec l'INC/INP, et co-financé par le CEA et le CNES.

Enfin, je tiens en particulier à remercier infiniment ma meilleure amie, Nirina, pour avoir embelli mon quotidien, pour m'avoir soutenu, encouragé, recadré tout au long de ces trois dernières années. Nos fous rires, nos discussions, nos délires, nos galères et nos découvertes, un grand merci à toi pour tous ces bons moments, et pour m'avoir aidé à tenir et à garder l'équilibre devant chaque défi de cette thèse.
Résumé

L'efficacité observée dans la formation des étoiles jusqu'au parsec ne dépasse pas quelques pour-
cents, sans oublier le décalage de la fonction initiale de masse (IMF) à seulement ~30% de la masse
du cœur pré-stellaire initial. Comprendre en détail le processus derrière d'aussi faibles efficacités
reste encore à ce jour une question ouverte. En outre, les simulations numériques les plus récentes
ont démontré que la turbulence et les champs magnétiques à eux seuls ne peuvent suffire à repro-
duire de telles valeurs. Elles montrent que la rétroaction des flots protostellaires joue un rôle crucial
en perturbant les écoulements d'accrétion, en évacuant la matière des cœurs, et/ou en maintenant
la turbulence. Malheureusement, que ce soit en termes de volume de nuage affecté, d'impulsion in-
jectée, de masse entraînée, ou d'impact sur le disque et l'enveloppe en effondrement. L'importance
de cette rétroaction dépend fortement de la géométrie sous-jacente du vent protostellaire. Cette
dernière reste encore débattue : "vent X grand angle" rapide, vent de disque MHD plus lent, ou
jet collimaté ? De toute évidence, afin d'évaluer fiablement l'impact de la rétroaction des flots sur
la formation stellaire, il est d'une importance cruciale de déterminer la géométrie de vent la plus
réaliste (et/ou les géométries que nous pouvons exclure). Pour apporter une nouvelle contribution
quant à cette question, nous présentons des simulations numériques de flots poussés par un jet pulsé
collimaté, lancé à travers un cœur pré-stellaire stratifié. Nous comparons nos simulations avec les
observations ALMA récentes, ainsi qu'avec les prédictions analogues pour un vent X grand angle.
Nos simulations sont les premières à combiner sur une échelle de 0.1 pc la variabilité du jet, la
stratification en densité de l'enveloppe et des échelles de temps de 10 000 ans comparables aux flots
jeunes observés. Les prédictions de nos simulations en termes de largeur de flot, de diagrammes
position-vitesse, et de distribution masse-vitesse, montrent une ressemblance frappante avec les ob-
servations ALMA de flots CO tels que HH46/47 et CARMA-7. L'accord est même plus prometteur
qu'avec les modèles de flots poussés par un "vent X grand angle". Ces résultats pourraient avoir
une implication majeure sur le rôle des flots dans la régulation de la formation stellaire.
<figure><img src="image_2.png" /></figure>
Abstract

A long-standing open question in star formation is the process responsible for its low efficiency on parsec scales (a few %), and for shifting down the Initial Mass Function (IMF) to only ~30% of the prestellar core mass distribution. The most recent numerical simulations show that neither turbulence nor magnetic fields can, alone, reproduce these low efficiencies, and that feedback by protostellar outflows must play a crucial role by disrupting accretion streams, expelling material from cores, and/or sustaining turbulence. Unfortunately, the magnitude of outflow feedback (affected cloud volume, injected momentum, entrained mass, impact on the disk and infalling envelope) depends strongly on the underlying protostellar wind geometry, which remains uncertain and heavily debated: a fast wide-angle "X-wind", a slower MHD disk wind, a narrow jet? Clearly, if we want to reliably assess the role of outflow feedback in star formation, it is of utmost importance to determine which wind geometry is the most realistic (and/or which one can be excluded). As a new contribution towards this goal, we present, for the first time, numerical predictions for outflows driven by a narrow pulsed jet in a stratified prestellar core. We compare our simulations against recent ALMA observations and analogous predictions for a wide-angle X-wind. Our simulations are the first to combine jet variability, ambient density-stratification, and long timescales up to 10 000 yrs (typical of young outflows) on scales up to 0.1 pc. We find that the predicted widths, position-velocity diagrams, and mass-velocity distribution, show striking resemblance with ALMA observations of CO outflows such as HH46/47 and CARMA-7, and in closer agreement than models based on a wide-angle "X-wind". The results obtained in this work could have major implications for the feedback of protostellar outflows on star formation.
<figure><img src="image_3.png" /></figure>
# TABLE DES MATIÈRES

<table><tr><td>Remerciements</td><td>2</td></tr></table>

## I Introduction 1

### 1 La formation d'une étoile de type solaire 1

1.1 Phase pré-stellaire 1

1.2 Phase protostellaire 2

1.2.1 Classe 0 2

1.2.2 Classe 1 4

1.2.3 Classe 2 4

1.2.4 Classe 3 5

### 2 Rôle des éjections protostellaires et questions ouvertes 5

2.1 Impact sur la formation des étoiles et sur le nuage moléculaire parent 5

2.2 Origine de l'éjection et moment cinétique 7

2.3 Impact sur la formation planétaire 7

2.4 Deux prescriptions de vent pour la rétroaction du flot protostellaire 7

2.5 Retour d'un scénario et objectifs de la thèse 9

2.6 Plan de la thèse 11

## II Modèles de flots protostellaires 13

### 1 Observables clés à reproduire 14

1.1 Morphologie des flots 14

1.1.1 Angle d'ouverture à la base du flot 14

1.1.2 Rapport d'aspect $q$ 15

1.1.3 Largeur totale de flot à $|v| \simeq 1$ km/s 16

1.2 Cinématique 17

1.2.1 Distribution masse-vitesse 17

1.2.2 "Loi de Hubble" 17

1.2.3 Faible superposition "bleu-rouge" 17

### 2 Modèle prédominant de flot entrainé par un vent grand angle 19

2.1 Présentation et succès du modèle 19

2.1.1 Modèle analytique de couche poussée par un vent radial (Shu et al. 1991) 19

2.1.2 Simulations MHD: effet du degré d'aplatissement et de la magnétisation du jet et du coeur 22
<table><tbody><tr><td>2.1.3</td><td>Modèle analytique simplifié avec vent-X "modifié" (Lee et al. 2001)</td><td>22</td></tr><tr><td>2.2</td><td>Critiques et limites : difficultés récentes du scénario</td><td>25</td></tr><tr><td>2.2.1</td><td>Des chocs internes pas assez incurvés</td><td>25</td></tr><tr><td>2.2.2</td><td>Un rapport d'aspect sous-estimé</td><td>25</td></tr><tr><td>2.2.3</td><td>Des âges trop courts ou des largeurs excessives</td><td>26</td></tr><tr><td>3</td><td>Flots poussés par un jet : un modèle trop rapidement critiqué</td><td>27</td></tr><tr><td>3.1</td><td>Principe général des chocs d'étrave</td><td>27</td></tr><tr><td>3.2</td><td>Modèles de chocs d'étrave dans un milieu homogène</td><td>28</td></tr><tr><td>3.2.1</td><td>Choc d'étrave à la tête du jet: modèle analytique d'Ostriker et al. (2001)</td><td>28</td></tr><tr><td>3.2.2</td><td>Le long d'un jet variable : modèle analytique de Tabone et al. (2018)</td><td>30</td></tr><tr><td>3.3</td><td>Succès du modèle de flot entraîné par un jet</td><td>34</td></tr><tr><td>3.3.1</td><td>Distribution masse-vitesse m ∝ v<sup>γ</sup></td><td>34</td></tr><tr><td>3.3.2</td><td>"Loi de Hubble" dans le diagramme position-vitesse</td><td>35</td></tr><tr><td>3.3.3</td><td>Morphologie des chocs en H<sub>2</sub></td><td>35</td></tr><tr><td>3.3.4</td><td>Evite les largeurs excessives</td><td>35</td></tr><tr><td>3.4</td><td>Critiques avancées contre le modèle (dans un milieu uniforme)</td><td>36</td></tr><tr><td>3.4.1</td><td>Collimation excessive</td><td>36</td></tr><tr><td>3.4.2</td><td>Vitesse trop faible dans les flots larges (W ≃ 10 000 au)</td><td>37</td></tr><tr><td>3.4.3</td><td>Superposition "bleu-rouge" excessive</td><td>37</td></tr><tr><td>3.5</td><td>Face aux critiques : idée d'un jet pulsé multimillénaire dans un cœur stratifié</td><td>37</td></tr><tr><td>3.5.1</td><td>Cœur stratifié: des rapports d'aspect plus ressemblants</td><td>37</td></tr><tr><td>3.5.2</td><td>Jet variable : des vitesses élevées même dans les flots âgés</td><td>38</td></tr><tr><td>3.5.3</td><td>Jet multimillénaire : une superposition "rouge-bleu" diminuée ?</td><td>38</td></tr><tr><td>III</td><td>Le code MPI-AMRVAC</td><td>39</td></tr><tr><td>1</td><td>Le code MPI-AMRVAC et équations résolues</td><td>39</td></tr><tr><td>1.1</td><td>Résumé des caractéristiques</td><td>39</td></tr><tr><td>1.2</td><td>Equations d'Euler et simulations hydrodynamiques</td><td>41</td></tr><tr><td>1.2.1</td><td>Des traceurs pour suivre la matière</td><td>42</td></tr><tr><td>1.3</td><td>Modularité du code : un grand potentiel d'applications</td><td>42</td></tr><tr><td>2</td><td>Méthodes, schémas de résolution numériques et machine Tycho</td><td>43</td></tr><tr><td>2.1</td><td>Approche des volumes finis et solutions numériques</td><td>43</td></tr><tr><td>2.2</td><td>Schémas numériques et limiteurs de pentes</td><td>44</td></tr><tr><td>2.3</td><td>Calculateur Tycho (MESOPSL)</td><td>47</td></tr><tr><td>IV</td><td>Application du code MPI-AMRVAC pour la modélisation de jets protostellaires variables dans un cœur stratifié</td><td>49</td></tr><tr><td>1</td><td>Modules utilisateurs</td><td>50</td></tr><tr><td>1.1</td><td>Travail de modélisation durant la thèse</td><td>50</td></tr><tr><td>1.2</td><td>Liste des modules utilisateurs</td><td>51</td></tr></tbody></table>
<table>
  <tbody>
    <tr>
      <td>2</td>
      <td>Modélisation des équations d'état hydrodynamiques</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.1</td>
      <td>Équations d'état</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.2</td>
      <td>Composition du gaz</td>
      <td>52</td>
    </tr>
    <tr>
      <td>2.2.1</td>
      <td>Formules générales</td>
      <td>53</td>
    </tr>
    <tr>
      <td>2.2.2</td>
      <td>Gaz parfait ionisé</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.2.3</td>
      <td>Gaz parfait doublement ionisé</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.2.4</td>
      <td>Gaz parfait complètement moléculaire</td>
      <td>54</td>
    </tr>
    <tr>
      <td>2.2.5</td>
      <td>Gaz parfait complètement atomique et neutre</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3</td>
      <td>Enveloppe environnante modélisée</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.1</td>
      <td>Liste des paramètres libres</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.2</td>
      <td>Conditions initiales</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.2.1</td>
      <td>Champs de densité et de vitesse du cœur</td>
      <td>55</td>
    </tr>
    <tr>
      <td>3.3</td>
      <td>Conditions aux bords de l'enveloppe</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3.3.1</td>
      <td>Cellules fantômes</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3.3.2</td>
      <td>Conditions aux bords en R = 0, R = R<sub>max</sub>, z = z<sub>max</sub>, z = 0</td>
      <td>56</td>
    </tr>
    <tr>
      <td>3.3.3</td>
      <td>Conditions aux bords pour le cœur en z = 0</td>
      <td>58</td>
    </tr>
    <tr>
      <td>4</td>
      <td>Modélisation du jet pulsé</td>
      <td>61</td>
    </tr>
    <tr>
      <td>4.1</td>
      <td>Liste des paramètres libres</td>
      <td>61</td>
    </tr>
    <tr>
      <td>4.2</td>
      <td>Valve d'injection: domaine et ouverture angulaire</td>
      <td>61</td>
    </tr>
    <tr>
      <td>4.3</td>
      <td>Champs de densité et de vitesse initiaux dans la valve d'injection</td>
      <td>64</td>
    </tr>
    <tr>
      <td>4.3.1</td>
      <td>Pour une valve rectangulaire</td>
      <td>64</td>
    </tr>
    <tr>
      <td>4.3.2</td>
      <td>Pour une valve trapézoïdale</td>
      <td>65</td>
    </tr>
    <tr>
      <td>4.4</td>
      <td>Variabilité du jet: un paramètre physique encore plus ou moins libre</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.4.1</td>
      <td>Profil de variabilité sinusoïdal</td>
      <td>67</td>
    </tr>
    <tr>
      <td>4.4.2</td>
      <td>Profil de variabilité en dent de scie (sawtooth)</td>
      <td>68</td>
    </tr>
    <tr>
      <td>5</td>
      <td>Termes sources</td>
      <td>68</td>
    </tr>
    <tr>
      <td>5.1</td>
      <td>Force de contrebalancement F&#773;<sub>p</sub> et terme source E<sub>p</sub></td>
      <td>69</td>
    </tr>
    <tr>
      <td>5.1.1</td>
      <td>Evaluation de F&#773;<sub>p</sub> : méthode numérique</td>
      <td>69</td>
    </tr>
    <tr>
      <td>5.1.2</td>
      <td>Evaluation de F&#773;<sub>p</sub> : méthode analytique</td>
      <td>70</td>
    </tr>
    <tr>
      <td>5.1.3</td>
      <td>Sens physique de l'énergie de contrebalancement E<sub>p</sub>(t)</td>
      <td>70</td>
    </tr>
    <tr>
      <td>5.1.4</td>
      <td>Évaluation de E<sub>p</sub></td>
      <td>72</td>
    </tr>
    <tr>
      <td>5.2</td>
      <td>Refroidissement -n<sup>2</sup><sub>H</sub> Λ(T) par désexcitations collisionnelles</td>
      <td>72</td>
    </tr>
    <tr>
      <td>6</td>
      <td>Tests de configuration pour notre problème physique</td>
      <td>73</td>
    </tr>
    <tr>
      <td>6.1</td>
      <td>Choix des schémas numériques et limiteurs de pente</td>
      <td>73</td>
    </tr>
    <tr>
      <td>6.1.1</td>
      <td>Critères de sélection</td>
      <td>74</td>
    </tr>
    <tr>
      <td>6.1.2</td>
      <td>Résultats</td>
      <td>74</td>
    </tr>
    <tr>
      <td>6.1.3</td>
      <td>Schéma d'intégration en temps</td>
      <td>75</td>
    </tr>
    <tr>
      <td>6.2</td>
      <td>Résolution spatiale: influence du niveau de raffinement par AMR</td>
      <td>75</td>
    </tr>
    <tr>
      <td>6.2.1</td>
      <td>Critères de sélection</td>
      <td>76</td>
    </tr>
    <tr>
      <td>6.2.2</td>
      <td>Résultats</td>
      <td>76</td>
    </tr>
    <tr>
      <td>6.3</td>
      <td>Paramètres de parallélisation avec MPI</td>
      <td>77</td>
    </tr>
    <tr>
      <td>6.3.1</td>
      <td>Test d'extensibilité forte et loi d'Amdhal</td>
      <td>77</td>
    </tr>
  </tbody>
</table>
<table><tr><td>6.3.2 Méthode</td><td>78</td></tr><tr><td>6.3.3 Résultats</td><td>78</td></tr></table>

## V Flots poussés par un jet variable: Influence de l'environnement et de l'âge

81

1 D'un environnement uniforme à un cœur stratifié

82

1.1 Paramètres du jet et du domaine
82

1.2 Choix des modèles de stratification du cœur
84

1.3 Résultats
85

1.3.1 Coupes de densité à t = 300 – 600 ans
85

1.3.2 Champs de vitesse
87

1.3.3 Diagrammes position-vitesse
87

1.3.4 Distributions masse-vitesse
91

1.4 Comparaison au flot poussé par un vent grand angle
92

2 Influence des paramètres libres dans le cas d'un jet conique

94

2.1 Configuration numérique
94

2.1.1 Dispersion angulaire
94

2.1.2 Refroidissement radiatif et grille adaptative
95

2.2 Grilles de modèles
95

2.3 Résultats : effet des paramètres libres
98

2.3.1 Cartes de vitesses
103

2.3.2 Forme de la cavité principale – méthode
103

2.3.3 Formes de la cavité principale – analyse
103

3 Simulations sur 10⁴ ans du flot poussé par un jet variable dans un cœur singulier aplati

104

3.1 Choix des paramètres de la simulation
104

3.2 Validation du calcul de la force de rappel
105

3.3 Comparaison avec les résultats aux temps courts
105

3.4 Evolution de la largeur et de l'ouverture des cavités
107

3.4.1 Décélération à long terme
107

3.4.2 Comparaison avec le modèle de flot poussé par un vent grand angle
109

3.4.3 Comparaison aux largeurs de flots observées
110

3.5 Superposition rouge-bleu et vitesse maximale
111

## VI Comparaison aux observations ALMA de flots typiques: HH46/47 et CARMA-7

115

1 Aniriniil : un code Python pour les cartes synthétiques

116

1.1 Cube de données
118

2 Comparaison avec le flot de HH46/47

120

2.1 Morphologie
120
<table><tr><td>2.2 Diagrammes PV transverses</td><td>121</td></tr><tr><td>2.2.1 Résultats</td><td>123</td></tr><tr><td>2.2.2 Comparaison au modèle de vent grand angle</td><td>123</td></tr><tr><td>2.3 Diagramme PV longitudinal</td><td>123</td></tr><tr><td>2.3.1 Résultats</td><td>123</td></tr><tr><td>2.3.2 Effet de la phase du jet</td><td>125</td></tr><tr><td>2.3.3 Comparaison au modèle de vent grand angle</td><td>125</td></tr><tr><td>2.4 Distribution masse-vitesse</td><td>125</td></tr><tr><td><b>3 Comparaison avec le flot de CARMA-7</b></td><td><b>128</b></td></tr><tr><td>3.1 Observations et paramètres adoptés</td><td>128</td></tr><tr><td>3.2 Diagrammes position-vitesse</td><td>129</td></tr><tr><td>3.2.1 Comparaison pour t = 500 ans</td><td>129</td></tr><tr><td>3.2.2 Facteurs favorisant un flot jeune</td><td>129</td></tr><tr><td>3.3 Quelques réserves et alternatives</td><td>132</td></tr><tr><td><b>4 Synthèse</b></td><td><b>133</b></td></tr><tr><td><b>VII Conclusion générale, réserves et perspectives</b></td><td><b>135</b></td></tr><tr><td><b>APPENDICES</b></td><td><b>149</b></td></tr><tr><td><b>A Force de contrebalancement d'un cœur singulier aplati</b></td><td><b>149</b></td></tr><tr><td><b>B Article publié en premier auteur</b></td><td><b>149</b></td></tr></table>
<figure><img src="image_4.png" /></figure>
# CHAPITRE I

## INTRODUCTION

### Sommaire

1 La formation d'une étoile de type solaire 1

1.1 Phase pré-stellaire 1

1.2 Phase protostellaire 2

1.2.1 Classe 0 2

1.2.2 Classe 1 4

1.2.3 Classe 2 4

1.2.4 Classe 3 5

2 Rôle des éjections protostellaires et questions ouvertes 5

2.1 Impact sur la formation des étoiles et sur le nuage moléculaire parent 5

2.2 Origine de l'éjection et moment cinétique 7

2.3 Impact sur la formation planétaire 7

2.4 Deux prescriptions de vent pour la rétroaction du flot protostellaire 7

2.5 Retour d'un scénario et objectifs de la thèse 9

2.6 Plan de la thèse 11

La signature la plus spectaculaire de la naissance d'une étoile est aussi souvent la première qui est observée, à savoir la formation d'un flot bipolaire de gaz moléculaire. Ce phénomène commence au début de la phase protostellaire, lors de l'assemblage de la masse stellaire (phase de Classe 0), persiste pendant la phase de dispersion de l'enveloppe (Classe 1) et est omniprésent au cours de la formation stellaire à toutes les masses (pour une revue, voir par exemple Frank et al. 2014).

## 1. La formation d'une étoile de type solaire

Dans cette section, je présente l'état de nos connaissances sur les étapes de la formation d'une étoile de type solaire telle que la nôtre. La Fig. 1 schématise les étapes de la formation d'une étoile de faible masse, telles qu'elles sont aujourd'hui consensuellement établies, avec les propriétés physiques caractéristiques et la durée de vie typique de chaque étape.

### 1.1. Phase pré-stellaire

Les étoiles se forment au sein des nuages moléculaires et plus précisément des "cœurs denses préstellaires". Avec des tailles de $r \sim 0.05$ pc et des densités typiques $n \sim 10^5$ cm$^{-3}$, ces structures de poussière et de gaz demeurent très froides, avec des températures avoisinant $T \sim 10 - 20$ K (Andre et al. 2000). Les forces thermiques, turbulentes et magnétiques en leur sein permettent initialement
de les maintenir en équilibre gravitationnel.

Une perturbation de cet équilibre entraîne l'effondrement d'une partie du nuage et une densification du cœur préstellaire. À un moment donné, le centre du cœur préstellaire devient optiquement épais et empêche au rayonnement thermique de s'en échapper, chauffant adiabatiquement le gaz comprimé. À ce stade, ce sont les gradients de pression thermiques qui compensent majoritairement la force gravitationnelle. Un premier cœur opaque hydrostatique de $r \sim 5$ au est formé.

Ensuite, la dissociation de l'hydrogène moléculaire brise l'équilibre du premier cœur et poursuit l'effondrement. Un second cœur hydrostatique de $r \sim 1-3$ au de rayon se forme alors, constituant ce que nous appelons une protoétoile et marquant par convention le début de la phase dite protostellaire (Larson 1969). L'âge $t = 0$ d'un objet prototellaire est compté à partir de ce moment.

## 1.2. Phase protostellaire

Cette phase, durant laquelle la protoétoile est enfouie à l'intérieur d'une enveloppe alimentant l'objet central, comporte 4 étapes (Classes 0 à 3). Elle se termine par le démarrage des réactions nucléaires dans le noyau de la protoétoile, qui marque le début de la séquence principale de l'étoile. Comme le cœur protostellaire initial possède un moment cinétique non-nul, le gaz de l'enveloppe finit par former un disque en rotation autour de la protoétoile (voir Ulrich 1976). Durant ma thèse, je ne me suis intéressé à l'enveloppe que jusqu'à $\sim 0.1$ pc, ainsi qu'au vent protostellaire "primaire" collimaté et qui permet de pousser le flot lent. Cependant, il ne faut pas oublier qu'enveloppe, protoétoile et disque protostellaire continuent d'interagir entre eux durant la phase protostellaire, via divers processus physiques (hydrodynamiques, gravitationnels, thermiques, radiatifs, magnétiques, chimiques).

### 1.2.1. Classe 0

La phase de Classe 0 est la phase principale d'accrétion durant laquelle l'enveloppe reste encore plus massive que le cœur hydrostatique enfoui, le rendant invisible à l'infrarouge proche et moyen ($\lambda < 30$ µm) en raison de l'extinction. La distribution spectrale d'énergie (SED) $\alpha = \frac{d \log(\lambda F_{\lambda})}{d \log \lambda}$ des Classes 0 est celle d'un corps gris (un corps noir approché) en raison de l'émission de la poussière chauffée de l'enveloppe, avec une température bolométrique<fnref n="1" /> valant $T_{bol} < 70$ K (Myers & Ladd 1993). Un autre paramètre intéressant est le rapport entre la luminosité sub-millimétrique $L_{subm}$ et la luminosité bolométrique $L_{bol}$, qui donne une idée de la différence de masse entre l'enveloppe et l'objet central accrété (Andre et al. 2000). Les Classes 0 sont caractérisées par un rapport supérieur à 0.5%. Selon des études statistiques sur ces objets et des estimations de taux d'accrétions, leur durée de vie typique est de $t \sim 10^4$ ans pour des taux d'accrétion $M_{acc} \sim 10^{-6} - 10^{-5} M_{\odot} \text{ yr}^{-1}$ (Bontemps et al. 1996; Andre et al. 2000).

Comme pour l'exemple de HH212 (Fig. 2), les observations de Classes 0 montrent la coexistence de deux composantes distinctes associées à un phénomène de perte de masse du système central. D'une part, nous voyons des flots (cavités) de gaz moléculaires lents ($\lesssim 30$ km/s projetées), dont la largeur augmente avec l'âge de la source. Ces cavités se resserrent au niveau de l'objet central avec une faible collimation, c'est-à-dire un angle d'ouverture de $\alpha_{proj} \sim 25^\circ - 100^\circ$ à la base du flot.

D'autre part, les interféromètres millimétriques et submillimétriques ont aussi montré l'existence de jets moléculaires plus rapides ($\ge 100$ km/s déprojetées) en CO, H₂, et SiO (Bachiller et al. 1991; Guilloteau et al. 1992). Ils sont aussi plus collimatés que les flots moléculaires ($\alpha_{proj} \sim$ quelques

<footnote n="1">1. Température de corps noir avec la même fréquence moyenne que la SED</footnote>
<figure><img src="image_5.png" /></figure>

FIGURE 1 – Etapes de la formation d'une étoile de faible masse, depuis l'effondrement du cœur préstellaire jusqu'à la phase de Classe 3. Les échelles de temps et de dimensions spatiales des systèmes, et leurs propriétés physiques typiques (température, densité, masses, distributions spectrales d'énergie ou SED) sont représentées. Adaptée de la thèse de Marc Joos (https://tel.archives-ouvertes.fr/tel-00766166/).

degrés). Cette émission présente souvent des formes de nœuds (knots), interprétées comme des chocs générés par la variabilité temporelle en vitesse de l'éjection (Raga et al. 1990).
<figure><img src="image_6.png" /><figcaption>FIGURE 2 – Émission $^{12}$CO(2–1) (rouge) et H$_2$ (vert) de la Classe 0 HH212 (Lee et al. 2015). Chaque canal (en bas à gauche de chaque panneau) montre une vitesse observée projetée différente.</figcaption></figure>

### 1.2.2. Classe 1

Durant la phase de Classe 1, la masse de l'enveloppe descend en dessous de celle de la protoétoile. Le système central devient alors observable en infrarouge dans la SED (Fig. 1), avec une température bolométrique comprise entre 70 K et 700 K. Pour cette étape, dont la durée typique est de $10^5$ans, le taux d'accrétion de l'enveloppe sur le disque est estimé à $M_{acc} \sim 10^{-7} - 10^{-6} M_{\odot} \text{yr}^{-1}$ (Bontemps et al. 1996; Antoniucci et al. 2008).

En prenant le cas de HH46/47 comme illustration (Fig. 3), jets et flots coexistent toujours durant cette étape. Plus précisément, à part en H$_2$, les jets de Classe 1 sont peu abondants en molécules. Ils sont observés en infrarouge et en H$\alpha$, [OI], [SII], et [NII]. Ces composantes atomiques forment une structure à haute vitesse (jusqu'à $\sim 150 \text{ km/s}$) qui est plus collimatée qu'en H$_2$ (Agra-Amboage et al. 2014, dans le cas de DG Tau B). Au-delà de $\sim 0.1 \text{ pc}$, le jet collimaté est moins brillant, mais des chocs d'étraves sont toujours visibles en H$_2$ et en optique jusqu'à plusieurs parsecs.

En revanche, les flots lents observés restent moléculaires, même si les cavités larges toujours détectables sont moins brillantes que durant la Classe 0 (comme pour HH46/47, Fig. 3).

### 1.2.3. Classe 2

Au moment où toute l'enveloppe est dispersée, la protoétoile entre en phase de Classe 2, devenant une TTauri classique (cTTs). Des jets atomiques très rapides (jusqu'à $\sim 400 \text{ km/s}$) sont toujours
HH46/47

<figure><img src="image_7.png" /><figcaption>FIGURE 3 – Image composite du flot et jet de Classe 1 HH46/47, enfoui dans son globule de Bok natal.</figcaption></figure>

observés, tandis que seules des cavités CO et H₂ évasées et à basse vitesse (~10 km/s) peuvent encore être observées autour de la base de certains jets comme DG Tau ou HH30 (Louvet et al. 2018). L'émission en infrarouge dans la SED est dominée par le disque d'accrétion, avec des taux d'accrétion estimés à $\dot{M}_{acc} \sim 10^{-8} M_{\odot} yr^{-1}$. Le temps de dissipation du disque nous donne la durée de vie de cette phase, estimée à quelques $10^6$ans.

### 1.2.4. Classe 3

Dernière étape de la formation stellaire, les objets de Classes 3, aussi appelés TTauris faibles (wTTs), voient la fin de l'accrétion sur ces objets et un disque protostellaire optiquement mince, voire absent. À ce stade, l'activité d'éjection s'est arrêtée, avec la disparition des cavités moléculaires et des jets observés dans les phases précédentes. D'ailleurs, l'absence de jet va dans le sens d'une connexion directe entre phénomène d'éjection et d'accrétion.

## 2. Rôle des éjections protostellaires et questions ouvertes

### 2.1. Impact sur la formation des étoiles et sur le nuage moléculaire par-ent

Compte tenu de leur grande taille, de leur flux de masse et d'impulsion typiques élevés, de leur ubiquité, et de leur durée, les flots moléculaires joueraient un rôle clé dans la formation d'étoiles à petite et à grande échelle : des simulations numériques récentes (voir par exemple, Krumholz & Federrath 2019, pour une revue) suggèrent qu'ils pourraient être le principal agent de rétroaction
qui fixe la masse stellaire finale et l'efficacité dans la conversion masse enveloppe-étoile (core-to-star efficiency, via l'élimination de la matière du nuage parent).

Ils réguleraient le pic dans la fonction initiale de masse (Initial Mass Function ou IMF) ou la fraction de multiplicité des étoiles (fraction de systèmes stellaires liés formés de plusieurs étoiles). Plus précisément, la distribution de masse des cœurs préstellaires (CMF) est décalée par rapport à celle de la masse des étoiles formées (IMF) (Könyves et al. 2015). Ce décalage peut se traduire par une efficacité de formation des étoiles, avec ~30% de la masse de l'enveloppe en effondrement qui arrive au final sur la protoétoile. Comprendre d'où vient une faible efficacité est une question qui reste ouverte, mais qui serait surtout connecté aux phases initiales de Classe 0 durant laquelle la majorité de la masse de l'enveloppe est accrétée.

Ils réguleraient aussi l'efficacité de formation des étoiles (star formation efficiency) à l'échelle des amas, via la perturbation des écoulements d'effondrement, ainsi que le réapprovisionnement de la turbulence. De plus, en emportant de la matière du système central jusqu'à plusieurs parsecs dans l'enveloppe environnante, les conditions et compositions physico-chimiques de cette dernière peuvent être modifiées dans le temps.

<figure><img src="image_8.png" /></figure>

FIGURE 4 – Cartes de densité de simulations d'évolution d'un nuage moléculaire géant, avec (à droite) et sans (à gauche) flots moléculaires. Les symboles représentent les protoétoiles en formation. Ici, les flots sont simulés suivant le modèle de vent X grand angle de Shu et al. (1991). Images tirées de Guszejnov et al. (2021).

La Fig. 4 montre un exemple de simulations de rétroaction des flots par Guszejnov et al. (2021), à l'échelle du nuage moléculaire. Selon leurs résultats, inclure les flots protostellaires diminue le taux de formation de protoétoiles dans le nuage d'un facteur 2 – 3, comparé au cas sans flots. En fonction du modèle de vent "primaire" utilisé pour pousser les flots (Sect. I.2.4 et Chapitre II), les cavités creusées et l'impact sur l'effondrement du nuage sont différents, et peuvent ainsi impacter différemment la formation stellaire.

Il est donc pertinent et important de comprendre en détail comment les jets et flots moléculaires sont produits, ainsi que leur rétroaction jusqu'à grande échelle. En effet, nous pouvons ainsi mieux saisir leur impact sur l'efficacité et le taux de formation stellaire, sur la masse finale de l'étoile formée, ainsi que sur l'évolution dynamique et chimique de chaque composante, à savoir : le disque
et la protoétoile en formation à petite échelle, et le nuage parent jusqu'à grande échelle.

## 2.2. Origine de l'éjection et moment cinétique

Les phénomènes d'éjection contribuent à emporter une partie du moment cinétique hérité de l'enveloppe parent en effondrement. Cependant, la compréhension exacte du mécanisme à l'origine de l'éjection et de l'extraction du moment cinétique reste sujet à débats.

Parmi les modèles avancés, nous avons celui du vent de disque MHD magnéto-centrifuge (Blandford & Payne 1982). L'ancrage du champ magnétique dans la matière le long du disque en rotation entraîne la torsion du champ et une éjection verticale du moment cinétique par accélération magnétocentrifuge. Cette éjection prend la forme d'un jet auto-collimaté par le champ magnétique toroïdal externe, ce qui permet de reproduire très efficacement la forme quasi-cylindrique des jets. Un autre modèle proposé par Shu et al. (1991) est celui du "vent X", qui est aussi magnéto-centrifuge. Mais contrairement aux vents de disque MHD, le vent X est lancé depuis le bord interne du disque. Les lignes de champ magnétiques ancrées dans la protoétoile peuvent aussi générer un vent stellaire (Decampli 1981). Enfin, l'interaction non-stationnaire entre magnétosphère protostellaire et champ magnétique du disque peut induire une variabilité dans le jet, tout en freinant la rotation de la protoétoile. Dans la réalité, il est possible que ces modèles contribuent simultanément, avec plus ou moins d'importance, à différentes échelles et étapes de la formation stellaire.

Chaque modèle d'éjection a un impact différent sur l'enveloppe à plus grande échelle (voir Chapitre II). En outre, malgré la résolution élevée des dernières observations (quelques 0.1"), la région d'éjection nous échappe encore. Ainsi, en comparant les simulations numériques et les observations récentes à l'échelle de l'enveloppe en interaction avec le vent "primaire", il est possible de mieux comprendre rétroactivement et en détails les mécanismes d'éjection (et les contributions de chaque modèle).

## 2.3. Impact sur la formation planétaire

Dans l'hypothèse où les jets et flots tirent leur réservoir de matière directement du disque en rotation, il est clair que la formation des planètes autour de la future étoile peut aussi être impactée par les phénomènes d'éjection. En effet, le réservoir de gaz des futures planètes est éjecté, tandis que les conditions en densité des disques sont affectées (Suzuki et al. 2016), impactant la migration des planètes (Ogihara et al. 2018). Pour mieux saisir en détail le scénario de formation des planètes, il est ainsi nécessaire de bien comprendre l'évacuation du gaz du disque entourant la protoétoile, et donc de mieux comprendre les phénomènes d'éjection.

## 2.4. Deux prescriptions de vent pour la rétroaction du flot protostellaire

L'effet exact de la rétroaction du flot dépend de la structure supposée du vent protostellaire "primaire" qui permet de pousser le flot lent. Dans la communauté scientifique, deux principaux modèles de vent sont actuellement utilisés. Tous deux impliquent un jet rapide et dense le long de l'axe du flot (comme couramment observé dans les traceurs atomiques ou moléculaires, voir Frank et al. 2014), mais diffèrent fortement par la quantité de mouvement injectée à des angles d'injection plus larges.

La première prescription de vent est celle de Shu et al. (1991, 1995) et Matzner & McKee (1999, ci-après MM99), qui est aussi la plus fréquemment utilisée dans les simulations incluant
la rétroaction sur l'enveloppe. Elle suppose un vent grand angle en expansion radiale à vitesse constante (≈ 100 km/s) à tous les angles, avec une forte augmentation de densité le long et près de l'axe. Celle-ci serait responsable de l'apparition du "jet" axial qui est observé. Cette structure asymptotique a d'abord été obtenue pour un "vent X" magnétiquement accéléré à partir du bord interne du disque (Shu et al. 1995) et ne s'applique qu'aux vents magnétohydrodynamiques (MHD) radialement lancés à partir d'une région étroite (MM99). En se rapprochant de l'équateur, le flux d'impulsion du vent constitue toujours une fraction importante de celui obtenu dans le cas d'un vent isotrope<fnref n="1" /> et peut avoir un impact direct sur l'effondrement équatorial de l'enveloppe.

Le deuxième modèle de vent, motivé par des simulations MHD et observations plus récentes, suppose que le jet axial rapide est entouré d'un vent de disque plus lent, éjecté dans un angle solide limité (Federrath et al. 2014; Rohde et al. 2019). La rétroaction est alors dominée par le jet, avec un impact bien moindre sur les régions équatoriales que dans la prescription de MM99. Il a également été démontré que la variabilité temporelle, sous la forme de sursauts épisodiques, affecte la rétroaction des flots (Rohde et al. 2019).

En principe, des simulations MHD réalistes de lancement de vent protostellaire devraient fournir la meilleure prescription de vent à adopter. Cependant, la structure du vent simulé dépend d'effets complexes qui font encore l'objet d'intenses recherches et débats, tels que l'interaction magnétosphère-disque, la distribution du flux magnétique retenu dans le disque longtemps après sa formation, la viscosité et la résistivité turbulente, ainsi que des effets de MHD non-idéale (voir par exemple, Ireland et al. 2021; Ferreira & Deguiran 2013; Béthune et al. 2017).

Une approche indépendante pour déterminer le modèle de vent le plus réaliste pour les études de rétroaction des flots consiste à simuler les propriétés de flots poussés jusqu'à des échelles atteignant celles du cœur protostellaire $≤ 0,1$ pc (où la structure de densité ambiante est dominée par l'auto-gravité et n'est pas encore perturbée par les inhomogénéités des nuages), pour ensuite voir quel modèle de vent reproduit le mieux la morphologie et la cinématique des flots observés. Durant cette thèse, nous avons aussi adopté cette approche au moment de confronter nos propres simulations sur l'un des modèles avec celles du modèle "concurrent". Une telle comparaison a été effectuée très tôt pour deux extrêmes dans la modélisation de la collimation du vent : le modèle de vent X grand angle (Shu et al. 1995; Matzner & McKee 1999) et celui d'un jet collimaté "pur" qui entraîne le flot le long de chocs d'étraves en expansion (Masson & Chernin 1993; Raga & Cabrit 1993). Des succès et des critiques, ainsi que des mises en garde, ont été identifiées dans chaque cas, sur la base des premières observations (Cabrit et al. 1997; Lee et al. 2001; Arce et al. 2007). Celles-ci sont brièvement résumées et remises à jour ci-dessous.

Premièrement, dans les modèles de flots entraînés par un vent X à grand angle, le milieu ambiant est supposé avoir une forte diminution de densité (en $1/r^2$), ainsi qu'un degré modéré d'aplatissement magnétique pour l'enveloppe de gaz environnante. En outre, le gaz de l'enveloppe est supposé se mélanger instantanément avec la matière du vent choqué (voir Chapitre II). La couche balayée se dilate ensuite radialement de façon auto-similaire, de façon à reproduire plusieurs caractéristiques observées des flots moléculaires : à savoir la cinématique en "loi de Hubble" $V ∝ z$ des flots, la distribution masse-vitesse avec une pente de loi de puissance $\gamma \simeq -2$ (avant correction d'opacité), ainsi que des formes paraboliques avec de larges angles d'ouverture à la base du flot (Shu et al. 1991; Li & Shu 1996; Lee et al. 2001; Shang et al. 2006, 2020). Cependant, une mise en garde intrinsèque vis-à-vis de ce modèle est qu'une vitesse de vent uniforme injectée à tous les angles prédit des chocs internes beaucoup plus plats que ceux observés dans le gaz en $H_2$ excité par des chocs le long de l'axe des flots. La forme incurvée observée en $H_2$ nécessite une chute brutale

<footnote n="1">1. la fraction est de $1/\ln(2/\theta_0) \simeq 1/5$ pour un angle de collimation du "jet" $\theta_0 \simeq 0.01$ rad, cf. Eq. (2) dans MM99</footnote>
de la pression du vent loin de l'axe (voir les discussions dans Lee et al. 2001; Arce et al. 2007).

En revanche, les simulations de chocs d'étraves poussés par des jets reproduisent la morphologie incurvée des fronts de chocs internes nettement observés en $H_2$ (e.g., Suttner et al. 1997; Völker et al. 1999). Elles reproduisent aussi les structures caractéristiques en forme "d'éperon" (spur-like) associées à ces fronts de chocs et observées dans les diagrammes position-vitesse en CO (Lee et al. 2000) (diagnostics décrits dans le Chapitre II). Ils peuvent également reproduire les relations masse-vitesse observées (Downes & Cabrit 2003; Moraghan et al. 2008) ainsi que les "lois de Hubble" apparentes observées lorsque le jet est variable et en précession (Völker et al. 1999; Rohde et al. 2019). Cependant, le modèle de flots poussés par un jet a été fortement critiqué pour avoir prédit des cavités trop allongées sur de longues échelles de temps (Ostriker et al. 2001), une superposition excessive des émissions décalées vers le bleu et vers le rouge sur une large gamme d'inclinaisons du flot (Lee et al. 2001), ainsi que des vitesses ~ 0.03 km/s trop faibles lorsque le choc d'étrave s'est étendu à des largeurs de flots typiques de 10 000 au (Arce et al. 2007). En conséquence, il est admis que les jets ne peuvent à eux seuls expliquer les flots avec un grand angle d'ouverture, comme cela a été avancé pour les flots de plusieurs Classes 1 évoluées.

## 2.5. Retour d'un scénario et objectifs de la thèse

Ainsi, une question générale qui motive le sujet de ma thèse est de savoir quel vent protostellaire "primaire" mettre pour décrire les flots observés à grande et à petite échelle : vent grand angle ou jet collimaté ? Et quels impacts sur l'enveloppe environnante nos résultats prédisent-t-ils ?

Afin de combiner les points forts de chaque modèle, une structure "à double vent" a été invoquée avec à la fois un jet rapide et un vent grand angle plus lent, où ce dernier dominerait progressivement à temps longs (Yu et al. 1999; Velusamy & Langer 1998; Arce et al. 2007; Zapata et al. 2014). Cependant, il existe plusieurs bonnes raisons de reconsidérer les chocs d'étraves poussés par un jet "pur" comme étant potentiellement à l'origine des flots moléculaires.

Tout d'abord, les critiques de Ostriker et al. (2001) et Lee et al. (2001) étaient basées sur des modèles de chocs d'étrave poussés par un jet dans un milieu ambiant uniforme ou quasi uniforme (avec au plus un facteur 2 de variation de la densité sur le domaine de calcul). En contraste, une forte diminution radiale de la densité (en $1/r^2$) est supposée dans les modèles de flots poussés par un vent grand angle, afin de produire l'accélération apparente en "loi de Hubble" (Shu et al. 1991; Lee et al. 2001). Une telle diminution est attendue aux échelles du cœur protostellaire $\le 0,1$ pc, en raison de l'autogravité. Raga & Cabrit (1993) et Cabrit et al. (1997) ont montré qu'il pouvait produire un angle d'ouverture plus large pour les chocs d'étraves poussés par les jets, et qui est davantage similaire aux observations. Ce résultat est confirmé par Moraghan et al. (2008) pour des jets stationnaires sur des échelles de temps < 1000 ans. Maintenant qu'une stratification à forte densité est largement confirmée par les observations de cœurs protostellaires (e.g., van Kempen et al. 2009, pour HH46/47) et que les capacités des simulations numériques se sont grandement améliorées, il est important d'explorer l'effet que peut avoir un espace de paramètres plus large et des échelles de temps plus longues sur les chocs d'étraves poussés par un jet, comparé aux possibilités au moment de la première étude de Moraghan et al. (2008).

Deuxièmement, le problème soulevé par Arce et al. (2007) lié à une vitesse de choc insuffisante à de grandes largeurs de chocs d'étrave ne s'applique plus avec un jet pulsé. De nouveaux chocs internes générés par la variabilité du jet réapprovisionnent le choc d'étrave externe lent (le premier formé) avec du gaz plus rapide (à des vitesses observables, voir Raga & Cabrit 1993; Völker et al.).
1999). De plus, l'interaction entre les chocs d'étrave successifs diminuera la vitesse transverse de ces derniers, pouvant atténuer la superposition bleu et rouge excessive prédit à des âges courts par Lee et al. (2001). Des observations récentes de l'interféromètre Atacama Large Millimeter Array (ALMA) ont révélé de multiples chocs d'étrave en H₂ le long de l'axe des flots, et qui sont connectés à des cavités en CO imbriquées le long des flancs des flots (par exemple, HH212 et HH46/47 dans Lee et al. 2015; Zhang et al. 2016). Bien que certains flots aient été modélisés comme des couches poussées par un vent grand angle (Zhang et al. 2019), il est tout aussi important que des prédictions similaires soient explorées pour des jets imbriqués sur des âges réalistes de $\simeq 10^4$ ans.

Troisièmement, des études observationnelles récentes montrent que les flots CO sont plus allongés et collimatés qu'on ne le pensait initialement : les cartes CO à l'échelle du nuage parent, ainsi que les relevés en imagerie optique et infrarouge, montrent par exemple que 40% des 20 flots observés par Narayanan et al. (2012), dans un échantillon non biaisé de sources, font plus d'un parsec de long (Frank et al. 2014, et les références qui s'y trouvent). Un exemple remarquable est le flot de Classe 1 B5-IRS1. Il a un angle d'ouverture total (et projeté) $\alpha \simeq 100°$ à la base, jusqu'alors avancé comme étant la preuve d'un vent grand angle. Cependant, des cartes CO récentes non biaisées du nuage parent révèlent que chaque lobe s'étend en réalité (au moins) jusqu'à 2.2 pc de la source centrale (Arce et al. 2010). Le rapport longueur-sur-largeur dépasse alors $q > 11$ (voir les cartes dans Frank et al. 2014). Cette nouvelle observation est incompatible avec les modèles actuels de flots poussés par un vent X grand angle (qui prédissent un rapport d'aspect $q < 3.5$ pour un angle d'ouverture de base $\alpha \ge 70°$, voir Table 2 de Shang et al. 2020). Cette taille de flots atteignant le parsec apparaît plutôt comme étant évocatrice de flots poussés par des jets.

Les estimations d'angles d'ouverture de flots sont également affectés par de nombreux biais : la résolution angulaire, la hauteur à laquelle ils sont mesurés (Velusamy & Langer 1998; Velusamy et al. 2014) et l'inclinaison du flot (les angles apparaissent plus larges dans les flots vus plus près par le pôle ou pole-on). Deux études récentes du nuage moléculaire d'Orion minimisent ces biais. En effet, ils fournissent des mesures uniformes à la même résolution linéaire (élevée) et à la même hauteur projetée, sur des échantillons de sources sélectionnés au hasard. La première est une étude ALMA sur les largeurs de 22 flots CO (principalement de Classe 0, Dutta et al. 2020). Elle montre 50% de flots décrivant des angles d'ouverture totaux et projetés contenus dans l'intervalle $\alpha = 25° - 65°$ à une altitude projetée de $z_{proj} = 800$ au. Dans un autre échantillon de 29 cavités de flots plus anciens (principalement de Classe 1) obtenues en lumière diffuse avec le Hubble Space Telescope (Habel et al. 2021), on voit que 50% sont dans l'intervalle $\alpha = 8° - 46°$ à $z_{proj} = 8000$ au, tandis que la fraction des sources ponctuelles (à savoir vues depuis la partie inférieure de l'intérieur de la cavité) suggèrent un angle d'ouverture maximal déprojeté $\alpha_{deproj} \le 70°$. Par conséquent, les flots CO semblent en moyenne être plus collimatés qu'on ne le croyait auparavant. Il devient nécessaire d'étudier si le scénario d'un jet pur (lancé dans une enveloppe stratifiée) pourrait reproduire les largeurs observées typiques, avant de conclure sur une possible contribution dominante d'un vent grand angle.

Au cours de ma thèse, nous avons examiné cette question en présentant les premières simulations à haute résolution de chocs d'étrave poussés par un jet "pur" à travers des enveloppes fortement stratifiées, jusqu'à des âges de simulation atteignant 10 000 ans, ainsi que sur une échelle de 0.1 pc. Pour la première fois dans les simulations de jets, nous considérons la même structure de cœur singulier aplati que dans le modèle concurrent de flot poussé par un vent X grand angle (Li & Shu 1996; Lee et al. 2001). Ainsi, le travail présenté dans ce mémoire se concentre sur une première exploration d'un modèle de flot poussé par un jet, tout en démontrant son potentiel et répondre aux critiques jusqu'ici avancées.
## 2.6. Plan de la thèse

Dans le chapitre II, nous décrivons plus en détails l'état de l'art sur les deux paradigmes (flots poussés par un vent grand angle ou par un jet collimaté), les modèles analytiques et les simulations déjà effectuées, ainsi que les observables clés permettant de caractériser chaque modèle. Ce sera aussi l'occasion de présenter en détail les succès et critiques avancés à ce jour vis-à-vis de chaque scénario et au regard des observations récentes, et qui ont motivé l'étude réalisée durant ma thèse.

Dans le chapitre III, nous présentons le code de simulation MPI-AMRVAC employé pour nos simulations numériques hydrodynamiques, ainsi qu'un résumé des méthodes et solveurs numériques qui sont disponibles et pertinents vis-à-vis de notre problème physique. J'y présente aussi un outil Python que j'ai programmé et avec lequel je réalise les diagrammes position-vitesse des flots simulés.

Le chapitre IV développe comment nous avons établi le modèle physique à simuler, ainsi que le choix des méthodes, solveurs et paramètres numériques qui forment la configuration générique commune à nos simulations. Par exemple, nous présentons ici comment nous simulons le lancement du jet, sa variabilité temporelle, les conditions aux bords adoptées, ainsi que la stratification de l'enveloppe environnante.

Le chapitre V présente les résultats de nos simulations en se concentrant sur l'effet de l'âge et d'une stratification abrupte en densité 1D et 2D de l'enveloppe sur la forme et la cinématique des chocs d'étraves poussés par le jet. Nous utilisons pour cela les largeurs de flots prédites, les diagrammes position-vitesse et la relation masse-vitesse. Nous comparons ensuite nos simulations avec celles d'un modèle de vent grand angle lancé dans la même enveloppe avec le même flux de masse. Nous introduisons un petit angle d'ouverture du jet (alors dit conique) pour ensuite explorer l'effet de divers paramètres libres sur la forme et la cinématique de la cavité. Si les simulations de la première partie de ce chapitre se concentrent sur des âges courts < 1000 ans, la seconde partie du chapitre explore l'influence de l'âge sur les mêmes diagnostics morphologiques et cinématiques, sur des âges atteignant 10 000 ans. C'est aussi l'occasion de discuter comment nos simulations montrent que les critiques précédemment avancées vis-à-vis du modèle de flots poussés par un jet ne sont plus valables au regard de nos simulations. Nous confrontons ces critiques à nos résultats, tels que le maintien de vitesses suffisamment importantes le long des flancs des chocs d'étrave, la décélération à temps long de ces derniers, ou encore la diminution avec l'âge de la superposition bleu et rouge.

Dans le chapitre VI, nous confrontons qualitativement les résultats de nos simulations à temps longs, et dans un cœur singulier aplati, avec des observations ALMA en CO de flots bipolaires : HH46/47 et CARMA-7. À l'issue de nos simulations, les diagnostics morphologiques et cinématiques montrent un excellent accord qualitatif avec ces récentes observations. Les chapitres 5 et 6 montrent ainsi l'intérêt de reconsidérer le modèle de flots poussés par un jet "pur" comme scénario possible pour rendre compte des observations de flots protostellaire.

Enfin, le Chapitre VII présente les conclusions de ma thèse, une présentation des limites et mises en garde par rapport à notre modèle et nos résultats, ainsi qu'une mise en perspective des possibilités d'amélioration et d'études ultérieures.
<figure><img src="image_9.png" /></figure>
# CHAPITRE

## II

## MODÈLES DE FLOTS
PROTOSTELLAIRES

### Sommaire

1 Observables clés à reproduire 14

1.1 Morphologie des flots 14

1.1.1 Angle d'ouverture à la base du flot 14

1.1.2 Rapport d'aspect $q$ 15

1.1.3 Largeur totale de flot à $|v| \simeq 1$ km/s 16

1.2 Cinématique 17

1.2.1 Distribution masse-vitesse 17

1.2.2 "Loi de Hubble" 17

1.2.3 Faible superposition "bleu-rouge" 17

2 Modèle prédominant de flot entraîné par un vent grand angle 19

2.1 Présentation et succès du modèle 19

2.1.1 Modèle analytique de couche poussée par un vent radial (Shu et al. 1991) 19

2.1.2 Simulations MHD: effet du degré d'aplatissement et de la magnétisation du jet et du coeur 22

2.1.3 Modèle analytique simplifié avec vent-X "modifié" (Lee et al. 2001) 22

2.2 Critiques et limites : difficultés récentes du scénario 25

2.2.1 Des chocs internes pas assez incurvés 25

2.2.2 Un rapport d'aspect sous-estimé 25

2.2.3 Des âges trop courts ou des largeurs excessives 26

3 Flots poussés par un jet : un modèle trop rapidement critiqué 27

3.1 Principe général des chocs d'étrave 27

3.2 Modèles de chocs d'étrave dans un milieu homogène 28

3.2.1 Choc d'étrave à la tête du jet: modèle analytique d'Ostriker et al. (2001) 28

3.2.2 Le long d'un jet variable : modèle analytique de Tabone et al. (2018) 30

3.3 Succès du modèle de flot entraîné par un jet 34

3.3.1 Distribution masse-vitesse $m \propto v^\gamma$ 34

3.3.2 "Loi de Hubble" dans le diagramme position-vitesse 35

3.3.3 Morphologie des chocs en $H_2$ 35

3.3.4 Evite les largeurs excessives 35

3.4 Critiques avancées contre le modèle (dans un milieu uniforme) 36

3.4.1 Collimation excessive 36
<table>
  <tbody>
    <tr>
      <td>3.4.2</td>
      <td>Vitesse trop faible dans les flots larges (W ≃ 10 000 au)</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.4.3</td>
      <td>Superposition "bleu-rouge" excessive</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.5</td>
      <td>Face aux critiques : idée d'un jet pulsé multimillénaire dans un coeur stratifié</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.5.1</td>
      <td>Coeur stratifié: des rapports d'aspect plus ressemblants</td>
      <td>37</td>
    </tr>
    <tr>
      <td>3.5.2</td>
      <td>Jet variable : des vitesses élevées même dans les flots âgés</td>
      <td>38</td>
    </tr>
    <tr>
      <td>3.5.3</td>
      <td>Jet multimillénaire : une superposition "rouge-bleu" diminuée ?</td>
      <td>38</td>
    </tr>
  </tbody>
</table>

# 1. Observables clés à reproduire

Pour savoir quel modèle est le plus apte à reproduire les flots protostellaires observés, nous disposons d'un certain nombre d'observables clés. Celles-ci sont essentiellement observées dans les raies de $^{12}$CO ($J=1-0$ ou $J=2-1$). Mais nous disposons à ce jour de très peu de cartes de flots disponibles dans les isotopes $^{13}$CO et $^{18}$CO, car ces espèces sont bien moins abondantes et donc émettent ainsi bien moins intensément.

## 1.1. Morphologie des flots

### 1.1.1. Angle d'ouverture à la base du flot

L'angle d'ouverture $\alpha$ mesure l'angle total (projeté sur le plan du ciel) délimitant les murs de la cavité d'un flot. Il est en général calculé comme l'angle soutenu (en projection) par la cavité depuis l'origine du flot, c'est à dire,

$$ \alpha(z) = 2 \arctan \left( \frac{W(z)}{(2z)} \right) \qquad (II.1) $$

où $W(z)$ est la largeur du flot à une distance projetée $z$ de la source centrale. Des mesures d'angle plus globales ont aussi été utilisées, à savoir : un angle pondéré par la brillance, mais qui sous-estime systématiquement l'angle géométrique soutenu (Velusamy et al. 2014); et un angle pondéré par le nombre de pixels, équivalent au rapport largeur/longueur, mais qui est très sensible à la troncation du flot, pour cause de carte trop petite (Offner et al. 2011).

Même si elle évite les écueils de ces définitions globales, la mesure d'angle d'ouverture selon l'Eq. (II.1) souffre quand même de plusieurs biais qui peuvent affecter les statistiques :

1. en résolution angulaire : une faible résolution aura tendance à élargir $W(z)$, et donc à surer-timer $\alpha$.

2. hauteur $z$ à laquelle $\alpha$ est mesurée : la forme souvent parabolique des cavités fait que $\alpha$ est maximum près de la base et diminue avec l'altitude (e.g. Velusamy & Langer 1998).

3. inclinaison : les flots peu inclinés sur la ligne de visée apparaissent plus ouverts que s'ils sont vus presque par la tranche, car la distance $z$ projetée dans le plan du ciel sera plus courte que la vraie distance à la source. Pour un flot conique, on aura ainsi $\tan \alpha_{obs} = \tan \alpha_{reel} / \sin i$, où $i$ est l'angle entre l'axe du flot et la ligne de visée.

Plusieurs études récentes ont permis de s'affranchir du premier biais en mesurant $\alpha$ à haute résolution angulaire, à partir de cartes CO obtenues avec ALMA (e.g. Dutta et al. 2020) ou d'images infrarouges ou optiques de la lumière diffusée par la poussière le long des parois des cavités (e.g. Velusamy et al. 2014; Habel et al. 2021).

Les mesures Spitzer de Velusamy et al. (2014) trouvent des angles d'ouverture projetés à la base du flot qui vont de $\alpha \simeq 40^\circ$ à $\alpha \simeq 130^\circ$. L'angle augmente avec la "température bolométrique" de la protoétoile, ce qu'ils interprètent comme un effet d'âge, avec une valeur qui stagnerait autour de
$\alpha \approx 110^{\circ}$ au-delà de $\approx 8000$ ans. Cependant, cette valeur n'est pas corrigée de l'inclinaison. C'est donc une limite supérieure. De plus, leur échantillon inclut des flots à des distances différentes, et donc les mesures ne sont pas toutes au même z.

Il a fallu attendre les études plus récentes de Dutta et al. (2020) avec ALMA et de Habel et al. (2021) avec le HST dans le complexe d'Orion, pour enfin mesurer $\alpha$ de façon homogène, dans un échantillon non-biaisé de sources situées à la même distance de nous (400 pc environ).

Sur les 22 flots CO (la plupart de classe 0) mesurés par Dutta et al. (2020), l'angle projeté d'ouverture à une distance projetée $z_{proj} = 800$ au varie entre $\alpha = 25^{\circ}$ et $115^{\circ}$, avec la moitié en dessous de $65^{\circ}$.

Sur les 29 cavités identifiées en lumière diffusée par Habel et al. (2021) (principalement autour de sources de Classe 1, donc plus âgées), 50 % ont un angle projeté d'ouverture à une distance projetée $z_{proj} = 8000$ au situé dans l'intervalle $\alpha \approx 8^{\circ} - 100^{\circ}$.

Ces deux dernières études confirment que l'angle d'ouverture augmente quand on se rapproche de la source. De plus, comme leurs échantillons (non-biaisés) ont plus de chance d'être inclinés aléatoirement que celui de Velusamy et al. (2014), on peut tenter de contraindre la valeur déprojetée de $\alpha$. Par exemple, pour reproduire la faible fraction de sources ponctuelles (c'est-à-dire de sources où la ligne de visée passe par l'intérieur du cône de la cavité) dans leur échantillon, Habel et al. (2021) estiment un angle maximal déprojeté de $\alpha_{deproj} \le 70^{\circ}$. Il est donc possible que les angles d'ouverture typiques des flots évolués de Classe 1 soient significativement moins ouverts que la valeur de $110^{\circ}$ observée par Velusamy et al. (2014), non-corrigée de l'inclinaison.

### 1.1.2. Rapport d'aspect $q$

Le rapport d'aspect

$$q = \frac{L}{W} \qquad (II.2)$$

est le rapport observé entre la longueur $L$ d'une cavité de flot et sa largeur maximale $W$, donnant ainsi une mesure du degré de collimation du flot.

Les observations en CO réalisées avant les années 1990 indiquaient d'abord, pour la plupart des flots moléculaires, des rapports d'aspect faibles $q \lesssim 5$ (Lada 1985; Margulis et al. 1988). Ensuite, avec l'augmentation du pouvoir de résolution, de la sensibilité et du champ de vue des instruments d'observations, un nombre croissant de flots avec des rapports d'aspect plus élevés $q > 10$ a depuis été observé (Bachiller et al. 1990; Andre et al. 1990; Bachiller et al. 1991; Richer et al. 1992; Bally et al. 1996; Yu et al. 1999; Arce et al. 2010)). Depuis les années 2000, le nombre de flots CO observés comme étant plus longs qu'un parsec est même en augmentation (Frank et al. 2014). C'est par exemple le cas de 40 % des flots protostellaires CO observés et sélectionnés de manière non-biaisée à l'échelle du nuage moléculaire Taurus par Narayanan et al. (2012). C'est aussi le cas par exemple le cas de B5-IRS1, un flot entrainé par une protoétoile de Classe 1, et dont le rapport d'aspect actuellement observé dépasse $q > 11$.

En lien avec la physique sous-jacente, la valeur de $q$ dépend du modèle de collimation de l'éjection qui entraîne le flot (Sects. II.2 et II.3). Ce qui en fait une observable d'intérêt pour contraindre et tester ces modèles.
<figure><img src="image_10.png" /></figure>

FIGURE 5 – Cartes de canaux de vitesse en $^{12}$CO(1 − 0) du flot L1228, tel qu'observé par Arce & Sargent (2004). La vitesse centrale (dans le référentiel au repos de la source protostellaire) de chaque intervalle de largeur $\Delta v = 0.65$ km/s est montrée en haut de chaque panneau.

1.1.3. Largeur totale de flot à $|v| \simeq 1$ km/s

Les cartes de canaux de vitesse (*channel maps*) désignent des cartes d'émission (ou d'une autre quantité, comme la densité) évaluées pour le gaz avec une vitesse située dans un intervalle $v \in [v - \delta v/2; v + \delta v/2]$. La vitesse $v$ montrée est généralement celle projetée le long de la ligne de visée. Selon les cartes de canaux de vitesse observées en CO, les flancs des couches de flots les plus externes contiennent encore une fraction non-négligeable de gaz accéléré à des vitesses locales appréciables, au-dessus de $|v| \gtrsim 1.00$ km/s. C'est par exemple le cas des flots B5-IRS1 (Velusamy & Langer 1998) et RNO91 (Lee et al. 2002). La Fig. 5 montre le cas des deux lobes du flot L1228. Arce & Sargent (2004) y mesurent, de chaque côté, un lobe bien défini/distinctif dont l'émission est détectée : en dessous de $v \le -1.10$ km/s pour le lobe décalé vers le bleu; au-dessus de $v \ge 1.80$ km/s
pour le lobe décalé vers le rouge.

## 1.2. Cinématique

### 1.2.1. Distribution masse-vitesse

Une distribution masse-vitesse (MV) représente la relation entre la masse $m = m(v)$ de chaque parcelle de gaz observée et sa vitesse $v$. La masse représentée dans les observations est celle inférée à partir des émissions de raies brillantes. Il s'agit donc surtout de la masse de la matière préalablement choquée, poussée, chauffée et comprimée par les chocs d'étraves dans le flot. La vitesse $v$ représentée est le plus souvent la vitesse observée $v_{obs}$ projetée le long de la ligne de visée.

Les vitesses dans le flot décroissent lorsque le gaz de l'enveloppe est poussé. Les observations de flots moléculaires montrent que cette distribution suit le plus souvent une loi de puissance $m(v) \propto v^{-\gamma}$ de pente $\gamma > 0$ (Fig. 6i). De nouvelles observations avec des instruments plus performants ont permis d'observer un changement brusque de pente $\gamma$ entre typiquement $v \sim 6$ km/s et 20 km/s (et parfois jusqu'à 30 km/s). La pente $\gamma$ augmente ensuite de manière abrupte à hautes vitesses (Ridge & Moore 2001; Downes & Cabrit 2003; Su et al. 2004; Arce et al. 2007; Qiu et al. 2007, voir Fig. 6ii). Cette rupture dans le diagramme MV reflète la chute d'efficacité dans l'entraînement de la masse poussée, lorsque la vitesse de la couche de flot augmente (Qiu et al. 2009). En outre, la pente de la distribution MV croît avec l'âge et l'énergie du flot (Richer et al. 2000).

L'intervalle des valeurs typiques de $\gamma$ observées est passée de $\gamma \simeq 1.3-2.1$ avant les années 2000 (Rodriguez et al. 1982; Lada & Fich 1996; Ridge & Moore 2001)<fnref n="1" />, à $\gamma \simeq 1-3.5$ plus récemment (Downes & Cabrit 2003; Arce et al. 2007; Liu et al. 2017, voir Fig. 16). Le diagramme MV peut servir d'outil diagnostic de l'interaction entre l'enveloppe ambiante et le flot poussé.

### 1.2.2. "Loi de Hubble"

Un diagramme position-vitesse (PV) représente la répartition d'une quantité physique (comme la masse simulée, ou l'intensité observée) dans l'espace ($\delta, v_{obs}$), où $\delta$ est la position repérée le long d'une direction projetée qui intercepte le flot. En pratique, le diagramme PV est évalué le long d'une coupe passant soit par l'axe et l'origine du flot, soit perpendiculairement à cet axe (et à n'importe quelle distance de l'origine du flot). Le premier cas est appelé **diagramme PV longitudinal**, et le second est nommé **diagramme PV transverse** (voir Fig. 7).

Les diagrammes PV des flots ainsi observés à ce jour mettent en évidence une relation dite de "Loi de Hubble" : la vitesse $v_{obs}$ des points brillants de ce diagramme augmente avec la distance à la source. Une relation linéaire entre vitesse et distance est même notamment observée pour beaucoup de flots raisonnablement collimatés (Lada & Fich 1996; Lee et al. 2000, 2001, 2002; Arce et al. 2007, 2010; Zhang et al. 2019).

### 1.2.3. Faible superposition "bleu-rouge"

Les diagrammes PV longitudinaux observés dans certains flots exhibent, dans chaque lobe, une faible superposition (dite "rouge-bleu") entre des composantes en vitesses décalées vers le rouge ($v_{obs} \ge 0$) avec des composantes décalées vers le bleu ($v_{obs} \ge 0$). Cette faible superposition est notamment mise en évidence dans les flots plus âgés de $> 10^3 - 10^4$ ans et dans les flots des

<footnote n="1">1. Certains flots ont même pu être observées à l'époque avec des lois de puissances très abruptes, allant jusqu'à $\gamma \sim 10$.</footnote>
<figure><img src="image_11.png" /></figure>

FIGURE 6 – Deux distributions MV observées avec ALMA : (i) HH46-47, observé en CO. Les carrés rouges et bleus montrent les distributions observées dans les lobes respectivement décalés vers le rouge et vers le bleu, tiré de Zhang et al. (2016). Les carrés évidés supposent une émission en CO optiquement mince, tandis que les cercles pleins incluent une correction de l’épaisseur optique en fonction de la vitesse. Les lignes tiretées montrent des lois de puissances de pentes $-\gamma = -3.4$ (rouge) et $-\gamma = -2.7$ (bleu). (ii) chacun des deux lobes décalés vers le rouge et vers le bleu du flot MM8/F, observé en CO et SiO, et tiré de Liu et al. (2017). Les masses dérivées sont normalisées. Les symboles colorés marquent les points observés, tandis que les lignes montrent les ajustements par des lois de puissance.

FIGURE 7 – Diagrammes PV du flot HH46-47, tirés de Zhang et al. (2016). *Panneau supérieur*: carte d'émission intégrée en CO(1–0) entre 2 et 10 km/s (contours rouges) et entre –10 et –2 km/s (contours bleus). Les rectangles montrent les découpes réalisées pour obtenir les diagrammes PV des panneaux inférieurs. *Panneau central*: Diagramme PV longitudinal découpé le long de l'axe du flot avec une épaisseur de 7''. Il faut remarquer qu'ici la vitesse est montrée sur l'axe des ordonnées. *Panneau du bas*: Diagrammes PV transverses découpés perpendiculairement à l'axe du flot, à différentes altitudes, et avec une épaisseur de 6''.

<figure><img src="image_12.png" /></figure>
protoétoiles de Classe 1. La Fig. 7 nous montre un exemple très marqué de cette caractéristique, avec le flot HH46-47, qui pourrait être âgé d'au moins plusieurs milliers d'années (Hartigan et al. 2005; Zhang et al. 2016, 2019).

# 2. Modèle prédominant de flot entrainé par un vent grand angle

## 2.1. Présentation et succès du modèle

### 2.1.1. Modèle analytique de couche poussée par un vent radial (Shu et al. 1991)

Ce modèle a initialement été proposé par Shu et al. (1991) afin d'expliquer plusieurs des caractéristiques observées dans les flots moléculaires (voir Sect. II.1), dont l'accélération apparente en "loi de Hubble" (voir Sect. II.1.2.2). Ce modèle s'est imposé comme le modèle prédominant outre-atlantique, et j'en présente ici les hypothèses et prédicitions.

#### Hypothèses

Dans ce modèle, un vent radial est injecté à une vitesse uniforme $v_w$ à travers une enveloppe statique, générant un front de choc dont la forme et la dynamique dépendent des stratifications en densité de l'enveloppe et du vent, et de l'efficacité du mélange entre les deux. En supposant que le vent et la matière ambiante se mélagent instantanément après le choc (supposé très fin à cause du refroidissement radiatif efficace), la conservation de la quantité de mouvement impose que la couche choquée — que nous appellerons "couche" et qui constitue le flot CO — est aussi en expansion radiale, à l'instar du vent.

Ce modèle de flot est alors simplement caractérisé par deux fonctions angulaires P et Q qui paramétrisent:

— d'une part, la distribution angulaire du flux d'impulsion $\dot{P}_W$ du vent par unité d'angle solide:

$$ \dot{P}_w(\mu) = r^2 \rho_w(\mu) v_w^2 = \frac{\dot{M}_w v_w}{4\pi} \times P(\mu), \qquad (II.1) $$

où $\dot{M}$ est le flux de masse total du vent, r est le rayon sphérique, et $\mu = \cos \theta$ avec $\theta$ l'angle polaire (notons que puisque $v_w$ ne varie pas avec la distance, $\rho_w(\mu)$ doit varier en $1/r^2$ pour conserver le flux d'impulsion).

— d'autre part, la distribution angulaire de densité dans l'enveloppe :

$$ \rho_a(r, \theta) = \frac{a^2}{2\pi G} \times \frac{Q(\mu)}{r^2}, \qquad (II.2) $$

où a est la vitesse du son dans l'enveloppe. Le cas $Q(\mu) = 1$ correspond à la sphère isotherme singulière, solution auto-similaire isotrope de l'équilibre hydrostatique, et la fonction $Q(\mu)$ décrit la dépendance angulaire de $\rho_a$, reliée au degré d'aplatissement du nuage (par exemple par les forces magnétiques ou la rotation).
## Vitesse d'expansion du flot et loi de Hubble

Pour déterminer la vitesse d'expansion radiale $v_s(\mu)$ de la couche dans le référentiel de l'étoile, considérons d'abord l'évolution de la masse de la couche par unité d'angle solide dans la même direction, $M(\mu) = dM(\mu)/d\Omega$. La variation de $M$ par unité de temps est donnée par la somme de la masse ambiante qui entre dans la couche à la vitesse relative $v_s$ d'une part et la masse apportée par le vent à la vitesse relative $v_w - v_s$ d'autre part, soit, par unité d'angle solide (nous omettons ici la variable $\mu$ pour plus de lisibilité) :

$$ \frac{dM}{dt} = r^2 \rho_a v_s + r^2 \rho_w (v_w - v_s) \qquad (II.3) $$

Le flux de masse apporté par le vent à la couche (voir deuxième terme de l'équation ci-dessus) lui donne aussi une impulsion radiale $v_w$ par unité de masse, dans le référentiel de l'étoile. En supposant un mélange instantané dans la couche, la conservation de l'impulsion radiale par unité d'angle solide nous donne donc :

$$ \frac{d(Mv_s)}{dt} = r^2 \rho_w (v_w - v_s) v_w. \qquad (II.4) $$

En combinant ces deux équations, on obtient :

$$ M \frac{dv_s}{dt} = \frac{d(Mv_s)}{dt} - v_s \frac{dM}{dt} \qquad (II.5) $$

$$ = r^2 \rho_w \left( (v_w - v_s)^2 - \frac{\rho_a}{\rho_w} v_s^2 \right) \qquad (II.6) $$

Cette égalité nous montre que la vitesse d'expansion radiale de la couche, $v_s$, tend rapidement vers l'attracteur $v_0$ tel que la parenthèse s'annule, c'est-à-dire tel que (nous avons ici aussi omis les dépendances en $\mu$):

$$ v_0 = \frac{v_w}{1 + \left(\frac{\rho_a}{\rho_w}\right)^{\frac{1}{2}}} = \frac{v_w}{1 + \sqrt{\eta^{-1}}}, \quad (\text{avec } \eta = \frac{\rho_w}{\rho_a}). \qquad (II.7) $$

La vitesse $v_s$ va tendre vers $v_0$ car

- si $v_s < v_0$, alors $\frac{dv_s}{dt} > 0$ (voir Eq. (II.6)) et donc $v_s$ tendra à augmenter et se rapprocher de $v_0$;

- inversement, si $v_s > v_0$ alors $\frac{dv_s}{dt} < 0$, et $v_s$ tendra à diminuer pour se rapprocher aussi de $v_0$.

Puisque les champs $\rho_a$ et $\rho_w$ varient tous deux avec le rayon en $r^{-2}$ (Eqs. (II.1) et (II.2)) leur rapport est indépendant de la distance; donc la vitesse $v_s = v_0$ de la couche dans une direction donnée (Eq. (II.7)) reste constante au cours de l'expansion. Cette propriété assure qu'à chaque instant $t$, les différentes parties de la couche suivent une "loi de Hubble" commune $\dot{v}_s(\theta) = \dot{r}_s(\theta)/t$ où $t$ est l'âge du flot. C'est le premier succès de ce modèle.

En utilisant les prescriptions des Eqs. (II.1) et II.2 on obtient, aux angles où $\rho_w \ll \rho_a$

$$ v_s(\mu) \approx (av_w)^{\frac{1}{2}} \left( \frac{\dot{M}_w}{2a^3/G} \right)^{\frac{1}{2}} \left( \frac{P(\mu)}{Q(\mu)} \right)^{\frac{1}{2}}. \qquad (II.8) $$
Les observations de jets suggèrent que si le jet est seulement la partie axiale d'un vent à large ouverture angulaire, le flux de masse total du vent, $\dot{M}_w$, serait une fraction importante du taux d'accrétion sur l'étoile. En supposant que ce dernier, dans la phase protostellaire classe 0/1, est proche en moyenne du taux d'effondrement d'une sphère singulière, $\simeq a^3/G$ (Shu 1977), le deuxième terme serait alors de l'ordre de l'unité.

Avec des valeurs typiques de $v_w \simeq 100$ km/s et $a \simeq 0.2-0.6$ km/s (pour une enveloppe à 10-100 K), on obtient donc un ordre de grandeur pour la vitesse d'expansion de $v_s \simeq (av_w)^{1/2} \simeq 4-8$ km/s, en bon accord avec les flots CO observés. C'est le deuxième succès de ce modèle.

## Application à un vent X dans un coeur magnétisé (Li & Shu 1996)

Un troisième succès obtenu par ce modèle analytique, et démontré dans l'étude de Li & Shu (1996), est qu'il reproduit aussi la collimation et la pente de la distribution masse-vitesse typique-ment observée dans les flots à cette époque, avec des hypothèses physiquement motivées pour les fonctions $P(\mu)$ et $Q(\mu)$.

Notons que l'expansion du flot produit dans ce modèle reste auto-similaire dans le temps selon $r_s(\mu) = tv_s(\mu)$. Sa forme et sa collimation sont donc définies uniquement par la variation angulaire de $v_s(\mu)$, contrôlée par le rapport $\rho_w/\rho_a \propto P(\mu)/Q(\mu)$ (voir Eq. (II.7)). La distribution de masse en fonction de la vitesse est de son côté fixée par l'Eq. (II.3), qui a une dépendance différente en $P(\mu)$ et $Q(\mu)$.

Pour la stratification angulaire du vent, Li & Shu (1996) ont adopté la solution des "vents X" radiaux lancés par accélération magnétocentrifuge depuis le bord interne du disque, calculée par Shu et al. (1995). La collimation magnétique y induit un champ de densité asymptotique quasiment cylindrique de la forme:

$$ \rho_w = \rho_{w0} \frac{r_0^2}{(r \sin \theta)^2} \qquad (II.9) $$

$$ \text{donc } P(\mu) \propto \frac{1}{\sin^2 \theta}. \qquad (II.10) $$

Pour la stratification du coeur environnant, Li & Shu (1996) ont calculé les solutions $Q(\mu)$ de l'Eq. (II.2) pour l'équilibre magnétostatique d'un coeur isotherme singulier, traversé par un champ poloidal.

Li & Shu (1996) montrent qu'il est possible de dériver une expression analytique approchée de $Q(\mu)$, en fonction d'un paramètre $n$ qui augmente avec la magnétisation.

$$ Q(\mu) = \Theta(\theta) \sin^n \theta \qquad (II.11) $$

où la fonction $\Theta$ peut être approchée par une expansion en série de la forme:

$$ \Theta = (a_0 + a_2 [\sin(\theta)]^2 + \dots) + [\sin(\theta)]^n (c_0 + c_2 [\sin(\theta)]^2 + \dots) + \dots \qquad (II.12) $$

Plus $n$ augmente, plus le champ est fort par rapport à la pression thermique, et plus la densité qu'il peut supporter dans le plan equatorial augmente, avec une région de plus en plus évacuée sur l'axe polaire. Le cœur enveloppant le flot poussé parait alors de plus en plus aplati. Li & Shu (1996) préconisent la solution avec $n = 2$, car c'est la plus compatible avec le degré moyen d'aplatissement moyen observé dans les cœurs préstellaires.
Avec ces choix de $P(\mu)$ et $Q(\mu)$, ils trouvent que la forme du flot prédite est approximativement parabolique, et la distribution masse-vitesse a une pente $\gamma \simeq 2$, en cohérence avec les observations de flots disponibles à l'époque.

Cependant, ces résultats sont approximatifs et pas entièrement exacts car l'Eq. (II.4) de conservation de l'impulsion néglige l'effet du champ magnétique (qui est supposé causer la collimation du vent X et l'aplatissement du coeur).

### 2.1.2. Simulations MHD: effet du degré d'aplatissement et de la magnétisation du jet et du coeur

Afin d'étudier l'effet du champ magnétique sur le flot prédit, l'équipe de H. Shang et al. a conduit une série de simulations magnétohydrodynamiques (MHD) de l'interaction d'un vent X magnétisé avec les solutions de cœurs pré-stellaires de Li & Shu (1996), sur des durées de 500-1000 ans.

Shang et al. (2006) ont étudié l'effet de la variation du degré d'aplatissement $n$ du coeur, mais sans y inclure de magnétisation. Wang et al. (2015a) ont ensuite ajouté le champ magnétique ambiant dans le cas $n = 2$. Enfin, Shang et al. (2020) ont présenté une très large grille de modèles qui résume comment la structure et cinématique de la couche sont affectées par les paramètres $n$, $M_A$ (nombre de Mach Alfénique du jet) et $\alpha_b$ (intensité du champ magnétique du coeur par rapport à la solution de Li & Shu (1996)).

Ces simulations confirment dans l'ensemble les prédictions analytiques de Li & Shu (1996) pour la cavité de flot, à savoir une expansion homothétique et proportionnelle à l'âge $t$ ("loi de Hubble"), des formes paraboliques à la base, et des angles d'ouverture $\alpha$ et des rapports d'aspect $q$ qui couvrent tout l'intervalle de valeurs observées dans les flots (voir la Table 2 de Shang et al. (2020)). Le champ magnétique ambiant n'affecte la morphologie globale de la cavité que pour des valeurs proches de l'équipartition (cas $\alpha_b = 1$ et $n \ge 4$); une approche purement hydrodynamique reste donc acceptable dans les autres cas.

### 2.1.3. Modèle analytique simplifié avec vent-X "modifié" (Lee et al. 2001)

Le modèle original de vent X a une vitesse $v_w \simeq 100$ km/s à tous les angles. La largeur du flot à l'équateur ne cesse donc de grandir sous l'effet de ce vent rapide, pour atteindre typiquement 100 au au bout de 1000 ans (Li & Shu 1996; Shang et al. 2006; Wang et al. 2015a; Shang et al. 2020). Cela est contraire aux flots observés, dont la largeur à la base reste inférieure à 100 au même à des âges de $10^4$ ans (voir eg. l'image de HH212 sur la Fig. 2).

Pour obtenir un meilleur ajustement des observations, Lee et al. (2001) ont proposé une variante dite de "vent-X modifié". Ce modèle rajoute une dépendance angulaire **ad-hoc** à la vitesse $v_w$ du vent, qui décroît maintenant vers l'équateur selon:

$$v_w(r, \theta) = v_{w0} \cos \theta \qquad (II.13)$$

Dans cette configuration, la largeur du flot à l'équateur (où $v_w = 0$) s'annule alors que la forme près de l'axe reste inchangée.

De plus, Lee et al. (2001) adoptent une forme simplifiée pour la stratification ambiante de la solution $n = 2$ favorisée par Li & Shu (1996), mais sous une forme simplifiée qui garde seulement le premier terme du développement limité de l'Eq. (II.12) :

$$\rho_a = \rho_{a0} \left(\frac{r_0}{r}\right)^2 \sin^2 \theta. \qquad (II.14)$$
Comme nous le montrons ci-dessous, ces deux modifications par rapport à la solution exacte de Li & Shu (1996) permettent d'obtenir une expression analytique simple pour la forme de la cavité, plus pratique pour ajuster les observations.

## Forme du flot en expansion

En combinant les considérations précédentes avec l'expression de la vitesse de l'expansion de la couche $v_0$ (Eq. (II.7)) et l'expression de $\rho_w$ dans le vent X (Eq. (II.10)) on obtient la morphologie de la couche de choc pour le modèle modifié:

$$r(\theta) = tv_s(\theta) = \frac{tv_w}{1 + \sqrt{\frac{\rho_a(\theta)}{\rho_w(\theta)}}} = \frac{tv_{w0} \cos \theta}{1 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}} \sin^2 \theta} \quad (II.15)$$

où $t$ est l'âge de la cavité.

En utilisant $z = r \cos \theta$, $R = r \sin \theta$ (rayon cylindrique) et $r^2 = R^2 + z^2$, l'Eq. (II.15) peut se réécrire sous la forme d'une ellipse passant par l'origine:

$$r^2 = \frac{tv_{w0}r \cos \theta}{1 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}} \sin^2 \theta} \quad (II.16)$$

$$\Rightarrow 1 = \frac{tv_{w0}r \cos \theta}{r^2 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}} r^2 \sin^2 \theta} = \frac{tv_{w0}z}{R^2 + z^2 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}} R^2} = \frac{ztv_{w0}}{z^2 + (1 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}})R^2} \quad (II.17)$$

$$\Rightarrow z^2 - ztv_{w0} + (1 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}})R^2 = 0 \quad (II.18)$$

$$\iff \left(z - \frac{tv_{w0}}{2}\right)^2 - \left(\frac{tv_{w0}}{2}\right)^2 + (1 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}})R^2 = 0 \quad (II.19)$$

$$\iff \frac{(z - \frac{L}{2})^2}{(\frac{L}{2})^2} + \frac{R^2}{(\frac{W}{2})^2} = 1, \quad (II.20)$$

avec un grand axe $L = tv_{w0}$ égal à la longueur du flot sur l'axe $z$ (où la densité ambiante est nulle), un centre en $(R_c, z_c) = (0, \frac{L}{2})$, et un petit axe $W = L/q$ avec un rapport d'aspect $q = L/W = \sqrt{1 + \sqrt{\frac{\rho_{a0}}{\rho_{w0}}}}$ fixe dans le temps (voir Fig. 8i).

## Approximation parabolique

Lee et al. (2001) montrent qu'on attend typiquement un rapport de densité coeur/vent très grand à l'équateur, $\rho_{a0}/\rho_{w0} = \eta^{-1} \sim 10^3$. Dans toute la région du flot où $\sin \theta \gg \eta^{1/4}$, nous aurons donc $\rho_a/\rho_w \gg 1$. Dans cette région, nous avons alors l'approximation :

$$r^2 = rtv_s(\theta) = \frac{rtv_w}{1 + \sqrt{\frac{\rho_a}{\rho_w}}} \approx \frac{rtv_w}{\sqrt{\frac{\rho_a}{\rho_w}}} = \frac{tv_{w0}r \cos \theta}{\sqrt{\frac{\rho_{a0}}{\rho_{w0}}} \sin^2 \theta} \iff z = r \cos \theta = CR^2, \quad (II.21)$$

$$\text{avec } C \equiv \frac{\sqrt{\rho_{a0}/\rho_{w0}}}{\tau v_{w0}} = C(\tau)) \equiv \frac{1}{R_0} \quad (II.22)$$
En d'autres termes, partout où le gaz du coeur ambiant est beaucoup plus dense que le X-wind modifié, la forme de la cavité de flot entrainée se rapproche d'une parabole de courbure C (voir Fig. 8ii).

<figure><img src="image_13.png" /><figcaption>FIGURE 8 – Forme du flot poussé dans le modèle de X-wind modifié, à deux âges successifs $\tau_1$ puis $\tau_2$: (i) forme elliptique exacte à toutes valeurs de $\rho_a, \rho_w$ (ii) forme approchée de parabole, de courbure $C(\tau) = \sqrt{\frac{\rho_{a0}}{\rho_{w0}}}(\tau v_{w0})$, là où $\rho_a \gg \rho_w$ (flots lents $v_w \gg v_s$) (iii) carte de densité à $t = 390$ ans de la simulation HD tirée de Lee et al. 2001.</figcaption></figure>

## Cinématique du flot en expansion

Que la couche poussée soit elliptique ou parabolique, sa vitesse d'expansion à tous les angles $\theta$ suit toujours une "loi de Hubble" $v = r/t$. Par conséquent, une fois la morphologie fixée, la cinématique ($v_z, v_R$) de la couche poussée peut, dans tous les cas, être décrite simplement avec pour seul paramètre l'âge $t$ de la couche:

$$v_z = \frac{z}{t}; \quad v_R = \frac{R}{t} \qquad (II.23)$$

## Simplicité de paramétrisation

La possibilité d'ajuster les observations de flots avec des équations analytiques très simples et seulement deux paramètres libres (C et t, voir Eqs. (II.22) et (II.23)) a contribué à l'attractivité de ce modèle de vent X modifié. Par exemple, les morphologies et cinématiques observées avec ALMA dans les deux lobes du flot HH46-47 ont été ajustées par Zhang et al. 2016, 2019 avec plusieurs couches de ce modèle. D'autres exemples d'ajustement ont été présentés par Lee et al. (2000).
Des hypothèses fortes derrière le modèle

Il faut toutefois retenir que ce modèle de X-wind modifié repose sur plusieurs hypothèses fortes,
les deux premières venant du modèle de Shu et al. (1991) et les deux dernières venant en surplus :

1. l'hypothèse d'une vitesse $v_s$ radiale suppose un mélange instantané entre le vent et le gaz choqué de l'enveloppe statique;

2. la "loi de Hubble" demande que les densités $\rho_w$ et $\rho_a$ aient la même dépendance en rayon, $r^{-2}$ (afin que $v_s$ ne change pas au cours du temps),

3. une dépendance ad-hoc de $v_w$ en cos($\theta$) est nécessaire pour garantir l'annulation de la largeur de flot à l'origine, mais elle n'est pas cohérente avec le modèle original de vent X (où la vitesse est la même à tous les angles, car tout le vent provient d'un même point).

4. La forme parabolique de la cavité est obtenue seulement pour une solution (approximative) de coeur aplati avec $n = 2$; elle ne serait plus valable pour des aplatissements différents (Shang et al. 2020).

2.2. Critiques et limites : difficultés récentes du scénario

Lorsqu'il a été proposé, le scénario décrit par le modèle de X-wind original et ses variantes a rapidement eu un grand succès, surtout outre-atlantique. Mais de nouvelles critiques à leur encontre et limitations ont aussi émergé avec les observations récentes, qui n'étaient pas disponibles dans les années 1990 et début 2000.

2.2.1. Des chocs internes pas assez incurvés

L'avènement d'images infrarouges à grand champ et haute résolution a révélé sur l'axe de nombreux flots des séries de chocs très incurvés dans les raies de H₂ (e.g. image de HH212 dans McCaughrean et al. (2002).

Ces formes incurvées ne sont pas bien reproduites par les modèles de choc internes dans un vent X variable en vitesse. En effet, malgré la stratification en densité, l'onde de choc se propage à la vitesse < $v_w$ > dans toutes les directions, ce qui crée un front trop aplati sur l'axe (de courbure égale à la distance à la source, voir simulation dans Lee et al. 2001).

2.2.2. Un rapport d'aspect sous-estimé

L'exemple du flot B5-IRS1 illustre le mieux les critiques sur cette question. L'observation de la base de ce flot à haute résolution (1") avec l'interféromètre du Owens Valley Radio Observatory (OVRO) a été la première à révéler un angle d'ouverture projeté très grand, $\alpha_{proj} \ge 90^\circ$. Cela a été considéré comme un argument très fort en faveur du modèle de flot poussé par vent grand angle, et de l'idée que la rétroaction du flot sur l'enveloppe serait capable de terminer l'accrétion sur la protoétoile (Velusamy & Langer 1998).

Mais suite à l'amélioration progressive de la sensibilité des radiotélescopes single-dish, les cartes observées en CO ont pu couvrir un champ de vue de plus en plus large. Grâce à celà, les observations de (Yu et al. 1999) puis de Arce et al. (2010) ont révélé que le flot B5-IRS1 est beaucoup plus long qu'initialement pensé dans les années 1990 (voir Fig. 9). Avec des lobes CO longs d'au moins 2.2 pc, le rapport d'aspect de ce flot dépasse à présent $q > 11$.

Or, les simulations magnétohydrodynamiques de flots poussés par un vent grand angle prédisent toujours des flots beaucoup plus courts ($q < 3.5$) pour ces grands angles d'ouverture $\alpha \ge 87^\circ$, quelle
FIGURE 9 – Evolution du champ de vue dans lequel le flot CO de B5-IRS1 a été cartographié : le champ s'agrandit successivement entre Bally et al. (1996), Yu et al. (1999) et Arce et al. (2010). Les contours bleus/rouges montrent l'extension de 2 pc des deux lobes du flot obtenue dans la dernière étude. (Image tirée de Frank et al. 2014.)

<figure><img src="image_14.png" /></figure>

que soit la stratification en densité du coeur ambiant (voir Table 2 de Shang et al. 2020). Ce modèle semble donc incompatible avec la valeur réévaluée de $q > 11$ pour B5-IRS1, appellant aujourd'hui un scénario avec un modèle différent.

### 2.2.3. Des âges trop courts ou des largeurs excessives

L'agrément des angles d'ouverture avec les observations n'est pas suffisant pour tester le modèle d'entrainement par un vent X à grande ouverture. En effet, sauf si le degré d'aplatissement du coeur évolue, cet angle ne change pas avec le temps (expansion homothétique).

Par contre, l'expansion a vitesse constante prédite par ce type de modèles mène à des largeurs de cavités trop grandes aux âges typiques des classes $0/1 (10^4 - 10^5 \text{ans})$. En effet, la largeur intrinsèque du flot augmente continûment avec l'âge (de façon linéaire à $\theta$ donné) et nous montrons ci-dessous que cela mène à un désaccord avec les observations, pour les paramètres actuellement revendiqués dans les modèles de Shang et al. (2020).

L'approximation parabolique (Eq. (II.22)) nous permet de décrire analytiquement l'évolution de la largeur totale $W_z(t)$ du flot à une hauteur $z$ donnée en fonction de l'âge $t$ :

$$W_z(t) = 2R(z) = 2\sqrt{\frac{z}{C}} = 2\sqrt{zt(\rho_{w0}\sqrt{\rho_{w0}/\rho_{a0}})} = 2\sqrt{ztV_0} \quad (II.24)$$

$$\text{où } V_0 \equiv v_{w0}\sqrt{\rho_{w0}/\rho_{a0}}.$$

Or les simulations numériques MHD de Shang et al. (2020), supposées correspondre à des conditions typiques, prédissent des largeurs de cavités $W \simeq 1000$ au à $z = 800$ au de la source pour un âge très court de 500 ans. C'est déjà égal à la largeur médiane mesurée par Dutta et al. (2020) à la même altitude projetée dans 22 flots de Classe 0 (l'intervale de largeurs allant de 400 à 2500 au). D'après la loi en $\sqrt{t}$ de l'Eq. (II.24), ce type de modèles mènerait donc à des largeurs de cavités 5-15 fois trop grandes en moyenne aux âges typiques des classes $0/1 (10^4 - 10^5 \text{ans})$.

L'ajustement de la morphologie et de la cinématique des flots observés par le modèle de vent X donne aussi souvent des âges trop courts par rapport aux âges typiques des classes $0/1 (10^4 - 10^5 \text{ans})$: Par exemple, Zhang et al. (2019) ajustent le modèle avec les signatures observées dans la cinématique en CO du flot HH46-47 et trouvent des âges oscillant entre $\simeq 10^2 - 10^3 \text{ans}$.

Nous reviendrons plus en détail sur ce point à la fin du Chapitre V, où nous comparerons ces prédictions avec les résultats de nos simulations dans le cas d'un flot poussé par un jet variable.
3. Flots poussés par un jet : un modèle trop rapidement cri-
tiqué

3.1. Principe général des chocs d'étrave

<figure><img src="image_15.png" /></figure>

FIGURE 10 – Schéma descriptif général de la structure d'un choc d'étrave, dans le référentiel du choc, tiré de Raga & Cabrit (1993). Lorsque le jet rencontre un obstacle plus lent (éjection précédente, comme ici, ou milieu ambiant si l'on se place à la tête du jet) une structure en double choc se forme dans le jet (appelée surface de travail, ou Working Surface = WS). La matière choquée en surpression en est éjectée latéralement et engendre un choc d'étrave dans le milieu ambiant, qui recourbe le choc en arrière sous l'effet de sa vitesse relative (ici, $-V_j$). Une couche de mélange se développe le long de l'interface, qui peut re-remplir partiellement la cavité.

Un choc d'étrave se forme lorsqu'un jet supersonique collimaté dans une direction *z* entre en collision avec un milieu plus lent. Il peut s'agir du nuage environnant au repos (à la tête du jet), ou bien d'une portion de jet éjectée précédemment à une vitesse plus faible (le long d'un jet variable, comme dans la Fig. 10). La structure de choc qui se développe est analogue dans les deux cas, et détaillée dans la Figure 10.

Tout d'abord, l'impact produit une "surface de travail" (working surface en anglais) dans le corps du jet, qui se propage à une vitesse $V_s \le V_j$ par rapport à la source. Elle est composée de deux chocs: un choc en retour ou "Mach disk" dans lequel le jet est freiné de $V_j$ à $V_s$ et un choc en avant dans lequel le gaz situé juste devant le jet est accéléré à $V_s$.
Ensuite, le gaz de la surface de travail, comprimé entre ces deux chocs, est latéralement éjecté hors du faisceau du jet, à l'intérieur d'un cône de Mach de demi-angle d'ouverture $\alpha$ (Falle & Raga 1993, voir Fig. 10). L'interaction de cet éjecta avec l'environnement (qui a une vitesse relative négative par rapport à la surface de travail) produit une structure recourbée vers l'arrière, plus ou moins épaisse, constituée en principe de deux chocs (voir Fig. 10). D'une part, un choc en avant (forward shock) le long de sa paroi externe, où le milieu ambiant est poussé latéralement par effet "chasse-neige" : il s'agit du choc d'étrave. D'autre part, un choc en retour (reverse shock) le long de sa paroi interne où le matériau éjecté depuis la surface de travail est déposé. Entre ces deux chocs courbes se trouve une discontinuité de contact où un mélange turbulent jet-environnement peut avoir lieu. Son sillage peut re-remplir partiellement la cavité ouverte par le choc d'étrave.

Du fait de la complexité de la structure de choc de la Fig. 10, des modèles analytiques de la forme et de la cinématique n'existent que pour un jet et un milieu ambiant de densités constantes. Ils nécessitent aussi d'autres approximations supplémentaires fortes (mouvement ballistique, refroidissement fort, mélange instantané), mais leurs résultats sont en bon accord avec les simulations numériques pour des jets protostellaires denses. Nous les présentons ci-dessous.

## 3.2. Modèles de chocs d'étrave dans un milieu homogène

### 3.2.1. Choc d'étrave à la tête du jet: modèle analytique d'Ostriker et al. (2001)

FIGURE 11 – Schéma simplifié du choc d'étrave formé à la tête d'un jet stationnaire, dans le modèle d'Ostriker et al. (2001). Le référentiel est en co-mouvement avec la surface de travail (w.s. = working surface), qui avance à la vitesse $V_s$ par rapport à l'enveloppe, supposée au repos. Le mélange post-choc entre le jet et le milieu ambiant est supposé instantané.

<figure><img src="image_16.png" /></figure>

Historiquement, un premier modèle analytique de flot poussé par la tête d'un jet collimaté a d'abord été présenté par Masson & Chernin (1993). Mais nous nous intéresserons ici au modèle analytique d'Ostriker et al. (2001), plus rigoureux et complet. En effet, les prédictions de ce modèle simplifié et des simulations qui l'ont accompagnées (Lee et al. 2001) ont servi de base à ces auteurs pour affirmer les critiques et limites de ce modèle dans la reproduction des observations de flots.

Comme les modèles de la Sect. II.2, ce modèle pose plusieurs hypothèses fortes.

Premièrement, il se restreint au cas analytiquement soluble où le jet et l'enveloppe environnante ont des densités homogènes $\rho_j$ et $\rho_a$. De même, ils ne considère que le cas où le jet est stationnaire
(en densité et vitesse). De plus, ils se restreint à la limite dite "ballistique" pour calculer la forme du choc d'étrave : en dehors de la surface de travail, les forces de gradient de pression thermique sont supposées négligeables par rapport à la pression dynamique et l'inertie du gaz.

Simultanément, ce modèle suppose que le refroidissement est suffisamment efficace pour que :

- le cône de Mach (angle $\alpha$ dans la Fig. 10) soit faiblement ouvert, de sorte que l'éjection de matière se fasse perpendiculairement au faisceau du jet,

- les deux chocs qui définissent la surface de travail (zone hachurée dans la Fig. 11) se maintiennent en équilibre de pressions dynamiques. Dans l'un, le jet est freiné de $V_j$ à $V_s$ et dans l'autre, le milieu ambiant est accéléré de 0 à $V_s$. Nous pouvons alors estimer la vitesse $v_s$ d'avancée de la tête du jet par:

$$ \rho_j(v_j - v_s)^2 \approx \rho_a v_s^2, \qquad (II.1) $$

donc

$$ v_s(\eta) \approx \frac{v_j}{1 + \left(\frac{\rho_a}{\rho_j}\right)^{\frac{1}{2}}} = \frac{v_j}{1 + \sqrt{\eta^{-1}}}, \quad (\eta^{-1} = \frac{\rho_a}{\rho_j}) \qquad (II.2) $$

Où $v_j$ est la vitesse du jet et $\eta = \frac{\rho_j}{\rho_a}$ est le rapport des densités du jet et du milieu ambiant.

Troisièmement, la géométrie de ce modèle est simplifiée par rapport à la structure décrite dans la Fig. 10. On suppose que toute l'impulsion du matériau ejecté latéralement est déposée immédiatement en sortie de la zone de travail ($R = R_j$) à la base du choc d'étrave, au lieu d'être éjectée le long de sa face interne ("jet shock" dans la Fig. 10 de Raga & Cabrit (1993)). On suppose aussi que seule une portion négligeable du matériau qui entre dans le choc d'étrave en ressortira dans le sillage turbulent. Ces hypothèses sont vérifiées par les simulations numériques de Lee et al. (2001) sur des jets denses dans un milieu homogène. De cette façon, le rapport des flux d'impulsion axiaux et radiaux $\dot{P}_z(z, R)/\dot{P}_R(z, R)$ reste constant le long du choc d'étrave, et égal à la pente $dz/dR$.

À ce stade, nous pouvons obtenir la forme de la surface du choc d'étrave, en chaque rayon cylindrique $R$ et en fonction des flux d'impulsion $\dot{P}_{0z}$, $\dot{P}_{0R}$ initiaux en sortie de la zone de travail ($R = R_j$):

$$ z = z_{\text{head}} + \frac{\dot{P}_{0z}(R - R_j) - \pi \rho_a v_s^2 (R^3/3 - RR_j^2 + 2R_j^3/3)}{\dot{P}_{0R}} \qquad (II.3) $$

où $z = z_{\text{head}}$ à la tête du jet. En notant que, dans le référentiel de la tête du jet, $\dot{P}_{0z} = 0$, et en écrivant $\dot{P}_{0R} \equiv \beta c_s \pi R_j^2 \rho_a v_s$ (avec $\beta = 3.8 - 4.4$ dans les simulations de (Lee et al. 2001) la forme du choc d'étrave dans ce modèle peut se simplifier comme:

$$ z = z_{\text{head}} - \left[ \frac{1}{3} \left( \frac{R}{R_j} \right)^3 - \frac{R}{R_j} + \frac{2}{3} \right] \frac{v_s}{\beta c_s} R_j \qquad (II.4) $$

où $z = z_{\text{head}}$ à la tête du jet. Avec ces hypothèses il est aussi possible d'estimer la vitesse moyenne $\bar{v}' = \bar{v}'_z + \bar{v}'_R$ en chaque point de la couche ainsi que la vitesse $\bar{u}' = \bar{u}'_z + \bar{u}'_R$ du gaz ambiant nouvellement choqué (Fig. 13i):
$$ \bar{v}'_z = \frac{R_j^2}{R^2} v_s , \qquad \bar{v}'_R = \frac{R_j^2}{R^2} \beta c_s \qquad (II.5) $$

$$ u'_z = \frac{v_s}{1 + \frac{v_s^2(R^2 - R_j^2)^2}{(\beta c_s R_j^2)^2}} , \qquad u'_R = \frac{\frac{v_s^2(R^2 - R_j^2)}{\beta c_s R_j^2}}{1 + \frac{v_s^2(R^2 - R_j^2)^2}{(\beta c_s R_j^2)^2}} \qquad (II.6) $$

Loin de la tête du jet ($R \gg R_j$), la forme de la couche approche une loi cubique plus simple $|z| \propto R^3$ (Fig. 12) et les vitesses décroissent en loi de puissance :

$$ z \approx z_{\text{head}} - (R/R_j)^3 (v_s R_j) / (3\beta c_s) \qquad (II.7) $$

$$ u'_R \rightarrow \beta c_s (R_j/R)^2 \qquad (II.8) $$

$$ u'_z \rightarrow (\beta c_s / v_s)^2 (R_j / R)^4 v_s \qquad (II.9) $$

(II.10)

Ces prédicitions analytiques pour un milieu homogène sont confirmées par les simulations numériques de Lee et al. (2001), pour une gamme de paramètres, et avec un refroidissement purement d'origine atomique. Comme le montrent la Fig. 12 et 13, la morphologie et la cinématique sont bien reproduites.

<figure><img src="image_17.png" /><figcaption>FIGURE 12 – Simulation hydrodynamique du choc d'étrave à la tête du jet dans un milieu homogène, à un âge de 650 ans, tirée de Lee et al. (2001). La courbe en blanc montre la forme prédite par le modèle ballistique analytique de Ostriker et al. (2001) pour les mêmes paramètres.</figcaption></figure>

### 3.2.2. Le long d'un jet variable : modèle analytique de Tabone et al. (2018)

La configuration différente de Tabone et al. (2018) calcule la forme et la cinématique d'un choc d'étrave causé par la variation en vitesse d'un jet plongé dans un vent de disque uniforme, de vitesse $v_w$ et de densité $\rho_w$ (Fig. 14).

À cause de la variabilité en vitesse supersonique du jet, les ejecta récents rattrapent les ejecta plus lents émis précédemment : il y a alors formation d'une zone de travail interne au jet (IWS = internal working surface), où le jet choque contre lui-même. Contrairement à la surface de travail à la tête du jet, dont la vitesse d'avancée $v_s$ dépend de la densité ambiante, la surface de travail interne
<figure><img src="image_18.png" /><figcaption>(i) Vitesses transversale $v_R$ et axiale $v_z$ en fonction de z.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>(ii) Diagrammes Vitesse projetée-Distance projetée (le long de l'axe du flot) pour 3 inclinaisons $i$ par rapport au plan du ciel.</figcaption></figure>

FIGURE 13 – Cinématique d'un choc d'étrave à la tête d'un jet stationnaire dans un milieu homogène: prédictions analytiques d'Ostriker et al. (2001) pour les vitesses moyennes $v'_z$ et $v'_R$ (Lignes pleines) et le gaz ambiant tout juste choqué, $u'_z$, et $u'_R$ (Lignes tiretées), comparées aux simulations numériques de Lee et al. (2001) (niveau de gris). Les vitesses sont ici évaluées dans le référentiel au repos de l'enveloppe et $i$.

se propage à une vitesse proche de la vitesse moyenne du jet, $v_j$ (Raga et al. 1990), indépendamment de $\rho_a$.

Le matériau éjecté latéralement de la IWS crée un choc d'étrave dans le vent qui l'entoure (voir Fig. 14). Ce choc d'étrave se trouve à l'intérieur de la cavité principale créée par la tête du jet. Si la variabilité du jet est périodique, il se formera une série de IWS qui produira une série de chocs d'étraves "imbriqués".
FIGURE 14 – Schéma d’un choc d’étrave interne créé par un jet variable, dans la configuration de Tabone et al. (2018). Le référentiel est en co-mouvement à la vitesse moyenne du jet, $v_j$. La matière du jet est éjectée latéralement à vitesse initiale $v_0$ et forme un choc incurvé dans le vent de disque plus lent qui entoure le jet, de densité uniforme $\rho_w$ et de vitesse relative $v_j - v_w$. La distance $x$ est comptée positivement vers la source du jet.

<figure><img src="image_20.png" /></figure>

Le modèle analytique de Tabone et al. (2018) fait les mêmes hypothèses que celui d'Ostriker et al. (2001), détaillées dans la section précédente, et suppose en plus une couche mince; la vitesse obtenue est alors l'équivalent de la vitesse moyenne dans le modèle d'Ostriker et al. (2001).

Le gaz est éjecté perpendiculairement de la IWS avec une vitesse $v_0$ et un flux de masse $\dot{m}_0$. La dynamique ($v_R, v_x$) et la position cylindrique ($R_b(x)$) de chaque point du choc d'étrave sont alors caractérisés par trois lois de conservations: celle de la masse ($m$), de l'impulsion selon l'axe $x$ ($P_x$) et de l'impulsion selon l'axe $R$ ($P_R$). Ainsi, dans le repère du choc d'étrave :

$$ \dot{m} = \dot{m}_0 + \int_{R_j}^{R_b} 2\pi R' \rho_w (v_j - v_w) dR' \quad (II.11) $$

$$ \dot{P}_x = \int_{R_j}^{R_b} 2\pi R' \rho_w (v_j - v_w)^2 dR' = \dot{m}v_x \quad (II.12) $$

$$ \dot{P}_R = \dot{m}_0 v_0 = \dot{m} v_x \quad (II.13) $$

avec $\dot{P}_x$ et $\dot{P}_R$ les flux d'impulsion selon $x$ et $R$. Comme $\rho_w$ ne dépend pas de $R$, les intégrales sont directes et la forme $R(x)$ du choc d'étrave peut être déterminée en intégrant $\frac{dR}{dx} = v_R/v_x = \dot{P}_R/\dot{P}_x$ :

$$ \frac{dR}{dx} = \frac{\dot{m}_0 v_0}{\pi \rho_w (R^2 - R_j^2) (v_j - v_w)^2} \quad (II.14) $$

Ainsi, le choc d'étrave décrit une forme de loi cubique, comme dans le modèle d'Ostriker et al. (2001) :

$$ \Rightarrow x(R) = \frac{R(R^2 - 3R_j^2) + 2R_j^3}{L_0^2} \quad \text{avec } L_0 = \sqrt{\frac{3\dot{m}_0 v_0}{\pi \rho_w v_j^2}} \quad (II.15) $$
<figure><img src="image_21.png" /></figure>

(i) Cartes de densité simulées de chocs d'étrave imbriqués créés par un jet variable dans une enveloppe statique homogène ($v_w = 0$ et $\rho_w = cte$), à différents âges. La courbe cyan dans le panneau de gauche montre un ajustement aux simulations de la forme analytique donnée par l'Eq. (II.15).

<figure><img src="image_22.png" /></figure>

(ii) Variation de $v_R$ et $v_z$ en fonction de l'altitude $z$ prédites dans l'hypothèse de mélange complet (courbes bleu foncé, Eqs. (II.17) & (II.18)) et pour le gaz immédiatement post-choc (courbes magenta). En couleur: résultats de la simulation du panneau de droite, codés par la pression, et tracés séparement pour le matériau du jet et de l'enveloppe.

FIGURE 15 – Forme et cinématique du choc d'étrave créé dans un environnement homogène statique par un jet variable; comparaison entre simulations numériques et prédictions analytiques (Tabone et al. 2018).

Dans la limite des "jets étroits" $R \gg R_j$, la couche tend vers une forme cubique plus simple :

$$x \approx \frac{R^3}{L_0^2}. \qquad (II.16)$$

Dans le repère de l'observateur, la cinématique le long du choc d'étrave peut aussi facilement être
<figure><img src="image_23.png" /><figcaption>FIGURE 16 – Distributions intensité-vitesse $I_{CO}(v)$ observées dans les flots L1448, Orion A, NGC2071, L1551, et Mon R2 comparées à celles simulées pour un flot poussé par un jet, en prenant en compte la dissociation et la température d'excitation (lignes pleines). Les lignes tiretées représentent les distributions masse-vitesse $m(v)$ simulées. Figure tirée de Downes & Cabrit (2003).</figcaption></figure>

retrouvée à partir des Eqs. (II.11), (II.12) et (II.13) :

$$v_R = \frac{\dot{P}_R}{\dot{m}} = v_0 \left( 1 + \frac{3R^2 v_0}{v_j L_0^2} \right)^{-1} \qquad (II.17)$$

$$v_z = v_j - v_x = v_w + (v_j - v_w) \left( 1 + \frac{3R^2 v_0}{v_j L_0^2} \right)^{-1} \qquad (II.18)$$

Pour tester ce modèle analytique, Tabone et al. (2018) ont lancé des simulations numériques de jet variable dans le cas d'une enveloppe homogène en mouvement ou statique ($v_w = 0$). Là aussi, la forme et la cinématique simulées s'accordent avec les prédictions analytiques du modèle, comme démontré dans la Fig. 15.

## 3.3. Succès du modèle de flot entraîné par un jet

### 3.3.1. Distribution masse-vitesse $m \propto v^\gamma$

Les simulations 2D hydrodynamiques de jets variables lancés à travers un milieu homogène montrent que le modèle de flot poussé par un jet collimaté peut remarquablement bien reproduire les valeurs des pentes observées pour les distributions masse-vitesse $m(v) \propto v^{-\gamma}$ des flots CO, comme on peut le voir par exemple dans la Fig. 16, tirée de Downes & Cabrit (2003).
### 3.3.2. "Loi de Hubble" dans le diagramme position-vitesse

Dès la fin des années 1990, les simulations numériques de flots poussés par un jet collimaté parviennent à reproduire les diagrammes position-vitesse de plusieurs flots collimatés observés (Völker et al. 1999; Lee et al. 2000, 2001). En particulier, elles reproduisent avec succès des structures de type "Loi de Hubble" à vitesses faibles ($|v_{obs}| \lesssim 30$ km/s), dans ces PV : avec une vitesse observée variant linéairement avec la distance. C'est par exemple le cas pour le flot HH212 comparé aux simulations de Lee et al. (2001), ou dans les simulations de jets variables de Völker et al. (1999) (voir Fig. 17).

<figure><img src="image_24.png" /></figure>

FIGURE 17 – *A gauche:* Diagramme Position-Vitesse (PV) en CO le long de l’axe du flot pour la simulation de flot poussé par un jet variable de Völker et al. (1999). *A droite:* PV du modèle ballistique de choc d’étrave de Ostriker et al. (2001) (courbes pleines), superposé sur le PV du flot HH212 observé par Lee et al. (2000) en CO (niveaux de gris). Figure tirée de Lee et al. (2001).

### 3.3.3. Morphologie des chocs en H₂

Les simulations numériques et modèles analytiques (Völker et al. 1999; Ostriker et al. 2001; Lee et al. 2001) montrent que les chocs d’étrave poussés par un jet collimaté peuvent, beaucoup mieux qu’un vent à large ouverture, reproduire la forme des chocs très incurvés observés en H₂ sur l’axe des flots, comme dans HH212 (voir Fig. 18).

### 3.3.4. Evite les largeurs excessives

Selon les Eqs. (II.5) et II.17, la vitesse transverse $v_R$ de la couche formée par un choc d’étrave dans le modèle ballistique d’Ostriker et al. (2001) et Tabone et al. (2018) décroît avec la distance à l’axe du flot en $v_R \propto R^{-2}$. Ainsi, ce modèle prédit que l’expansion latérale décélère au fur et à mesure que le flot s’élargit. Cette propriété permet à ce modèle d’éviter d’atteindre des largeurs excessives à temps longs $t \gtrsim 10^4 – 10^5$ ans (contrairement aux modèles de flots autosimilaires poussés par un vent grand angle de la Sect. II.2).
<figure><img src="image_25.png" /></figure>

FIGURE 18 – Trois panneaux supérieurs : Images de l'émission H₂ issues de simulations 3D de flots poussés par 3 configurations de jets variables de Völker et al. (1999). Panneau inférieur: image H₂ de HH212 (nuances de gris, Zinnecker et al. 1998) + contours d'émission CO (Lee et al. 2000) + modèles ballistiques de chocs d'étrave d'Ostriker et al. (2001) en courbes noires. Tiré de Lee et al. (2001).

## 3.4. Critiques avancées contre le modèle (dans un milieu uniforme)

### 3.4.1. Collimation excessive

Selon l'Eq. (II.4) du modèle d'Ostriker et al. (2001), le rapport d'aspect $q = |z - z_{head}|/2R$ du choc d'étrave à une distance $|z - z_{head}|$ de la tête du jet est, à grande distance, égal à :

$$q = \frac{|z - z_{head}|}{2R} \approx \frac{1}{2} \left( \frac{v_s}{3\beta c_s} \right)^{1/3} \left| \frac{z - z_{head}}{R_j} \right|^{2/3} \qquad (II.19)$$

En posant $|z - z_{head}| = tv_s$ (où t est l'âge du flot) et en utilisant des valeurs de références appropriées, Ostriker et al. (2001) obtiennent :

$$q = 17.9 \left( \frac{t}{10^3 \text{ ans}} \right)^{2/3} \left( \frac{\beta c_s}{32 \text{ km/s}} \right)^{-1/3} \left( \frac{R_j}{100 \text{ au}} \right)^{-2/3} \left( \frac{v_s}{100 \text{ km/s}} \right) \qquad (II.20)$$

Par conséquent, à l'exception des âges les plus courts, un choc d'étrave dans un milieu homogène produit une couche de flot extrêmement collimatée, avec un rapport q typique supérieur à ~20 à partir de ≳ 1000 ans. Les auteurs concluent que ces rapports d'aspect sont trop élevés par rapport à ceux mesurés à l'époque, q~5 (voir Sect. II.1.1.2), faisant émerger la première critique avancée contre ce modèle. De même, ils concluent que leurs modèles ne sont pas capables de produire les
angles d'ouverture très larges ≃ 90° et les formes paraboliques observés par exemple dans le flot de B5-IRS1.

### 3.4.2. Vitesse trop faible dans les flots larges ($W \simeq 10\,000$ au)

Dans le choc d'étrave d'un jet stationnaire, la majorité de l'impulsion est communiquée au voisinage de la tête du jet, et la vitesse chute ensuite rapidement avec la distance à l'axe, en $1/R^2$ (voir Eq. (II.5)). Arce et al. (2007) a noté que cela donnait des vitesses d'expansion incompatibles avec celles observées dans les flots larges : Prenons par exemple un flot avec une largeur totale de 10 000 au, poussé par un jet stationnaire de 150 au de rayon, et en supposant une vitesse $\beta c_s = 32$ km/s d'éjection latérale du gaz en sortie de la zone d'interaction en tête de jet. La vitesse moyenne transverse correspondante prédite par l'Eq. (II.5) n'est que de ~0.03 km/s. Des valeurs aussi faibles sont incompatibles avec les observations : dans des flots pourtant encore plus larges, la vitesse des parois de cavité bien identifiées dépassent > 1 km/s.

### 3.4.3. Superposition "bleu-rouge" excessive

Comme nous pouvons le voir dans la Fig. 13, les simulations numériques de Lee et al. (2001) montrent que les vecteurs vitesse deviennent essentiellement transverses ($v_R \gg v_z$) dans les ailes du choc d'étrave. Par conséquent, les diagrammes PV longitudinaux obtenus montrent au même endroit des composantes décalées vers le rouge ($v_{obs} \ge 0$) et vers le bleu ($v_{obs} \le 0$), à toutes les inclinaisons du flot. Cependant, cette superposition "bleu-rouge" prédite n'est pas observée dans la majorité des flots. C'est la troisième critique exprimée contre le scenario d'entrainement des flots CO par un jet.

## 3.5. Face aux critiques : idée d'un jet pulsé multimillénaire dans un coeur stratifié

### 3.5.1. Cœur stratifié: des rapports d'aspect plus ressemblants

Les critiques de Ostriker et al. (2001), Arce et al. (2007), et Lee et al. (2001) étaient basées sur des modèles de flots poussés par des jets à travers une enveloppe uniforme, ou quasi-uniforme. Or, une décroissance abrupte de densité avec le rayon est prédite à l'échelle des coeurs protostellaires $\le 0.1$ pc, en raison de l'autogravité. Raga & Cabrit (1993) et Cabrit et al. (1997), une telle structure aiderait à obtenir des angles d'ouverture plus importants à la base des flots poussés par un jet. Cette prédiction a été une première fois testée par Moraghan et al. (2008) pour des jets stationnaires sur des âges < 1000 ans, mais avec une stratification en densité $\rho_a \propto [1 + z/z_c]^{-\kappa}$ ($\kappa > 0$), mais qui est différente de celle du coeur singulier isotherme $\rho_a \propto (\sin(\theta)/r)^2$ préconisé par Lee et al. (2001) (Eq. (II.14)).

Depuis, les observations indépendantes à haute définition des nuages moléculaires ont confirmé la stratification abrupte en densité des coeurs protostellaires (e.g. van Kempen et al. 2009 pour le coeur autour du flot HH46-47). En sommes, prendre en compte la stratification du coeur environnant apporterait une solution à la fois plus cohérente et plus réaliste aux problèmes de largeurs totales et/ou de rapports d'aspect, jusqu'alors posés par les modèles de flots poussés par un jet dans un coeur homogène.
### 3.5.2. Jet variable : des vitesses élevées même dans les flots âgés

L'insuffisance des vitesses des chocs d'étraves dans les flots larges (Arce et al. 2007, Sect. II.3.4.2) n'est plus une difficulté valable du modèle de flot poussé par un jet, dès lors que celui-ci est pulsé. En effet, dans ce cas, nous pourrions nous attendre à ce que la couche de flot poussée par la tête du jet soit réapprovisionnée en matière plus rapide ($v_{obs}$ plus important) par les chocs d'étraves internes générés par la variabilité du jet (Raga & Cabrit 1993; Völker et al. 1999).

### 3.5.3. Jet multimillénaire : une superposition "rouge-bleu" diminuée ?

Alors que le choc d'étrave produit à la tête d'un jet ralentit avec le temps (Sect. II.3.3.4), cette décélération transverse pourrait restreindre l'expansion latérale des ailes des chocs d'étraves imbriquées successifs dans la couche principale. Cette interaction réduirait alors la vitesse transverse des chocs d'étraves successifs, les orientant préférentiellement vers l'avant du flot. En parallèle, la décélération du choc d'étrave à temps longs diminuerait progressivement la vitesse d'expansion locale du gaz contenu dans la couche principale. Par conséquent, sur des échelles de temps longs atteignant typiquement plusieurs milliers d'années, ce double phénomène pourrait atténuer la superposition "rouge-bleu" excessive prédite et simulée à temps courts par Lee et al. (2001) (Sect. II.3.4.3). En parallèle, nous avons déjà soulevé dans la Sect. II.1.2.3 qu'une faible superposition rouge-bleu semble observée dans plusieurs flots qui seraient âgés d'au moins quelques millénaires.

Pour résumer, les critiques énoncées à l'encontre du modèle de flot poussé par un jet collimaté devraient ainsi être reconsidérées, non seulement à l'aune des nouvelles découvertes observationnelles depuis les années 1990, mais aussi de l'amélioration des capacités et performances des (super)calculateurs ainsi que des méthodes numériques. D'autant que le scénario du modèle opposé de flot poussé par un vent grand angle rencontre des difficultés récentes vis-à-vis de ces observations. En contraste, les arguments de cette Sect. II.3.5 justifient fortement l'intérêt de comparer ces observations à des simulations numériques inédites de jets variables, lancés à travers un coeur stratifié en densité, et avec un temps d'intégration allant jusqu'à $t\sim10^4$ ans.
# CHAPITRE

## III

## LE CODE MPI-AMRVAC

### Sommaire

<table><tr><td>1</td><td>Le code MPI-AMRVAC et équations résolues</td><td>39</td></tr><tr><td>1.1</td><td>Résumé des caractéristiques</td><td>39</td></tr><tr><td>1.2</td><td>Equations d'Euler et simulations hydrodynamiques</td><td>41</td></tr><tr><td>1.2.1</td><td>Des traceurs pour suivre la matière</td><td>42</td></tr><tr><td>1.3</td><td>Modularité du code : un grand potentiel d'applications</td><td>42</td></tr><tr><td>2</td><td>Méthodes, schémas de résolution numériques et machine Tycho</td><td>43</td></tr><tr><td>2.1</td><td>Approche des volumes finis et solutions numériques</td><td>43</td></tr><tr><td>2.2</td><td>Schémas numériques et limiteurs de pentes</td><td>44</td></tr><tr><td>2.3</td><td>Calculateur Tycho (MESOPSL)</td><td>47</td></tr></table>

### 1. Le code MPI-AMRVAC et équations résolues

#### 1.1. Résumé des caractéristiques

Le projet Message Passing Interface-Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC, Keppens et al. 2021) est une boîte à outil générale permettant spécifiquement de résoudre des équations aux dérivées partielles hyperboliques et paraboliques, sur $N_{dim} = 1, 2$ ou 3 dimensions spatiales, dans des coordonnées cartésiennes, cylindriques, polaires ou sphériques. Il s'appuie sur une approche Volumes Finis centrée sur les cellules d'une grille eulérienne adaptative à chaque pas de temps de la simulation. Ces applications très variées en (Astro)physique se concentrent sur les problèmes dominés par les chocs, et il est utilisé pour résoudre la dynamique des fluides et des plasmas magnétisés telles que décrites par des équations (magnéto)hydrodynamiques dans le cadre newtonien ou de la relativité restreinte (liste non-exhaustive de publications liées au code : http://amrvac.org/md_doc_publications.html).

L'Adaptive Mesh Refinement (AMR) de ce code permet d'ajuster en temps réel la maille de calcul pour en augmenter ou en diminuer la résolution, dans les zones où cela est jugé nécessaire d'après l'évolution des solutions aux équations résolues. Le maillage s'adapte ainsi continuellement autant aux phénomènes localisés (e.g. en présence de fronts de chocs) qu'aux structures à grandes échelles (e.g. configurations de tubes de flux magnétiques) par le biais d'un schéma numérique de blocs (van der Holst & Keppens 2007) : le domaine total simulé est regroupé en plusieurs sous-grilles (les blocs) de cellules, où chaque bloc a un niveau *level*. Lorsqu'une cellule d'un bloc de niveau *level* a besoin d'être raffinée d'un niveau, chaque cellule de ce bloc est raffinée en $2^{N_{dim}}$ nouvelles
<figure><img src="image_26.png" /></figure>

FIGURE 19 – Exemple d’une partie d’une grille rectangulaire adaptative en deux dimensions dans MPI-AMRVAC, raffinée ici avec des blocs allant du niveau 1 à 3 ($2^2 = 4$ cellules par niveau d’un bloc). Ici chaque bloc peut contenir $block_n x^D=2$ cellules de même niveau suivant chaque dimension. Le bloc de niveau 1 en haut à droite est d’abord raffiné une fois en 4 blocs de niveau 2. Puis, l’un deux (le bloc inférieur gauche) est ensuite raffiné en 4 blocs de niveau 3.

cellules (voir Fig. 19). Parmi toutes les techniques du code pour définir le critère de raffinement-déraffinement de chaque bloc, nous avons sélectionné celle basée sur des estimateurs d’erreurs de Löhner (Löhner 1987) de quantités physiques X choisies<fnref n="1" /> parmi le taux de refroidissement, la densité, la vitesse, l’énergie spécifique, etc....

Enfin, il est aussi construit pour pouvoir être lancé efficacement en parallèle sur plusieurs cœurs de calcul d'un ordinateur ou d'un superordinateur, grâce à une parallélisation de la structure en blocs du domaine et basée sur le standard Message Passing Interface (MPI). Continuellement mis à jour depuis deux décennies par une équipe internationale dont un de nos collaborateurs (Z. Meliani) fait partie, le code MPI-AMRVAC est gratuit et libre d'accès à l'adresse : https://github.com/amrvac/amrvac. Au début de ma thèse, j'utilise d'ailleurs une version variante adaptée et développée du code. Elle est maintenue par l'un de nos collaborateurs, Z. Meliani, qui l'a déjà validée, utilisée et publiée dans de précédentes études (Millas et al. 2017; Hervet et al. 2017; Meliani & Hervet 2017; Monceau-Baroux et al. 2015)

<footnote n="1">1. Critère où sont essentiellement quantifiées les dérivées secondes locales pondérées de X en chaque cellule.</footnote>
## 1.2. Equations d'Euler et simulations hydrodynamiques

Parmi d'autres applications et problèmes, le module d'hydrodynamique de MPI-AMRVAC en particulier permet de résoudre les équations d'Euler avec termes sources traduisant :

— la conservation de la masse :

$$ \partial_t \rho + \nabla \cdot (\rho \vec{v}) = M_p, \qquad (III.1) $$

— la conservation de l'impulsion :

$$ \partial_t (\rho \vec{v}) + \nabla \cdot (\rho \vec{u}\vec{u} + p \vec{I}) = \vec{F}_p, \qquad (III.2) $$

— la conservation de l'énergie :

$$ \partial_t (e) + \nabla \cdot ((e + p) \vec{v}) = E_p, \qquad (III.3) $$

où $\rho$ est la densité totale du fluide, $p$ est la pression thermique; $\vec{v}$ est le champ eulérien de vitesse du fluide; $M_p$, $\vec{F}_p$, et $E_p$ traduisent respectivement les termes sources des trois équations de conservations; $\vec{I}$ est le tenseur identité; et $e = \frac{p}{\gamma_{ad}-1} + \frac{\rho u^2}{2}$ est la densité totale en énergie du gaz, avec $\gamma_{ad}$ l'indice adiabatique. Dans un gaz parfait avec trois degrés de liberté, $\gamma_{ad} = 5/3$ (gaz monoatomique); tandis que dans un gaz isotherme, $\gamma_{ad} = 1$.

Le module d'hydrodynamique permet aussi de fermer le système d'équations de conservation précédent par une équation d'état du gaz $p = p(\rho, T)$ ($T$ étant la température du gaz).

La version de MPI-AMRVAC que j'utilise comporte des modules déjà écrits, validés et publiés (en partie par l'un de nos collaborateurs, Z. Meliani). Ainsi, au-delà des Eqs. (III.1)-(III.3), le code comporte plus généralement des modules permettant de résoudre les équations écrites sous la forme conservative générale :

$$ \frac{\partial \mathbf{U}}{\partial t} + \nabla \cdot \mathbf{F}(\mathbf{U}) = \mathbf{S}_{\text{phys}} \left( \mathbf{U}, \left\{ \frac{\partial \mathbf{U}}{\partial x_i} \right\}, \left\{ \frac{\partial}{\partial x_i} \frac{\partial \mathbf{U}}{\partial x_j} \right\}, \mathbf{x}, t \right) \qquad (III.4) $$

à l'aide de schémas (ainsi que de limiteurs et de méthodes) numériques déjà programmés. Dans les Eqs. (III.1)-(III.3):

$$ \vec{U} = \begin{bmatrix} \rho \\ \rho \vec{v} \\ e \end{bmatrix}, \qquad (III.5) $$

est le vecteur d'état résolu, et

$$ \vec{F} = \begin{bmatrix} \rho \vec{v} \\ \rho \vec{v} \vec{v} + p \vec{I} \\ e + p \end{bmatrix} \qquad (III.6) $$

est le flux associé à $\vec{U}$.

Les termes sources physiques $S_{phys}$ (refroidissement par collisions radiatives, champ gravitationnel, etc...) peuvent être résolus suivant différentes méthodes numériques dans le code. De plus, des modules de termes sources physiques ont déjà été écrits et validés par les développeurs (comme le refroidissement par les atomes).
### 1.2.1. Des traceurs pour suivre la matière

En plus du gaz advecté suivant les équations (III.1)-(III.3), il est possible d'ajouter des traceurs physiques ou purement numériques {$\theta_{tr} i$}$_{i=1,2,...,N_{tracer}}$ des différentes parties de ce gaz qui n'affectent pas la dynamique de ce dernier. Pour ce faire, chaque traceur rajoute simplement une équation de conservation :

$$ \partial_t(\rho\theta_{tr}) + \nabla \cdot (\rho\theta_{tr}\vec{v}) = M_p, \qquad (III.7) $$

Cette forme garantit d'une part de pouvoir traiter la quantité $D_{tr} = \rho\theta_{tr}$ de manière similaire à l'équation de conservation de la masse (Eq. (III.1)), et d'autre part d'utiliser les Eqs. (III.1) et (III.7) pour obtenir les valeurs du traceur $\theta_{tr} = \frac{D_{tr}}{\rho}$, dont le développement de l'Eq. (III.7) montre à son tour qu'il obéit simplement à l'équation d'advection pure :

$$ \partial_t \theta_{tr} + \vec{v} \cdot \nabla (\theta_{tr}) = 0, \qquad (III.8) $$

## 1.3. Modularité du code : un grand potentiel d'applications

J'ai eu l'occasion, durant la première année de ma thèse, de participer à l'école thématique "astrosim2019 : Ecole AstroSim pour l'astrophysique numérique" organisée au Centre Informatique National de L'enseignement Supérieur (CINES) de Montpellier. C'est ainsi que j'ai pu manipuler, durant deux semaines, une panoplie de codes numériques dédiés à l'Astrophysique. Alors que je ne suis qu'à quelques mois de thèse, ce qui m'a rapidement frappé est le potentiel de modularité de MPI-AMRVAC, comparé aux codes qui m'ont été présentés, tels que RAMSES, PLUTO ou MagIC (même si leurs domaines d'applications et leurs techniques numériques diffèrent, avec leurs propres forces).

Premièrement, la modularité de ce code est basée sur une stratégie de programmation tournée vers une pensée versatile, avec des options pour traiter des problèmes variés suivant un schéma de conception dit de *Bridge design pattern*. Sans entrer dans les détails, les avantages principaux de cette approche générale sont : une maintenance et une compatibilité simplifiées du code, entre différents modules physiques du code, voire avec d'autres codes compilés (e.g. gravitation, refroidissement par collisions, transfert de rayonnement, traitement de la chimie, etc...); et enfin une extension automatique des applications existantes du code par de nouvelles fonctionnalités. Par exemple, en s'appuyant sur cette modularité, Moens et al. (2022) a rajouté un module d'hydrodynamique supplémentaire pour le traitement du transfert radiatif.

Deuxièmement, cette modularité a rapidement diversifié les applications du code en Astro-physique, avec aussi bien : de l'hydrodynamique relativiste dans les rémanents de sursauts gamma (Meliani et al. 2007); l'étude de phénomènes de déplétion dans les disques circumstellaires jeunes (Robert et al. 2020); les phénomènes d'accrétion autour de supergéantes binaires X (El Mellah et al. 2017); ou encore des simulations de jets MHD relativistes (Meliani & Keppens 2009; Millas et al. 2017; Hervet et al. 2017; Meliani & Hervet 2017) et extragalactiques (Keppens et al. 2008) générés par des Noyaux Actifs de Galaxies, voire dans des quasars (Monceau-Baroux et al. 2015).

Troisièmement, l'école thématique m'a permis de constater que la barrière sélective quant à la prise en main de la programmation dans MPI-AMRVAC par un étudiant, sans forcément une formation stricte/complète en sciences numériques, est moins élevée que dans des codes comme
RAMSES, PLUTO ou MagIC. Selon moi, cette stratégie modulaire a grandement aidé dans mon apprentissage efficace et rapide du fonctionnement du code, pour son application aux problèmes que j'étudie au cours de cette thèse.

Ainsi, dans la pratique, l'intérêt de la stratégie de ce code est que :

- sauf si cela est vraiment nécessaire (e.g. pour résoudre un bug), je ne suis en général pas obligé de modifier les fichiers déjà développés à la racine du code, afin de programmer le problème physique qui m'intéresse. Je peux même facilement m'inspirer des modules déjà écrits pour écrire mes routines et des modules supplémentaires,

- la compilation de la racine du code est séparée dans une bibliothèque (library) compilée en premier avant les modules spécifiques (dits utilisateurs) que j'ai modifié pour mon problème physique. Ainsi, en général, seuls ces modules utilisateurs sont recompilés à chaque modification (économie de temps considérable à long terme);

- dans mon travail côté numérique, je peux en général me concentrer sur l'amélioration de la paramétrisation du problème physique spécifique dans les modules utilisateurs, sans affecter l'intégrité et la validité du reste du code.

# 2. Méthodes, schémas de résolution numériques et machine Tycho

## 2.1. Approche des volumes finis et solutions numériques

L'approche des Volumes Finis (VF) est largement utilisée dans la résolution des problèmes décrits par des équations aux dérivées partielles (PDEs), comme dans le cas des équations (méhéto)hydrodynamiques. Dans cette approche, les systèmes d'équations résolus, tels que les Eqs. (III.1)-(III.3), sont réécrites sous une forme conservative. Dans le cas hydrodynamique, celle-ci est donnée par les Eqs. (III.4)-(III.6) (stricto sensu, la forme conservative traite $S_{phys} = 0$).

Du point de vue numérique, les Eqs. (III.1)-(III.3) sont discrétisées dans l'espace et dans le temps, tel que :

$$ \vec{U}^{n+1} = \vec{U}^n + \Delta t^n \left[ \sum_{idim=1}^{idim=ndim} \frac{\vec{F}_{idim}^{i+\frac{1}{2}} - \vec{F}_{idim}^{i-\frac{1}{2}}}{\Delta x_{idim}} \right], \quad (III.1) $$

où $U^n = (\rho^n, (\rho\vec{v})^n, e^n)^T$ et $U^{n+1} = (\rho^{n+1}, (\rho\vec{v})^{n+1}, e^{n+1})^T$ sont les variables dites conservatives, évaluées respectivement aux temps $t^n$ et $t^{n+1} = t^n + \Delta t^n$, avec $\Delta t^n$ le pas de temps. Les flux moyennés (dans le temps) $\vec{F}_{idim}^{i+\frac{1}{2}}$ et évalués pour l'intervalle de temps $[t^n, t^{n+1}]$ au niveau de l'interface située entre une cellule d'indices i et ses voisins i + 1 dans la direction idim, sont calculées à partir des solutions aux problèmes de Riemann<fnref n="1" />. Cette approche est illustrée sur la Fig. 20.

En HD, le problème de Riemann est décrit par une structure à 3 ondes caractéristiques. Ces trois valeurs propres correspondent aux ondes de choc (propagées à l'avant) et de raréfaction (se déplaçant à l'arrière du choc), avec entre les deux une onde entropique (discontinuité de contact). Ces ondes sont définies dans la direction idim par :

<footnote n="1">1. Un problème de Riemann est un problème à valeur initiale pour une équation de dynamique des fluides conservative. Elle implique une distribution discontinue des variables conservatives. Ce problème amène au calcul des caractéristiques ou des valeurs propres correspondant à différentes ondes solutions du problème.</footnote>
<figure><img src="image_27.png" /><figcaption>FIGURE 20 – Schématisation de l’approche Volume Finis, avec des variables initialisées au centre de chaque cellule et des flux calculés au centre des faces</figcaption></figure>

$$ 
\begin{aligned}
\lambda_{1,3}^{\text{idim}} &= v^{\text{idim}} \mp c_s \\
\lambda_2^{\text{idim}} &= v^{\text{idim}}
\end{aligned}
\qquad (\text{III.2})
$$

où $c_s = \sqrt{\left(\frac{\partial p}{\partial \rho}\right)}$ est la vitesse locale du son. Nous pouvons ici noter que les ondes HD caractéristiques, décrites sur la Fig. 21, sont indépendantes de la direction.

## 2.2. Schémas numériques et limiteurs de pentes

Les schémas ou solveurs de Riemann sont un outil fondamental dans le développement de l'approche Volumes Finis. Ils sont basés sur une simplification des équations hyperboliques. En outre, à chaque pas de temps, ils simplifient l'état physique en des valeurs constantes par morceaux, avec des sauts de discontinuités au niveau de certains ou de tous les vecteurs propres (avec leurs valeurs propres associées) qui caractérisent les équations hyperboliques. Dans le cas des équations HD, chaque onde est associée à une valeur propre (de vitesses caractéristiques Eq. (III.2)) : $\lambda_{1,3}$ sont associés aux ondes de chocs et de raréfaction, tandis que $\lambda_2$ est associée à une discontinuité de contact (voir Fig. 21).

Il existe plusieurs types de solveurs de Riemann (voir e.g. LeVeque et al. 1998):
<figure><img src="image_28.png" /><figcaption>FIGURE 21 – *Eventail de Riemann* : Solutions d'ondes caractéristiques au problème de Riemann hydrodynamique (*Riemann fan*) en 1 dimension.</figcaption></figure>

— Un solveur de Riemann exact requiert une méthode itérative et n’est donc pas pratique pour un code HD. Cependant, les résultats du problème de Riemann exact sont utilisés comme solutions de références pour vérifier la précision numérique et les performances des solveurs de Riemann approchés.

— Un solveur de Riemann linéarisé tel que le solveur de Roe (Roe 1981) requiert une décomposition en caractéristiques des vecteurs propres à gauche et à droite de chaque interface de celulle. Ce solveur est complexe et très consommateur en temps de calcul dans les cas HD et MHD. En outre, le solveur de Roe peut mener à des solutions avec pressions négatives non-physiques (Einfeldt et al. 1991).

— Un solveur de Riemann basé sur des estimations, tel que le solveur HLL (Harten-Lax-van Leer), introduit une estimation de la vitesse du son ainsi qu'une solution moyennée sur l'ensemble de l'éventail de Riemann (Fig. 21).

Durant cette thèse, parmi tous les solveurs disponibles dans MPI-AMRVAC, nous nous sommes concentrés sur l'usage des solveurs de Riemann de types basés sur des estimations, les plus rapides parmi les trois types de solveurs précédemment décrits. Dans cette catégorie, MPI-AMRVAC possède ainsi quatre solveurs principaux parmi les plus utilisés :

— le solveur de Riemann TVDLF (Total Variation Diminishing Lax-Friedrich Rusanov) (e.g. Tóth & Odstrčil 1996; Yee 1994) est construit en supposant un état moyen donné par l'onde caractéristique la plus rapide (Fig. 22). En effet, le flux à l'interface $i + 1/2$ est donné par :

$$F^{i+1/2} = \frac{1}{2} [F(U_{L,i+1/2}) + F(U_{R,i+1/2}) - S_{\max}(U_{R,i+1/2} - U_{L,i+1/2})], \quad (III.3)$$

où la vitesse de propagation la plus rapide est :

$$S_{\max} = \max(S_{L,\max}, S_{R,\max}),$$

avec $S_{L,\max}$, $S_{R,\max}$ les valeurs propres associées à l'onde se propageant le plus rapidement, respectivement à gauche et à droite de l'interface $i + 1/2$. Elles sont définies comme :

$$S_{Lmax} = \max_{(1 \le k \le 3)} (|\lambda_{k,L}|), \quad S_{Rmax} = \max_{(1 \le k \le 3)} (|\lambda_{k,R}|). \quad (III.4)$$

— Le solveur HLL (Harten-Lax-van Leer) proposé par Harten et al. (1983) est construit en supposant un état intermédiaire moyenné entre l'onde la plus lente et la plus rapide. Le
<figure><img src="image_29.png" /><figcaption>(i) Ondes résolues par les solveurs HLL,HLLC et HLLCD(*: lorsque la sous-région traitée applique HLLC)</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>(ii) Ondes résolues par les solveurs TVDLF et HLLCD(**: lorsque la sous-région traitée applique TVDLF). $c_{max}$: vitesse de la caractéristique la plus rapide</figcaption></figure>

FIGURE 22 – Solutions au problème de Riemann hydrodynamique (Riemann fan) avec en plus les sous-structures résolues par les solveurs HLL,HLLC,TVDLF, et HLLCD

solveur HLL traite ainsi uniquement une partie des ondes caractéristiques de l'éventail de Riemann (Fig. 22i): les ondes de choc (à l'avant) et de raréfaction (se déplaçant à l'arrière du choc).

$$F^{i+1/2} = \begin{cases} F(U_{L,i+1/2}) & S_{L,i+1/2} > 0 \\ F(HLL, i + 1/2) & S_{L,i+1/2} \le 0 \le S_{R,i+1/2} \\ F(U_{R,i+1/2}) & S_{R,i+1/2} < 0 \end{cases} \quad (III.5)$$

où:

$$F(\text{HLL}, i + 1/2) = \frac{S_{R,i+1/2} U_{R,i+1/2} - S_{L,i+1/2} U_{L,i+1/2} + F(U_{L,i+1/2}) - F(U_{R,i+1/2})}{S_{R,i+1/2} - S_{L,i+1/2}} \quad (III.6)$$

Les vitesses des ondes les plus externes de l'éventail de Riemann, $S_{L,i+1/2}$ et $S_{R,i+1/2}$, sont estimées en utilisant les états à gauche et à droite :

$$\begin{aligned} S_L &= \min_{(1 \le k \le 3)} (\lambda_{k,L}, \lambda_{k,R}) \\ S_R &= \max_{(1 \le k \le 3)} (\lambda_{k,L}, \lambda_{k,R}) \end{aligned} \quad (III.7)$$

— Le solveur HLLC (Harten-Lax-van Leer-Contact) proposé par Toro (1997) est un solveur HLL
à deux états. Il introduit en effet des sous-structures associées à la discontinuité de contact dans le sous-état lent du solveur HLL (qui ne tenait pas compte de cette discontinuité).

— Le solveur de Harten-Lax-van Leer-Contact-Diffused (HLLCD) (Meliani et al. 2008) : HLLC est partout appliqué, avec cependant un basculement sur TVDLF à l'intérieur de blocs où suivant une direction donnée, un changement de signe a lieu entre les flux numériques de l'interface gauche et droite d'une cellule. Cela permet de réduire les bruits et instabilités numériques potentielles pouvant apparaître avec HLLC.

Dans les solveurs de Riemann présentés ci-dessus, les estimations numériques des flux à l'interface $F^{i+1/2}$ requièrent les valeurs des variables conservées à gauche $U_{L,i+1/2}$ et à droite $U_{R,i+1/2}$ de l'interface de la cellule située en $x_{i+1/2}$. Ces variables sont reconstruites à partir des valeurs $U$ évaluées au centre des cellules adjacentes à l'interface (voire sur davantage de cellules voisines). Ces reconstructions spatiales peuvent être réalisées à l'aide d'un limiteur de pente. Il existe une grande variété de limiteurs de pentes, tels que MINMOD (Roe 1986, limiteur de type Total Variation Diminishing symétrique et classique d'ordre 2), la Piecewise Parabolic Method (PPM, Colella & Woodward 1984, limiteur TVD d'ordre 3 basé sur un stencil à quatre points) ou Koren (Koren 1993, limiteur de type TVD, asymétrique et spatialement précis à l'ordre 3).

La résolution numérique des équations aux dérivées partielles (III.1)-(III.3) requiert aussi une précision d'ordre élevé en temps. Ceci est par exemple permis par l'usage d'un schéma d'ordre deux de type prédicteur-correcteur, même si des méthodes d'ordres plus élevés existent dans le code (e.g. le schéma de Runge-Kutta, très robuste pour préserver la stabilité des solutions).

## 2.3. Calculateur Tycho (MESOPSL)

Grâce à la parallélisation déjà faite sur la racine du code MPI-AMRVAC, nos simulations ont pu être lancées en parallèle sur les nombreux processeurs de la gruppe de calcul interne Tycho de l'Observatoire de Paris. Nous avons décidé d'y employer la version Intel® de la librairie OpenMPI (²) qui permet spécifiquement de tirer parti des performances des compilateurs des processeurs Intel® de Tycho. À titre de comparaison, les compilateurs GNU retrouvés sur les ordinateurs de bureau classiques sont a priori moins performants, mais ont l'avantage de la portabilité : ces compilateurs sont retrouvés sur tous les systèmes GNU/Linux, avec de nombreux codes fournis avec des fichiers Makefile écrits pour ces compilateurs.

Pour toutes les raisons ci-dessus, et sous la tutelle de Z. Meliani en tant qu'expert de ce code, MPI-AMRVAC est vite apparu comme un code plus qu'approprié pour ma thèse, pour ma formation personnelle en sciences numériques, et pour la modélisation des problèmes physiques des Chapitres IV, V, et VI.

<footnote n="2">2. voir https://docs.hpc.shef.ac.uk/en/latest/sharc/software/parallel/openmpi-intel.html</footnote>
<figure><img src="image_31.png" /></figure>
# CHAPITRE

## IV

# APPLICATION DU CODE MPI-AMRVAC POUR LA MODÉLISATION DE JETS PROTOSTELLAIRES VARIABLES DANS UN CŒUR STRATIFIÉ

## Sommaire

1 Modules utilisateurs 50

1.1 Travail de modélisation durant la thèse 50

1.2 Liste des modules utilisateurs 51

2 Modélisation des équations d'état hydrodynamiques 52

2.1 Équations d'état 52

2.2 Composition du gaz 52

2.2.1 Formules générales 53

2.2.2 Gaz parfait ionisé 54

2.2.3 Gaz parfait doublement ionisé 54

2.2.4 Gaz parfait complètement moléculaire 54

2.2.5 Gaz parfait complètement atomique et neutre 55

3 Enveloppe environnante modélisée 55

3.1 Liste des paramètres libres 55

3.2 Conditions initiales 55

3.2.1 Champs de densité et de vitesse du cœur 55

3.3 Conditions aux bords de l'enveloppe 56

3.3.1 Cellules fantômes 56

3.3.2 Conditions aux bords en $R = 0, R = R_{max}, z = z_{max}, z = 0$ 56

3.3.3 Conditions aux bords pour le cœur en $z = 0$ 58

4 Modélisation du jet pulsé 61

4.1 Liste des paramètres libres 61

4.2 Valve d'injection: domaine et ouverture angulaire 61

4.3 Champs de densité et de vitesse initiaux dans la valve d'injection 64

4.3.1 Pour une valve rectangulaire 64

4.3.2 Pour une valve trapézoïdale 65

4.4 Variabilité du jet : un paramètre physique encore plus ou moins libre 67
<table><tbody><tr><td>4.4.1</td><td>Profil de variabilité sinusoïdal</td><td>67</td></tr><tr><td>4.4.2</td><td>Profil de variabilité en dent de scie (sawtooth)</td><td>68</td></tr><tr><td><b>5</b></td><td><b>Termes sources</b></td><td><b>68</b></td></tr><tr><td>5.1</td><td>Force de contrebalancement <b>F&#773;<sub>p</sub></b> et terme source <b>E<sub>p</sub></b></td><td>69</td></tr><tr><td>5.1.1</td><td>Evaluation de <b>F&#773;<sub>p</sub></b> : méthode numérique</td><td>69</td></tr><tr><td>5.1.2</td><td>Evaluation de <b>F&#773;<sub>p</sub></b> : méthode analytique</td><td>70</td></tr><tr><td>5.1.3</td><td>Sens physique de l'énergie de contrebalancement <b>E<sub>p</sub>(t)</b></td><td>70</td></tr><tr><td>5.1.4</td><td>Évaluation de <b>E<sub>p</sub></b></td><td>72</td></tr><tr><td>5.2</td><td>Refroidissement -<b>n<sup>2</sup><sub>H</sub>Λ(T)</b> par désexcitations collisionnelles</td><td>72</td></tr><tr><td><b>6</b></td><td><b>Tests de configuration pour notre problème physique</b></td><td><b>73</b></td></tr><tr><td>6.1</td><td>Choix des schémas numériques et limiteurs de pente</td><td>73</td></tr><tr><td>6.1.1</td><td>Critères de sélection</td><td>74</td></tr><tr><td>6.1.2</td><td>Résultats</td><td>74</td></tr><tr><td>6.1.3</td><td>Schéma d'intégration en temps</td><td>75</td></tr><tr><td>6.2</td><td>Résolution spatiale : influence du niveau de raffinement par AMR</td><td>75</td></tr><tr><td>6.2.1</td><td>Critères de sélection</td><td>76</td></tr><tr><td>6.2.2</td><td>Résultats</td><td>76</td></tr><tr><td>6.3</td><td>Paramètres de parallélisation avec MPI</td><td>77</td></tr><tr><td>6.3.1</td><td>Test d'extensibilité forte et loi d'Amdhal</td><td>77</td></tr><tr><td>6.3.2</td><td>Méthode</td><td>78</td></tr><tr><td>6.3.3</td><td>Résultats</td><td>78</td></tr></tbody></table>

Dans ce chapitre, je présente le travail d'adaptation des modules du code MPI-AMRVAC que j'ai réalisé en collaboration avec Zakaria Meliani en vue d'effectuer la modélisation d'un jet variable lancé à travers une enveloppe environnante stratifiée. En particulier, je montrerai comment l'enveloppe entourant le jet est construite, le domaine d'injection de ce dernier, le maillage de la simulation, les conditions aux bords adoptées, les conditions initiales, ainsi que l'inclusion des divers termes sources (refroidissement atomique, gravitation effective, etc...) pour un gaz atomique faiblement ionisé.

Ensuite, je montrerai une série de tests numériques m'ayant permis de sélectionner la configuration de base pour le lancement des simulations des Chapitres V et VI en termes de : schémas et limiteurs de pente numériques, paramètres d'AMR, paramètres de parallélisation de l'exécution du programme.

# 1. Modules utilisateurs

## 1.1. Travail de modélisation durant la thèse

J'ai débuté ma thèse avec une version particulière de MPI-AMRVAC, maintenue par Z. Meliani, ainsi que six modules compilables fournis, dits "utilisateurs". La compilation de ces modules externes intervient après celle de la racine du code, et ils servent à poser et caractériser les problèmes physiques qui m'intéressent en utilisant MPI-AMRVAC. Les modules utilisateurs fournis au début de ma thèse ont déjà été utilisés dans des études précédentes sur d'autres problèmes physiques (par
exemple dans les jets d'AGN, Hervet et al. 2017, Millas et al. 2017, Marcowith et al. 2020). Cependant, en tant que tels, ces modules ne permettaient pas de résoudre spécifiquement les situations physiques auxquelles je me suis intéressé. Par conséquent, la première étape de ma thèse a été de me familiariser avec le code. En collaboration avec Z. Meliani, j'ai consacré une bonne partie de mon travail de thèse à comprendre, adapter, et vérifier la version du code qui m'a été fournie, afin de répondre aux objectifs fixés tout au long de la thèse. Tout en analysant le code, j'ai aussi apporté mes propres modifications et implémentations à la version initiale : certaines modifications ont été moins importantes et en marge de la version de base ; d'autres, bien plus importantes et volumineuses.

## 1.2. Liste des modules utilisateurs

Voici la liste des six modules utilisateurs sur lequel Z. Meliani et moi avons travaillés :

1. `mod_usr`, le module général utilisateur. Ce centre névralgique et polyvalent des modules utilisateurs effectue la connexion avec la bibliothèque compilée à la racine du code. Il définit aussi des routines permettant d'évaluer l'ajout de termes sources (champ gravitationnel, etc...), le pas de temps qui en découle, des procédures d'écriture/de lecture de fichiers durant la simulation, ou encore des opérations sur la grille de calcul (e.g. manipulation de l'AMR) et/ou sur les champs manipulés à différentes étapes de l'exécution du code (e.g. avant/après évaluation des conditions aux bords, avant/après résolutions des champs à $t + \Delta t$, etc...).

2. `mod_obj_usr_unit`, le module de normalisation. Un grand nombre de modules programmés par les développeurs du code travaillent avec des variables physiques ou des constantes fondamentales dont la valeur numérique est normalisée par des valeurs de références, identiques et inchangées à travers tout le code (e.g. température normalisée par une même valeur de référence). Le but principal est de minimiser les erreurs produites par les opérations sur la virgule flottante, quand les valeurs numériques s'étendent sur plusieurs ordres de grandeurs, soit trop grandes (e.g. $1.2345 \times 10^{10} \text{ cm}^{-3}$), soit trop petites (e.g. $1.2345 \times 10^{-31} \text{ g cm}^{-3}$), et dont les combinaisons peuvent créer des erreurs d'arrondis. L'utilisateur du code est libre de choisir les valeurs de normalisation (définissant automatiquement le facteur inverse permettant de redonner les valeurs "physiques" originales). Avec Z. Meliani, j'ai choisi de relier les valeurs de normalisation entre elles, suivant les mêmes équations physiques que dans le problème physique (e.g. Eqs. (IV.1)-(IV.3)), afin que ces normalisations en reflètent les échelles caractéristiques.

3. `mod_obj_cla jet`, le module du jet. Il regroupe les paramètres, variables, champs, et routines principales associées au traitement de la zone d'injection du jet (paramétrisation, initialisation, mise à jour à chaque pas de temps).

4. `mod_obj_ism`, le module de l'enveloppe. Comme pour le module du jet, il regroupe les variables, paramètres, champs, et routines principales associées à l'enveloppe environnante.

5. `mod_obj Mat`, le module utilitaire et de mathématique. Ce module particulier implémente à la fois d'une part des conditions aux bords additionnelles à celles de la racine du code dans la version reçue au début de la thèse, et d'autre part des routines mathématiques pour réaliser des opérations et appliquer des formules mathématiques susceptibles d'être utilisées plus d'une fois dans les autres modules utilisateurs (calcul de l'angle polaire entre les cellules et l'axe du jet; profil de la variabilité de la vitesse ET de la densité dans le jet; etc...)

6. `mod_obj_global_parameters`, le module des paramètres globaux utilisateurs. Il regroupe les variables, structures et types dérivés partagés, connus par l'ensemble des 5 autres modules
utilisateurs (mais pas accessibles par les modules de la racine du code).

## 2. Modélisation des équations d'état hydrodynamiques

### 2.1. Équations d'état

Pour compléter les Eqs. (III.1)–(III.3), il nous faut une relation de fermeture : une équation d'état reliant la pression $p$ du gaz, sa densité $\rho$, son poids moléculaire moyen $\mu$, et sa température $T$. Je me suis aussi occupé de vérifier et d'adapter le code moi-même, afin que ces équations soient respectées à toutes les étapes du code : relation entre les paramètres et leurs valeurs de normalisations, initialisation, mise en place/à jour des conditions aux bordes, relation entre les champs de $p$, $\rho$, $T$ et $\mu$ à chaque pas de temps $t > 0$. Que ce soit pour le gaz dans l'enveloppe ou dans la zone d'injection du jet, il a donc fallu adapter les routines des modules utilisateurs aux équations d'un gaz parfait :

1. sur la relation entre $p$, $\rho$ et $T$:

$$p = \frac{\rho}{\mu m_H} k_B T \qquad (IV.1)$$

où $m_H$ est la masse de l'atome d'hydrogène; $k_B$ est la constante de Boltzmann; $T$ est la température du gaz; et $\mu$ est le poids moléculaire moyen reliant $\rho$ au nombre volumique total $n_{tot}$ de particules, et dont nous supposons qu'elles contribuent ici toutes à la pression thermique:

$$\mu \stackrel{\text{def}}{=} \frac{\rho}{n_{tot} m_H} \qquad (IV.2)$$

2. sur les relations entre la vitesse du son $c_{s,a}$ du gaz parfait dans le cœur et sa pression thermique $p$ ou sa température $T$ :

$$c_{s,a}^2 = \frac{\gamma p}{\rho} = \frac{\gamma k_B T}{\mu} \qquad (IV.3)$$

où $\gamma = C_P/C_V$ est le rapport des capacités calorifiques spécifiques. Le poids $\mu$ dépend de la composition chimique du gaz (Sect. IV.2.2).

### 2.2. Composition du gaz

Dans la version de MPI-AMRVAC reçue au début de ma thèse, le module traitant du cas hydrodynamique, avec refroidissement par désexcitations collisionnelles des atomes, n'étaient écrits que dans l'hypothèse d'un gaz composé uniquement d'hydrogène atomique et sans hélium. En effet, ils considéraient uniquement le cas $\mu = 1$ et où le poids moyen $\bar{m}$ en équivalent d'hydrogène du gaz vaut aussi $\bar{m} = 1$ (voir Sect. IV.2.2.1 pour la définition de $\bar{m}$). Ces quantités sont impliquées dans les relations entre la densité, la température, la pression et la vitesse du son. Sans entrer dans les détails, il a ainsi fallu rajouter un champ supplémentaire de poids moléculaire moyen $\mu$ ainsi qu'un paramètre $\bar{m}$, tous deux manquants dans la version originale. Bien que l'une des quantités physiques reliées à la pression et à la densité soit le nombre volumique total $n_{tot}$ de particules :

$$n_{tot} = n(H) + n(H^+) + n(H_2) + \underbrace{n(He) + n(He^+) + n(He^{2+})}_{\equiv n_{He}} + n(e^-), \qquad (IV.4)$$
les études précédentes sur les jets et flots protostellaires utilisent et montrent plutôt plus souvent
le nombre volumique $n_H$ en nucléons d'hydrogène du gaz :

$$n_H \stackrel{\text{def}}{=} n(\text{H}) + n(\text{H}^+) + 2n(\text{H}_2) \qquad (\text{IV.5})$$

Cette quantité pratique est par exemple utilisée dans l'évaluation des abondances relatives à l'hydrogène des espèces chimiques qui peuplent un gaz (par observations ou simulations numériques). Dans le module traitant du refroidissement, il nous a fallu ajouter la prise en compte de $\mu$ dans $n_H$. Je me suis occupé de tester les modules modifiés et même d'en corriger des bugs affectant l'énergie. J'ai aussi écrit une nouvelle routine permettant d'évaluer la valeur de $\mu$ et de $\bar{m}$ en fonction de la composition supposée pour le gaz.

L'évaluation de la pression $p$ à partir de l'énergie totale $e$ d'un gaz parfait suit la relation $p = \rho[e - 0.5v^2]/(\gamma - 1)$ et n'implique pas d'évaluer ni $n_H$, ni $\mu$, ni $\bar{m}$. Par contre, la température $T$ à chaque pas de temps est reliée à la densité $\rho$, à la pression thermique $p$ et au poids moléculaire moyen $\mu$ (qu'il a fallu inclure) :

$$T = \frac{m_H \mu p}{k_B \rho} \qquad (\text{IV.6})$$

### 2.2.1. Formules générales

De manière générale, pour un gaz composé d'hydrogène atomique et moléculaire neutre et ionisé (à savoir H, H+, H₂), ainsi que d'hélium 4 neutre et ionisé (He, He+, He²⁺), et enfin d'électrons (e⁻), nous avons les décompositions et définitions suivantes :

* le nombre volumique $n(X_i)$ de particules de l'espèce $X_i$ dans le gaz en cm⁻³.

* la fraction massique $x(He)$ d'hélium 4 par rapport à l'hydrogène :

$$x(He) \stackrel{\text{def}}{=} \frac{n(He) + n(He^+) + n(He^{2+})}{n_H} = \frac{n_{He}}{n_H} \qquad (\text{IV.7})$$

* le poids $\bar{m}$ en équivalent d'hydrogène du gaz :

$$\bar{m} = \frac{\rho}{n_H m_H} = \frac{[n(\text{H}) + n(\text{H}^+) + 2n(\text{H}_2) + 4n(\text{He}) + 4n(\text{He}^+) + 4n(\text{He}^{2+})] m_H}{n_H m_H} \qquad (\text{IV.8})$$

Ici, puisqu'un noyau d'hydrogène est ~1836 fois plus lourd qu'un électron, nous supposons pouvoir négliger l'apport de la densité massique des électrons à la densité $\rho$, ce qui permet de simplifier l'expression de $\rho$ et de ne pas y inclure la masse de l'électron. Cependant, le nombre volumique d'électrons $n(e⁻)$ du gaz n'est pas négligé dans le bilan Eq. (IV.4) du nombre volumique total $n_{tot}$, servant derrière à évaluer $\mu$ :

$$\bar{m} = \frac{\rho}{n_H m_H} = \mu \frac{n_{tot}}{n_H} \qquad (\text{IV.9})$$

Les valeurs de $\bar{m}$ et de $\mu$ dépendent des abondances chimiques des espèces précédentes dans le gaz. Même si dans les Chapitres V et VI, je suppose que le gaz maintient une fraction d'hélium 4 constante partout et dans le temps, valant $x(He) = 0.1$, j'ai programmé le code pour prendre un paramètre $x(He)$ quelconque. En outre, j'ai en pratique implémenté dans la partie utilisateur du code MPI-AMRVAC plusieurs compositions possibles du gaz (ci-après).
### 2.2.2. Gaz parfait ionisé

Supposons que le gaz n'est composé que d'hydrogène sous forme ionisée H⁺, d'hélium 4 simplement ionisé He⁺, et des électrons arrachés à ces deux espèces (de sorte que la charge totale reste nulle) :

$$n(H^+) = n_H; n(He^+) = x(He)n_H; n(H), n(H_2), n(He), n(He^{2+}) \ll n_H \quad (IV.10)$$

$$n(e^-) = n(H^+) + n(He^+) = (1 + x(He))n_H \quad (IV.11)$$

Par conséquent, selon les équations de la Sect. IV.2.2.1 :

$$\bar{m} = \frac{[0 + 0 + 0 + 0 + 4 \times x(He)n_H + 0] m_H}{n_H m_H} = 1 + 4x(He) \quad (IV.12)$$

$$n_{tot} = 0 + n_H + 0 + 0 + x(He)n_H + 0 + (1 + x(He))n_H = (2 + 2x(He))n_H \quad (IV.13)$$

$$\implies \mu(\vec{r}, t) = \frac{\bar{m}}{\frac{n_{tot}}{n_H}} = \frac{1 + 4x(He)}{2 + 2x(He)} \quad (IV.14)$$

### 2.2.3. Gaz parfait doublement ionisé

Supposons que le gaz n'est composé d'hydrogène sous forme ionisée H⁺, d'hélium 4 cette fois-ci doublement ionisé He²⁺, et des électrons arrachés à ces deux espèces (de sorte que la charge totale reste nulle) :

$$n(H^+) = n_H; n(He^{2+}) = x(He)n_H; n(H), n(H_2), n(He), n(He^+) \ll n_H \quad (IV.15)$$

$$n(e^-) = n(H^+) + 2n(He^{2+}) = (1 + 2x(He))n_H \quad (IV.16)$$

Par conséquent, selon les équations de la Sect. IV.2.2.1 :

$$\bar{m} = 1 + 4x(He), \quad n_{tot} = (2 + 3x(He))n_H \quad \implies \mu(\vec{r}, t) = \frac{1 + 4x(He)}{2 + 3x(He)} \quad (IV.17)$$

### 2.2.4. Gaz parfait complètement moléculaire

Supposons que le gaz n'est composé que d'hydrogène moléculaire H₂ et d'hélium 4 neutre He :

$$n(H_2) = n_H/2; n(H), n(H^+), n(He^+), n(He^{2+}), n(e^-) \ll n_H \quad (IV.18)$$

$$n_{He} = n(He) = x(He)n_H \quad (IV.19)$$

Par conséquent, selon les équations de la Sect. IV.2.2.1 :

$$\bar{m} = 1 + 4x(He), \quad n_{tot} = (0.5 + x(He))n_H \quad \implies \mu(\vec{r}, t) = \frac{1 + 4x(He)}{0.5 + x(He)} \quad (IV.20)$$
### 2.2.5. Gaz parfait complètement atomique et neutre

Supposons que le gaz n'est composé que d'hydrogène atomique neutre H et d'hélium 4 neutre He :

$$n(H) = n_H; n(H^+), n(He^+), n(He^{2+}), n(e^-) \ll n_H; n(H_2) \ll n_H \quad (IV.21)$$

$$n_{He} = n(He) = x(He)n_H \quad (IV.22)$$

Par conséquent, selon les équations de la Sect. IV.2.2.1 :

$$\bar{m} = 1 + 4x(He), \quad n_{tot} = (1 + x(He))n_H \quad \implies \mu(\vec{r}, t) = \frac{1 + 4x(He)}{1 + x(He)} \quad (IV.23)$$

**Remarque :** Même si en pratique, le code modifié peut maintenant traiter les différentes compositions de gaz précédentes, dans les Chapitres V et VI, nous supposerons que le gaz simulé reste complètement atomique et neutre avec $x(He) = 0.1$.

## 3. Enveloppe environnante modélisée

Dans cette section, je résume la liste des paramètres d'entrées lus par le code; puis, je décris les conditions initiales (à $t = 0$) en densité et en vitesse dans l'enveloppe entourant le jet; enfin, je décris comment les conditions aux bords du domaine de simulation sont traitées (à chaque $t \ge 0$).

### 3.1. Liste des paramètres libres

La Table 1 résume la liste des paramètres libres principaux pour les enveloppes adoptées dans les Chapitres V et VI.

### 3.2. Conditions initiales

#### 3.2.1. Champs de densité et de vitesse du cœur

Différentes distributions d'environnement sont écrites pour la densité et la pression (Chapitres V et VI). En exploitant la linéarité entre la densité $\rho$ et la pression $p$ (Eq. (IV.1)), un champ uniforme de densité $\rho_{a0}$ et de pression $p_{a0}$ peut être transformé en champs stratifiés en les multipliant par un même champ adimensionné, qui rend compte de la stratification de l'enveloppe :

Partout dans le domaine simulé, le champ de densité de l'enveloppe environnante est initialisé à $t = 0$ par un champ

$$\rho_a(R, Z, t = 0), \quad (IV.1)$$

préalablement implémenté et choisi dans le code par l'utilisateur. Par exemple, dans le Chapitre V, nous utiliserons 3 expressions de $\rho_a(R, Z, t = 0)$, pour étudier l'influence de la stratification de l'enveloppe environnante sur les flots poussés par un jet collimaté.

Avec une température de l'enveloppe supposée uniforme à $t = 0$ (choix adopté par simplicité) et donnée par l'utilisateur, le champ de pression $p(R, Z, t = 0)$ est ensuite déduit pour un gaz parfait à l'aide de la relation Eq. (IV.1), en fonction de la composition du gaz.
<table><caption>TABLE 1 – Liste des paramètres libres principaux et initiaux ($t = 0$) de l'enveloppe environnante utilisée dans les Chapitres V et VI</caption><thead><tr><th></th><th>Paramètre</th><th>Intervalle de valeurs<br>utilisé durant la thèse</th></tr></thead><tbody><tr><td></td><td>Stratification z<sub>c</sub> ou r<sub>0</sub></td><td>z<sub>c</sub> = 10<sup>16</sup>cm et r<sub>0</sub> = 2.5 × 10<sup>15</sup> cm</td></tr><tr><td>Densité ρ<sub>a0</sub> en (R, z) = (R<sub>j</sub>, 0) (<sup>a</sup>) ou en (r, θ) = (r<sub>0</sub>, π/2) (<sup>b</sup>)</td><td></td><td>1.6 × 10<sup>-20</sup> g cm<sup>-3</sup> – 1.6 × 10<sup>-18</sup> g cm<sup>-3</sup></td></tr><tr><td>Nombre volumique en nucléons d'hydrogène n<sub>H</sub> (<sup>c</sup>)</td><td></td><td>~7 × 10<sup>3</sup> cm<sup>-3</sup> – ~7 × 10<sup>5</sup> cm<sup>-3</sup></td></tr><tr><td>Température du gaz</td><td></td><td>30 K – 100 K</td></tr><tr><td>Vitesse moyenne du gaz</td><td></td><td>0 km/s</td></tr><tr><td>Composition du gaz</td><td></td><td>complètement atomique<br>m̄ = 1.4, μ = 1.27</td></tr></tbody></table>

a. si l'enveloppe est stratifiée suivant z

b. si l'enveloppe est stratifiée suivant (r, θ)

c. toujours en (R, z) = (R<sub>j</sub>, 0) <sup>a</sup> ou en (r, θ) = (r<sub>0</sub>, π/2) <sup>b</sup>.

## 3.3. Conditions aux bords de l'enveloppe

### 3.3.1. Cellules fantômes

La Fig. 23 schématise comment la totalité de la grille de calcul $\mathbb{D}$ est structurée, sur un exemple où des rangées $\mathbb{B}$ de 2 cellules dites fantômes (ghost cells) sont placées le long de chaque bord du domaine de calcul dit intégré $\mathbb{M} = \mathbb{D} \setminus \mathbb{B}$. Ce dernier désigne les cellules dont les valeurs sont avancées dans le temps par la résolution numérique des équations Eq. (III.1)-(III.3).

L'approche Volumes Finis de MPI-AMRVAC basée sur une description eulérienne (maillage fixe) justifie la nécessité de cellules fantômes : la condition appliquée au bord du domaine $\mathbb{M}$, pour chacun des champs résolus, est déterminée par la valeur que celle-ci prend dans les cellules fantômes. En outre, le code doit par exemple calculer des dérivées et gradients numériques, ou appliquer un limiteur de pente pour la reconstruction spatiale de ces dérivées à l'interface des cellules. Réaliser ces opérations dans $\mathbb{M}$ nécessite aussi que chaque cellule $(i, j)$ ait accès aux valeurs des champs de cellules voisines ...,$(i - 2, j)$, $(i - 1, j)$, $(i + 1, j)$, $(i + 2, j)$,... suivant $R$ et ...,$(i, j - 2)$, $(i, j - 1)$, $(i, j + 1)$, $(i, j + 2)$,... suivant $z$. Les cellules situées aux bords de $\mathbb{M}$ ont ainsi besoin des cellules fantômes pour réaliser ces opérations avec la même précision que pour les cellules plus internes.

### 3.3.2. Conditions aux bords en $R = 0, R = R_{max}, z = z_{max}, z = 0$

Pour les bords situés au sommet du domaine ($z = z_{max}^{(M)}$) et sur le bord droit ($R = R_{max}^{(M)}$), j'ai toujours utilisé des conditions aux bords ouvertes, à savoir telles que le gaz (venant de l'enveloppe ou du jet) puisse librement s'échapper du domaine de calcul.

L'axe $z$ du domaine de calcul ($R = 0$) étant l'axe de symétrie du système physique, aligné sur l'axe du jet, j'ai choisi d'y appliquer des conditions aux bords axisymétriques, c'est-à-dire :

— Des champs $\rho, \theta_{trjet}, \theta_{trenv}, v_z$, et $p$ symétriques par rapport à l'axe $z$. Si $(R, \theta) = (R_{gc}, -\theta_{gc})$
<figure><img src="image_32.png" /><figcaption>FIGURE 23 – Grille "totale" $\mathbb{D}$ de calcul (jaune+gris, indice (G)) avec des rangées de 2 cellules fantômes sur chaque bord : Grille de calcul intégrée/interne $M = \mathbb{D} \setminus \mathbb{B}$ (jaune, indice (M)) + rangées de cellules fantômes $\mathbb{B}$ (gris, indicé (B)).</figcaption></figure>
($\theta_{gc} > 0$) est la position d'une cellule fantôme au bord de $R = 0$ et $\theta$ l'angle polaire, alors :

$$W(R = R_{gc}, \theta = -\theta_{gc}, t) = W(R_{gc}, +\theta_{gc}, t), \quad \forall t \ge 0 \qquad (IV.2)$$

pour $W = \rho, \theta_{trjet}, \theta_{trenv}, v_z, p$.

— Des champs $v_R$ et $v_\phi$ antisymétriques par rapport à l'axe $z$ :

$$W(R_{gc}, -\theta_{gc}, t) = -W(R_{gc}, +\theta_{gc}, t), \quad \forall t \ge 0 \qquad (IV.3)$$

pour $W = v_R$ et $v_\phi$.

### Cas $z = 0$: des conditions "classiques" inadaptées

Pour nos simulations, je cherche des conditions aux bords en $z = 0$ qui limitent les effets de bords. D'après mes recherches, j'ai surtout réfléchi à gérer les problèmes d'aspiration par effet Venturi (par exemple à la frontière de la valve d'injection, cf. Sect. IV.4.2), à la formation de zones de dépression du gaz, ou encore à l'évidage progressif du réservoir de gaz dans les cellules fantômes.

J'ai ainsi d'abord testé la condition aux bords dite miroir en $z = 0$ (Fig. 24), c'est-à-dire où les champs scalaires et vectoriels sont symétriques par rapport au plan perpendiculaire à l'axe $z$ et passant par $z = 0$. En d'autres termes, j'ai pris des champs $\rho, v_R, v_\phi, p$ et des traceurs symétriques par rapport à $z = 0$; tandis que le champ $v_z$ est antisymétrique par rapport à $z = 0$ (il s'annule en $z = 0$).

Cependant, la condition aux bords miroir s'est vite avérée inadaptée dans nos simulations. En effet, dans la gamme de paramètres testés durant la thèse, nous avons vite pu voir qu'elle peut créer au cours du temps des zones de dépression et/ou d'aspiration (effet Venturi) du gaz le long de $z = 0$. En effet, les vitesses dans les cellules symétriques par rapport à $z = 0$ s'opposent au début de chaque pas de temps. Par exemple, le phénomène d'aspiration du gaz (Bernoulli grip) au voisinage de la partie du jet collée au bord peut éloigner le gaz de part et d'autre de $z = 0$ (Fig. 24). Une zone de dépression persistante apparaît derrière le gaz aspiré. Dans nos simulations avec cette condition aux bords, ces zones de dépressions apparaissent et grandissent au cours des 10 000 ans que nous cherchons à simuler.

Pour empêcher au réservoir de gaz des cellules fantômes d'être progressivement vidé, nous avons d'abord pensé à une autre condition aux bords "classique" consistant à bloquer l'entrée de matière dans le domaine intégré (voir Fig. 25). Cependant, ce choix s'est avéré d'autant plus inapproprié, en suivant le même raisonnement qu'avec la condition miroir : le gaz des cellules fantômes n'est pas aspiré (étant artificiellement "bloqué") lorsque du gaz monte depuis le bord inférieur, favorisant toujours la création de zones de dépression en $z = 0$. Dans les Chapitres V et VI, nous avons donc dû adopter d'autres conditions aux bords plus élaborées.

### 3.3.3. Conditions aux bords pour le cœur en $z = 0$

Devant les problématiques posées par les conditions précédentes, nous avons décidé d'utiliser deux conditions: les conditions figées et limit inflow. Je me suis occupé de les tester dans le cadre notre problème physique.

#### Conditions figées

Selon cette condition, à chaque pas de temps, les cellules fantômes du bord inférieur sont fixées à leurs valeurs initiales, c'est-à-dire telles que celles-ci sont évaluées à $t = 0$. Ainsi, si $(R, \theta) =$
<figure><img src="image_33.png" /><figcaption>FIGURE 24 – Schématisation de deux conditions aux bords (miroir et figées) sur un exemple en $z = 0$ et avec deux rangées de cellules fantômes par bord. Les couleurs indiquent la valeur en densité $\rho$ dans chaque cellule du domaine ($\rho$ augmente en allant successivement du cyan au bleu foncé, puis au violet et au rouge). Les flèches jaunes et cyan montrent les vecteurs des vitesses totales $\vec{v} = \vec{v}_R + \vec{v}_z$ dans ces mêmes cellules et dans le référentiel au repos de l'enveloppe. La première et la seconde ligne confrontent l'état des cellules dans les deux conditions aux bords, respectivement à l'instant $t = 0$ et $t > 0$.</figcaption></figure>
$(R_{gc}, \frac{\pi}{2} + \theta_{gc})$ est la position de chaque cellule fantôme du bord inférieur, alors pour tous $t \ge 0$ :

$$W(R_{gc}, \frac{\pi}{2} + \theta_{gc}, t) = W(R_{gc}, \frac{\pi}{2} + \theta_{gc}, t = 0) \qquad (IV.4)$$

où $W = \rho, \vec{v}, p, \theta_{trjet}, \theta_{trism}$.

La différence fondamentale avec la condition dite miroir standard (Fig. 24) est que cette dernière copie l'état changeant des champs à l'intérieur de $M$, à chaque instant $t \ge 0$.

De plus, cette dernière a un double avantage comparé à la condition miroir. D'abord, au cours des $t = 10 000$ ans, elle permet d'éviter un possible évidage du réservoir de gaz entraîné et aspiré dans le domaine intégré, puisque la valeur des champs résolus est remise à sa valeur initiale dans les cellules fantômes, au début de chaque pas de temps. Dans ce cas, le maintien d'une densité fixe à chaque pas de temps permet de remplir les zones de dépressions qui tendraient à apparaître. En parallèle, nous sélectionnerons pour les Chapitres V et VI des champs initiaux $\rho, p$ et $\vec{v}$ qui s'avèreront être symétriques par rapport à $z = 0$. Ainsi, utiliser la condition figée permet de simuler une condition de bord symétrique par rapport à $z = 0$, en ce qui concerne la partie de l'enveloppe qui n'est pas encore perturbée par le jet ou par les chocs d'étraves.

### Condition limit inflow

<figure><img src="image_34.png" /><figcaption>FIGURE 25 – Schématisation des conditions aux bords bloquantes en entrée et limit inflow (avec $\chi_B = 50\%$) sur un exemple en $z = 0$ et avec deux rangées de cellules fantômes par bord. Les couleurs indiquent la valeur en densité $\rho$ dans chaque cellule du domaine ($\rho$ augmente en allant successivement du cyan au bleu foncé, puis au violet et au rouge). Les flèches jaunes montrent les vecteurs des vitesses totales $\vec{v} = \vec{v}_R + \vec{v}_z$ dans ces mêmes cellules et dans le référentiel au repos de l'enveloppe.</figcaption></figure>

Cette autre condition généralise l'idée de la condition bloquante en entrée sous une forme plus
faible. La matière peut effectivement entrer à l'intérieur du domaine intégré $\mathbb{M}$, mais seulement partiellement. Plus précisément, la condition dite *limit inflow* impose que $\forall t, \Delta z \ge 0$:

$$v_z(R_{gc}, \frac{\pi}{2} + \theta_{gc}, t) = \begin{cases} \chi_B v(R_{gc}, \frac{\pi}{2} - \theta_{gc}, t, t) & , \text{ si } v_z(R_{gc}, \frac{\pi}{2} - \theta_{gc}, t) > 0, & (1) \\ v_z(R_{gc}, \frac{\pi}{2} - \theta_{gc}, t) & , \text{ si } v_z(R_{gc}, \frac{\pi}{2} - \theta_{gc}, t) < 0, & (2) \end{cases} \quad (IV.5)$$

tandis que les champs $W = \rho, v_R$ ou $p$ dans les cellules fantômes vérifient :

$$W(R, \frac{\pi}{2} + \theta_{gc}, t) = W(R_{gc}, \frac{\pi}{2} + \theta_{fc}, t), \quad \forall R \ge 0, \quad (IV.6)$$

où pour chaque rayon $R: \frac{\pi}{2} + \theta_{fc}$ est la position de la première cellule du domaine intégré accolée au bord $z = 0$ (Fig. 25); $\chi_B \in [0; 1]$ est un coefficient traduisant la perméabilité du bord $\mathbb{B}(z = 0)$. Le cas $\chi_B = 1$ correspond à des conditions ouvertes, tandis que $\chi_B = 0$ revient aux conditions bloquantes en entrée.

D'un côté, augmenter $\chi_B$ permet d'augmenter la quantité de matière venant des cellules fantômes et entrant dans le domaine intégré, pour remplir les zones de dépression ; d'un autre côté, diminuer $\chi_B$ permet de contrôler et diminuer la quantité de gaz s'échappant des cellules fantômes, les évidant. Nous pouvons ainsi jouer sur la valeur de $\chi_B$ pour équilibrer et limiter l'importance de ces deux effets problématiques. Mes tests ont montré que $\chi_B = 50\%$ constitue une valeur optimale au cours des 10 000 ans simulés. Par conséquent, la condition aux bords *limit inflow* est employée dans les simulations à temps longs des Chapitres V et VI.

# 4. Modélisation du jet pulsé

## 4.1. Liste des paramètres libres

Dans le cadre de nos simulations numériques de la propagation d'un jet pulsé à travers une enveloppe stratifiée, j'ai eu besoin de choisir un domaine de paramètres pour les paramètres libres des simulations. Le but est de permettre la comparaison de nos résultats avec les simulations de référence réalisées par Lee et al. (2001); puis une comparaison qualitative aux observations des flots HH46/47 et CARMA-7. La Table 2 résume la liste des paramètres libres principaux du jet que j'ai ainsi utilisés dans les Chapitres V et VI.

## 4.2. Valve d'injection: domaine et ouverture angulaire

Comme cela a été fait dans les simulations d'études précédentes sur les jets, toutes nos simulations s'intéressent **uniquement à la région de propagation du jet une fois déjà collimaté** (sans simuler les mécanismes de lancement sous-jacents). Le jet est ainsi lancé à partir d'un sous-domaine de cellules (appelé **valve d'injection du jet J**) localisé sur le bord inférieur de la grille de calcul $\mathbb{D}$, en $z = 0$ (Fig. 26 et Fig. 23, simulant une condition de débit entrant (*inflow*).

La valve J peut prendre plusieurs formes. Au cours de ma thèse, j'en ai implémenté deux, en fonction de la valeur de l'angle d'ouverture $\theta_j$ du jet que je souhaite simuler (Fig. 26): une **valve rectangulaire** ($\theta_j = 0$, Fig. 26i); ou une **valve trapézoïdale** ($\theta_j > 0$, Fig. 26ii).
<table><caption>TABLE 2 – Liste des paramètres principaux du jet pulsé utilisés dans les Chapitres V et VI</caption><thead><tr><th>Paramètre</th><th>Intervalle de valeurs<br>utilisée durant la thèse</th></tr></thead><tbody><tr><td>Rayon du jet<br>en (R, z) = (R<sub>j</sub>, 0)</td><td>50 au – 167 au</td></tr><tr><td>Densité de base<br>(en (R, z) = (R<sub>j</sub>, 0))</td><td>1.6 × 10<sup>-20</sup> g cm<sup>-3</sup> – 1.8 × 10<sup>-15</sup> g cm<sup>-3</sup></td></tr><tr><td>Nombre volumique en nucléons d'hydrogène<br>(aux mêmes points)</td><td>~7 × 10<sup>3</sup> cm<sup>-3</sup> – ~7 × 10<sup>8</sup> cm<sup>-3</sup></td></tr><tr><td>Température<br>Vitesse moyenne v<sub>j0</sub> du gaz<br>(en (R, z) = (0, z<sub>j</sub>))</td><td>100 K – 270 K<br>100 km/s – 120 km/s</td></tr><tr><td>Forme de la valve</td><td>rectangulaire ou trapézoïdale</td></tr><tr><td>Angle d'ouverture à la base</td><td>0° – 7°</td></tr><tr><td>Amplitude de variabilité Δv</td><td>60 km/s – 90 km/s</td></tr><tr><td>Période de variabilité P</td><td>115 ans – 300 ans</td></tr><tr><td>Profil de variabilité v<sub>j</sub>(t)/v<sub>j0</sub></td><td>sinus ou dent de scie</td></tr><tr><td>Niveau d'AMR minimum<br>demandé à l'intérieur du jet</td><td>1 – 5</td></tr><tr><td>Composition du gaz</td><td>complètement atomique<br>m&#773; = 1.4, μ = 1.27</td></tr></tbody></table>
### Valve rectangulaire

<figure><img src="image_35.png" /></figure>

<figure><img src="image_36.png" /></figure>

<figure><img src="image_37.png" /></figure>

### Valve trapézoïdale

<figure><img src="image_38.png" /></figure>

<figure><img src="image_39.png" /></figure>

<figure><img src="image_40.png" /></figure>

(i) Formes physiques possibles.

(ii) Implémentation sur la grille de (iii) Forme physique reconstruite en trois dimensions.

<caption>FIGURE 26 – Valve d'injection du jet</caption>

1. VALVE RECTANGULAIRE

Comme montré sur la Fig. 26i, le domaine d'une valve d'injection rectangulaire est donné par :

$$ \mathcal{J} = \{(R, z) \mid 0 \le R \le R_j \text{ et } z \le z_j\} \qquad (IV.1) $$

Reconstruit en trois dimensions, cette valve prend la forme d'un cylindre (Fig. 26iii).

2. VALVE TRAPÉZOÏDALE

Comme montré sur la Fig. 26ii, le domaine d'une valve d'injection trapézoïdale est donné par :

$$ \mathcal{J} = \{(R(z), z) \mid 0 \le R(z) \le R_j + z \tan(\theta_j), \forall z \le z_j\} \qquad (IV.2) $$
Reconstruit en trois dimensions, cette valve prend la forme d'un **cône tronqué** (Fig. 26iii).

!

Pour résoudre de manière précise l'injection du gaz à travers l'enveloppe en-
vironnante, nous avons imposé comme contrainte sur la résolution spatiale
de la simulation (et donc du niveau de raffinement des cellules de la valve)
que la largeur en z = 0 de la valve d'injection doit comporter au moins 10
cellules, comme cela est préconisé pour les valves d'injection de simulations
hydrodynamiques de jet dans des études précédentes (Wagner et al. 2012;
Gojon & Bogey 2017; Bosch-Ramon 2022). Cependant, la hauteur z_j n'est ja-
mais fixée à plus de ≃ 10 cellules, pour minimiser le décalage entre la hauteur
z_j d'injection et l'origine (R, z) = (0, 0) du repère utilisé pour initialiser les
champs initiaux de l'enveloppe.

4.3. Champs de densité et de vitesse initiaux dans la valve d'injection

Les champs de densité, de vitesses, de pression et de traceurs sont initialisés dans la valve à t = 0.

4.3.1. Pour une valve rectangulaire

Nous pouvons montrer que le flux de masse **hémisphérique** $\dot{\Phi}_m(z, t = 0) [^1]$ à travers la section transverse d'une valve rectangulaire $\mathcal{J}$ est conservée lorsque nous prenons des champs $(\rho, \vec{v})$ tels que (voir Fig. 27) :

— la vitesse du gaz dans toute la valve est uniforme et axiale :

$$
\vec{v}(R, z, t = 0) = v_{j0}\vec{u}_z, \quad \forall (R, z) \in \mathcal{J} \qquad (IV.3)
$$

— la densité dans la valve est uniforme :

$$
\rho(R, z, t = 0) = \rho_{j0}, \quad \forall (R, z) \in \mathcal{J} \qquad (IV.4)
$$

Par conséquent, en supposant que la température $T(R, z, t = 0) = T_j$ dans la valve est initialement uniforme, la pression du gaz dans $\mathcal{J}$ est alors elle aussi uniforme :

$$
p_{j0} = \frac{\rho_{j0} k_B T_j}{\mu m_H} \qquad (IV.5)
$$

Le jet lancé à partir d'un champ de vitesse axial donné par l'Eq. (IV.3) est appelé un **jet cylindrique** (Fig. 26).

Une valve rectangulaire a l'avantage d'être la forme la plus simple à implémenter pour lancer un jet cylindrique, ce qui justifie son utilisation récurrente dans plusieurs simulations d'études précédentes sur des jets protostellaires (e.g. Raga et al. 1990; Cabrit et al. 1997; Suttner et al. 1997; Lee et al. 2001; Moraghan et al. 2006; Tabone et al. 2018; Cerqueira et al. 2021).

<footnote n="1">1. Cette précision est nécessaire, car contrairement à notre étude, certains articles de simulations numériques (tels que Lee et al. (2001)) ou d'observations (e.g. Zhang et al. (2019)) donnent la somme des flux de masse des deux jets (ou des deux flots) de part et d'autre de l'objet central (par exemple en raison d'une résolution spatiale insuffisante pour discerner les flux hémisphériques).</footnote>
Valve rectangulaire

<figure><img src="image_41.png" /><figcaption>FIGURE 27 – Champs dans la valve rectangulaire.</figcaption></figure>

Pour la conservation du flux de masse **hémisphérique** $\dot{\Phi}_m(z, t = 0)$, nous pouvons en effet montrer qu’à chaque hauteur $z$ :

$$
\begin{align}
\dot{\Phi}_m(z, t=0) &= \int_{R'=0}^{R'=R_j} \int_{\varphi'=0}^{\varphi'=2\pi} [\rho(R(R', \varphi'), z(R', \varphi'), t=0) \nonumber \\
&\qquad \times \vec{v}(R(R', \varphi'), z(R', \varphi'), t=0) \cdot (R' d\varphi' dR' \vec{u}_z)] \nonumber \\
&= \rho_{j0} v_{j0} \pi R_j^2 = \dot{\Phi}_m(z=0, t=0) \tag{IV.6}
\end{align}
$$

4.3.2. Pour une valve trapézoïdale

Pour simuler le lancement d'un jet avec un angle d'ouverture non-nul, nous montrons facilement
que le flux de masse obtenu à partir des Eqs. (IV.3) et (IV.4) n'est pas conservé dans une valve
trapézoïdale :

$$
\begin{align}
\dot{\Phi}_m(z, t=0) &= \int_{R'=0}^{R'=R_j+z \tan(\theta_j)} \int_{\varphi'=0}^{\varphi'=2\pi} [\rho(R(R', \varphi'), z(R', \varphi'), t=0) \nonumber \\
&\quad \times \vec{v}(R(R', \varphi'), z(R', \varphi'), t=0) \cdot (R' d\varphi' dR' \vec{u}_z)] \nonumber \\
&= \rho_{j0} v_{j0} \pi (R_j + z \tan(\theta_j))^2 > \dot{\Phi}_m(z=0, t=0) \tag{IV.7}
\end{align}
$$

Pour conserver le flux de masse hémisphérique à travers une valve trapézoïdale, nous utilisons
des champs (ρ, v) tels que :

— la vitesse du gaz dans J soit orientée le long de lignes de courant divergentes à partir d’un
point virtuel (Fig. 28) :

$$
(R, z) = (0, -z_0) = (0, -\frac{R_j}{\tan(\theta_j)}) \qquad (\text{IV.8})
$$

et de sorte que la norme du vecteur vitesse reste uniforme ($v_{j0}$) dans la valve, ce qui donne
Valve trapézoïdale

<figure><img src="image_42.png" /><figcaption>FIGURE 28 – Champs dans la valve trapézoïdale.</figcaption></figure>

au final un champ de vecteurs vitesses initiaux :

$$ \vec{v}(R, z, t = 0) = v_{j0}\vec{u}_r \left( \Theta \stackrel{\text{def}}{=} \arctan\left(\frac{R}{z+z_0}\right) \right), \quad \forall (R, z) \in \mathbb{J} \qquad (\text{IV.9}) $$

$$ = v_{j0} \cos(\Theta) \vec{u}_z + v_{j0} \sin(\Theta) \vec{u}_R \qquad (\text{IV.10}) $$

$$ = v_{j0} \frac{z+z_0}{\sqrt{(z+z_0)^2 + R^2}} \vec{u}_z + v_{j0} \frac{R}{\sqrt{(z+z_0)^2 + R^2}} \vec{u}_R, \quad \forall (R, z) \in \mathbb{J} \qquad (\text{IV.11}) $$

,où $\vec{u}_r(\Theta = \arctan(\frac{R}{z+z_0}))$ est le vecteur unitaire radial et sphérique évalué à une distance $r = \sqrt{R^2 + (z+z_0)^2}$ de $(R, z) = (0, -z_0)$.

- la densité $\rho$ (et par conséquent la pression $p$) du gaz dans $\mathbb{J}$ prend des valeurs constantes sur chaque arc de cercle de rayon $r = z + z_0$ par rapport à $(0, -z_0)$ (ainsi que sur les calottes sphériques $C_S(r = z + z_0)$ de mêmes rayons $r$, en 3D).

- Nous pouvons montrer que pour conserver le flux de masse à travers chaque calotte sphérique $C_S(r = z + z_0)$ de la valve trapézoïdale, le champ de densité adapté s'écrit alors :

$$ \rho_j(R, z, t = 0) = \rho_{j0} \left( \frac{z_0^2 + R_j^2}{(z+z_0)^2 + R^2} \right), \quad \forall (R, z) \in \mathbb{J} \qquad (\text{IV.12}) $$

Le jet lancé à partir d'une valve $\mathbb{J}$ dont le champ de vitesse est donné par l'Eq. (IV.11) est appelé un **jet conique** (Fig. 26).
## 4.4. Variabilité du jet : un paramètre physique encore plus ou moins libre

Une caractéristique du modèle de jet est la variabilité temporelle de ce dernier (Raga et al. 1990; Cabrit et al. 1997; Lee et al. 2001, voir Chapitre II). Pour induire cette variabilité, le code met à jour à chaque pas de temps les champs de densité, de pression et de norme du champ de vitesses dans le domaine $J$ de la valve d'injection, ainsi que son prolongement dans les cellules fantômes :

$$v_j(t) = v_{j0} \times h(t) \qquad (IV.13)$$

$$\rho_j(t) = \frac{\rho_j(\vec{r}, t)}{h(t)} \qquad (IV.14)$$

où $h(t) = \frac{v_j(t)}{v_{j0}} = \frac{\rho_j(\vec{r}, t=0)}{\rho_j(\vec{r}, t)}$ pour tout $t > 0$. Le champ de densité initial $\rho_j(\vec{r}, t = 0)$ dans la valve est donné par l'Eq. (IV.4) ou (IV.12). Ici, la norme de la vitesse est proportionnelle à $h(t)$, tandis que la densité est inversement proportionnelle à $h(t)$, car nous avons choisi de garder un flux de masse constant dans le temps (comme pour Raga & Cabrit 1993; Lee et al. 2001) : $\dot{\Phi}_m(z, t > 0) = \int_{R'=0}^{R'=R_j} \int_{\varphi'=0}^{\varphi'=2\pi} [\rho(t > 0) \times \vec{v}(t > 0) \cdot (R'd\varphi'dR'\vec{u}_z)] = \dot{\Phi}_m(z, t = 0)$. À ce jour, nous disposons d'une certaine latitude dans les choix faits pour modéliser la variabilité temporelle, du fait des difficultés observationnelles et théoriques pour contraindre de manière précise l'historique en densité et en vitesses de ce phénomène au niveau des jets/chocs d'étraves. Par exemple, les durées estimées entre les épisodes d'éjection protostellaires vont de quelques dizaines d'années pour des flots jeunes comme CARMA 7 (Plunkett et al. 2015; Vorobyov et al. 2018) à plusieurs centaines d'années pour des flots a priori plus évolués comme HH 46/47 (Zhang et al. 2019) ou HH 111 (Erkal et al. 2021). Ainsi, en plus des limites de résolution angulaire des instruments, il est pour l'instant encore difficile de suivre le mouvement des nœuds brillants le long des flots, sur d'aussi longues durées.

Les simulations numériques sont une des solutions nécessaires à cette contrainte. Au cours des 30 dernières années, plusieurs profils de variabilité temporelle $h(t)$ ont été proposés pour $v_j(t)$ et leurs résultats comparés aux observations. Chacun a ses atouts et ses inconvénients, en fonction de la source observée qui est comparée aux simulations (e.g. Lee et al. 2001), ainsi que des caractéristiques du modèle simulé (e.g. Suttner et al. 1997). Durant ma thèse, je me suis intéressé à deux profils particuliers.

### 4.4.1. Profil de variabilité sinusoïdal

L'idée selon laquelle les objets de Herbig-Haro (HH) comme HH46/47 ou CARMA-7 peuvent génériquement être décrits par des nœuds internes successifs associés à une variabilité sinusoïdale de la vitesse du jet est bien établie depuis longtemps (voir par exemple Raga et al. 1990; Downes & Ray 1999; Reipurth & Bally 2001). Ainsi, dans un modèle de jet pulsé suivant un profil sinusoïdal, la vitesse $v_j(t)$ du jet dans la valve est variée suivant la loi :

$$v_j(t) = ||\vec{v}_j(t)|| = v_0 \left( 1 + \frac{\Delta V}{v_0} \sin \left( \frac{2\pi t}{P} \right) \right) \quad \left( h(t) = 1 + \frac{\Delta V}{v_0} \sin \left( \frac{2\pi t}{P} \right) \right), \qquad (IV.16)$$
où $v_0$ est la vitesse médiane du jet, $\Delta V$ est la variation maximale de la vitesse du jet $v_j(t)$, P est la période de variabilité associée;

Certaines études ont modélisé la variabilité observée dans les jets par une superposition de plusieurs termes sinusoïdaux avec des périodes $P_i$ plus ou moins différentes. Ces termes ont par exemple été abordés dans les modèles simulés par Downes & Ray (1999); Downes & Cabrit (2003) et Castellanos-Ramírez et al. (2018). Ce dernier a comparé les résultats obtenus avec les observations du jet HH 1, montrant qu'une variabilité multimodale de la vitesse du jet peut jouer un rôle important pour rendre compte de la morphologie et de la cinématique de certains flots observés.

Cependant, dans les Chapitres V et VI, nous faisons le choix simplificateur d'adopter un profil unimodal de variabilité, donné par l'Eq. (IV.16), mettant de côté le cas multimodal.

### 4.4.2. Profil de variabilité en dent de scie (sawtooth)

Dans le Chapitre VI, j'ai choisi d'employer le même profil en dent de scie que dans les simulations hydrodynamiques de jets moléculaires de Suttner et al. (1997):

$$v_j(t) = \|\vec{v}_j(t)\| = v_0 \left( 1 + \frac{\Delta V}{v_0} \left\{ 1 - 2 \cdot \mod \left( \frac{t}{P}, 1 \right) \right\} \right) \quad (IV.17)$$

$$\left( h(t) = 1 + \frac{\Delta V}{v_0} \sin \left( \frac{2\pi t}{P} \right) \right) \qquad (IV.18)$$

Suttner et al. (1997) emploient ce profil au lieu d'un profil sinusoïdal, car dans leurs simulations et avec leur résolution spatiale, ce dernier profil ne produisait pas de chocs internes au jet suffisamment forts pour dissocier la molécule de dihydrogène $H_2$, dissociation qu'ils cherchaient à reproduire car jouant un rôle dans leur traitement de la chimie pour cette molécule.

Les Chapitres V et VI ne simulent qu'un jet entièrement atomique, mais comme la variabilité temporelle reste un paramètre plus ou moins libre, j'ai voulu m'intéresser aux résultats de simulations utilisant un profil de variabilité en dent de scie. De plus, le Chapitre V étudiera l'influence du choix du profil de variabilité (sinusoïdal ou en dent de scie) sur la forme et la cinématique des flots prédits.

La Fig. 29 compare le profil de variabilité sinusoïdal avec celui en dent de scie.

## 5. Termes sources

Durant ma thèse, j'ai pris en compte et implémenté plusieurs termes sources physiques $S_{phys}$ dans les équations hydrodynamiques de la forme des Eqs. (III.4).

Les équations les plus générales que j'ai résolu dans les simulations des Chapitres V et VI sont:

$$\text{Conservation de la masse: } \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (IV.1)$$

$$\text{Conservation de l'impulsion: } \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot (\rho\mathbf{v}\mathbf{v}) = -\nabla[p] + \mathbf{F}_p, \qquad (IV.2)$$

$$\text{Conservation de l'énergie: } \frac{\partial e}{\partial t} + \nabla \cdot (e\mathbf{v}) = -\nabla \cdot (p\mathbf{v}) + E_p - n_H^2 \Lambda(T), \qquad (IV.3)$$
<figure><img src="image_43.png" /><figcaption>FIGURE 29 – Comparaison entre le profil de variabilité temporelle sinusoïdal et celui en dent de scie (tel qu’adopté par Suttner et al. 1997), tel qu’adopté dans les Chapitres V et VI.</figcaption></figure>

## 5.1. Force de contrebalancement $\vec{F}_p$ et terme source $E_p$

Appliqués à chaque pas de temps $t \ge 0$, les termes sources $\mathbf{F}_p(t) = \nabla [p(t=0)]$ et $E_p(t) = \nabla \cdot (p(t) \mathbf{v}(t))$ sont respectivement les forces et l'énergie qui appliquent d'une part une compensation par le gradient de la pression initiale $p(t=0)$ des forces de pression thermique dans l'enveloppe statique, et d'autre part la divergence du champ $p(t)\vec{v}(t)$ de l'enveloppe à chaque pas de temps :

$$ -\nabla [p(t=0)] + \mathbf{F}_p(t) \stackrel{\text{def}}{=} \mathbf{0}, \quad \forall t \ge 0 \qquad (\text{IV.4}) $$

$$ -\nabla \cdot (p(t)\mathbf{v}(t)) + E_p(t) \stackrel{\text{def}}{=} 0, \quad \forall t \ge 0 \qquad (\text{IV.5}) $$

Ces termes sources visent à maintenir en équilibre hydrostatique le gaz de l'enveloppe qui n'est pas encore perturbé par le jet et le flot qu'il creuse. Ils visent à annuler l'effet des forces de gradient de pression crées dans une enveloppe de gaz parfait stratifiée en densité.

### 5.1.1. Evaluation de $\vec{F}_p$ : méthode numérique

Une première manière pour évaluer la force $\vec{F}_p$ est d'approximer numériquement le gradient $\nabla [p(t=0)]$ par la méthode des différences finies centrées, spatialement précise à l'ordre deux (ci-après en coordonnées cartésiennes et en deux dimensions) :

$$ \nabla(p) = \frac{p(x_{i+1,j}, y_{i,j}) - p(x_{i-1,j}, y_{i,j})}{2\Delta x} \mathbf{e}_x + \frac{p(x_{i,j}, y_{i,j+1}) - p(x_{i,j}, y_{i,j-1})}{2\Delta y} \mathbf{e}_y, \qquad (\text{IV.6}) $$
où : $\Delta x = x_{i+1,j} - x_{i,j}$ et $\Delta y = y_{i,j+1} - y_{i,j}$.

Enfin, c'est cette méthode numérique approchée qui a été utilisée pour maintenir l'enveloppe en équilibre hydrostatique dans Rabenanahary et al. (2022).

### 5.1.2. Evaluation de $\vec{F}_p$ : méthode analytique

Après la rédaction de Rabenanahary et al. (2022) j'ai réussi à évaluer plus exactement la force $F_p$. Cette méthode (que j'emploie dans les Chapitres V et VI) évalue analytiquement l'expression de $\nabla [p(t = 0)]$, en utilisant les expressions des champs initiaux de densité et de pression de l'enveloppe. Les composantes de la force analytiquement évaluée sont calculées dans l'Appendice VII.A, dans le cas d'un cœur singulier aplati.

Indépendamment de la méthode d'évaluation de $\vec{F}_p$, la contribution $F_{p,k}$ de chacune de ses composantes à l'équation de conservation de l'impulsion est ajoutée suivant un même schéma d'intégration temporelle. Afin d'inclure ce terme source dans la résolution de l'équation de conservation de l'impulsion, le code offre une variété de choix de schémas. Ainsi, nous utilisons une méthode d'Euler explicite en temps, de type prédicteur-correcteur en deux étapes, et précis à l'ordre deux (voir justification à la Sect. IV.5.1.4) :

$$ (\rho v_k)^{n+1/2} = (\rho v_k)^n + \frac{\Delta t}{2} F_{p,k} (t, (\rho v_k)^n) \qquad (IV.7) $$

$$ (\rho v_k)^{n+1} = (\rho v_k)^n + \Delta t F_{p,k} \left( t + \frac{\Delta t}{2}, (\rho v_k)^{n+1/2} \right) \qquad (IV.8) $$

où $(\rho v_k)(t) = (\rho v_k)^n$ est la densité d'impulsion suivant la direction $k = x, y$ évaluée au pas de temps $t$, et $(\rho v_k)(t + \Delta t) = (\rho v_k)^{n+1}$ est évaluée au pas de temps $t + \Delta t$.

### 5.1.3. Sens physique de l'énergie de contrebalancement $E_p(t)$

Sens physique de $-\nabla \cdot (pv)$

En appliquant la règle du produit sur les opérateurs différentiels, $-\nabla \cdot (pv)$ s'écrit:

$$ -\nabla \cdot (pv) = \underbrace{-\nabla [p] \cdot \mathbf{v}}_{\equiv W_{-\vec{\nabla}_p}} \underbrace{-p\nabla \cdot (\mathbf{v})}_{\equiv E_{el}(t)} \qquad (IV.9) $$

À chaque instant $t \ge 0$, l'énergie $-\nabla \cdot (pv)(t)$ se décompose ainsi comme la somme de deux énergies $-\nabla \cdot (pv)(t) = W_{-\vec{\nabla}_p} + E_{el}(t)$, où :

- $W_{-\vec{\nabla}_p} = -\nabla [p(t)] \cdot \mathbf{v}$ est le travail des forces de gradients de pression $-\nabla [p(t)];$

- En utilisant la règle des produits sur l'opérateur divergence :
$$-p(t)\nabla \cdot (\vec{v}(t)) = -p(t) \left[ \frac{1}{\rho(t)} \rho(t) \nabla \cdot (\vec{v}(t)) \right] \qquad (IV.10)$$

$$= -p(t) \left[ \frac{1}{\rho(t)} \{ \nabla \cdot (\rho(t)\vec{v}(t)) - \vec{v}(t) \cdot \nabla [\rho](t) \} \right] \qquad (IV.11)$$

$$= -\frac{p(t)}{\rho(t)} \left[ -\frac{\partial \rho(t)}{\partial t} - (\vec{v}(t) \cdot \nabla) \rho(t) \right] \qquad (IV.12)$$

$$= \frac{p(t)}{\rho(t)} \left( \frac{\partial}{\partial t} + (\vec{v}(t) \cdot \nabla) \right) \rho(t) \qquad (IV.13)$$

$$= \frac{p(t)}{\rho(t)} \frac{D\rho(t)}{Dt} = \frac{c_s^2(t)}{\gamma} \frac{D\rho(t)}{Dt} \qquad (IV.14)$$

(IV.15)

où j'ai utilisé la conservation de la masse (IV.1) dans un champ de vitesse $\vec{v} = \vec{v}(t)$. Main-tenant, lorsque le gaz est localement comprimé (resp. dilaté) alors qu'il est advecté par un champ $\vec{v}(t)$, la dérivée particulière $\frac{D\rho}{Dt}$ devient positive (resp. négative). Par conséquent, $-p(t)\nabla \cdot (\mathbf{v}(t)) > 0$ (resp. $< 0$): le gaz gagne de l'énergie (resp. il en perd). Par conséquent, la contribution $-p(t)\nabla \cdot (\mathbf{v}(t))$ peut être interprétée comme la **variation d'énergie du gaz provoquée par sa compression ou sa dilatation locale.**

L'énergie $E_p$ vient donc compenser d'une part le travail des forces de gradient de pression par le travail de la force $\vec{F}_p$ et d'autre part la variation d'énergie $E_{el}$ due à la compressibilité du gaz.

<figure><img src="image_44.png" /></figure>

### Cas particulier

Dans la limite des écoulements incompressibles, nous avons $\nabla \cdot (\vec{v}) = 0$, donc $E_{el}(t) = 0$.

### Intérêt de compenser $-\nabla \cdot (pv)$

Dans les Chapitres V et VI ainsi que dans l'article Rabenanahary et al. (2022), le terme source $E_p(t)$ amortit l'éventuelle variation en énergie (et à postériori en vitesses) induite par le gradient de pression local et la mise en mouvement du gaz de l'enveloppe par ce gradient, à chaque pas de temps. En effet, nous avons pu voir que l'application de la force $\vec{F}_p$ n'empêche pas parfaitement au gaz de l'enveloppe de se mouvoir, surtout sur des échelles de temps de $10^4$ans. Il peut s'agir d'erreurs de précision dans l'évaluation de $\vec{F}_p$ ou à la méthode numérique employée. La variation en vitesse crée alors une variation d'énergie dans le terme $\nabla \cdot (p(t)\mathbf{v}(t)) \neq 0$. $E_p$ agit ainsi pour compenser cette variation et renforcer le maintien du gaz de l'enveloppe à l'équilibre hydrostatique. Nous pouvons noter que sans toujours être expliqué et précisé explicitement dans les publications scientifiques, ce terme est systématiquement inclus dans les études précédentes de lancement de jet à travers une enveloppe (e.g. Lee et al. 2001).
### 5.1.4. Évaluation de $E_p$

À chaque pas de temps, l'énergie $E_p$ est approximée à l'aide d'un opérateur divergence numérique pré-implémenté dans le code, évalué par une méthode de différences finies centrées, et spatialement précis à l'ordre deux.

$$E - p = \nabla \cdot (p\vec{v}) = \frac{p(x_{i+1,j}, y_{i,j})v_x(x_{i+1,j}, y_{i,j}) - p(x_{i-1,j}, y_{i,j})v_x(x_{i-1,j}, y_{i,j})}{2\Delta x} + \frac{p(x_{i,j}, y_{i,j+1})v_y(x_{i,j}, y_{i,j+1}) - p(x_{i,j}, y_{i,j-1})v_y(x_{i,j}, y_{i,j-1})}{2\Delta y}, \quad (IV.16)$$

où : $\Delta x = x_{i+1,j} - x_{i,j}$ et $\Delta y = y_{i,j+1} - y_{i,j}$.

Comme pour les composantes de $\vec{F}_p$, $E_p$ est pris en compte pour résoudre les équations de conservation de l'énergie en utilisant un schéma d'intégration temporelle de type prédicteur-correcteur en deux étapes :

$$e^{n+1/2} = e^n + \frac{\Delta t}{2} E_p(t, e^n) \quad (IV.17)$$

$$e^{n+1} = e^n + \Delta t E_p \left( t + \frac{\Delta t}{2}, e^{n+1/2} \right) \quad (IV.18)$$

où $e(t) = e^n$ est l'énergie totale évaluée au pas de temps $t$, et $e(t + \Delta t) = e^{n+1}$ est évaluée au pas de temps $t + \Delta t$.

## 5.2. Refroidissement $-n_H^2 \Lambda(T)$ par désexcitations collisionnelles

La version de base du code fournie au début de la thèse permet déjà de traiter le refroidissement du gaz par désexcitations collisionnelles atomiques (voir Sect. IV.2.2).

Pour évaluer la variation d'énergie due au taux de refroidissement $\Lambda(T)$ du gaz (erg s⁻¹ cm³), nous avons choisi la méthode d'intégration exacte déjà disponible dans le code et développée par Townsend (2009). Cette méthode de traitement du refroidissement a l'avantage d'être plus précise, numériquement plus stable et plus rapide (Townsend 2009) que les méthodes eulériennes explicites, implicites et semi-implicites aussi disponibles dans le code. Cette méthode résout la température $T^{n+1}$ au pas de temps $t + \Delta t$ en fonction du taux de refroidissement $\Lambda(T^n)$ au pas de temps $t$, à l'aide de "fonctions d'évolution temporelle" (TEF) propres à cette approche. Le taux $\Lambda(T^n)$ est quant à lui évalué à partir de différents choix de tables ayant leurs propres domaines de validité en température. Le terme source $-n_H^2 \Lambda(T)$ est ensuite ajouté à l'énergie totale $e(t)$ suivant un schéma d'intégration temporelle d'Euler explicite en temps (précis à l'ordre un) :

$$e^{n+1} = e^n - \Delta t n_H^2 \Lambda(T) \quad (IV.19)$$
# 6. Tests de configuration pour notre problème physique

## 6.1. Choix des schémas numériques et limiteurs de pente

Dans cette section, je teste les choix de schéma numérique et de limiteur de pente (Sect. III.2.2) sur l'angle d'ouverture $\theta$ de la base de la cavité poussée par un jet.

<figure><img src="image_45.png" /><figcaption>FIGURE 30 – Comparaison des cartes de densité à $t = 300$ ans, pour une simulation de jet identique lancée avec différentes combinaisons de schémas numériques de flux et de limiteurs de pente.</figcaption></figure>