La Fig. 30 nous montre une capture à $t = 300$ ans des cartes de densité d'une même simulation de jet cylindrique pulsé. Nous adoptons une grille uniforme de résolution $\Delta R = 16$ au. La formation de chocs hydrodynamiques étant le phénomène au cœur de notre problème physique, je me suis intéressé aux solveurs de Riemann (Sect. III.2.2). Plus particulièrement ceux dits approchés, car plus rapides et moins gourmands en temps de calcul que des solveurs de Riemann exacts ou linéarisés tels que le solveur de Roe (Roe 1981)). Ainsi, chaque capture de la Fig. 30 montre un schéma numérique différent parmi TVDLF, HLL, HLLC et HLLCD; en même temps qu'un limiteur de pente parmi Koren, MINMOD et PPM.

L'enveloppe environnante est uniforme et 10 fois moins dense ($\rho_{a0} = 1.6 \times 10^{-19}$ g cm$^{-3}$) que le gaz injecté en $(R,z) = (R_j, 0)$ ($\rho_{j0} = 1.6 \times 10^{-18}$ g cm$^{-3}$). Le jet est pulsé sur une période $P = 50$ ans, avec une vitesse moyenne de $v_0 = 100$ km/s, pour une demi-amplitude $\Delta v = 30$ km/s.

### 6.1.1. Critères de sélection

Conformément à nos objectifs et aux données observationnelles employées pour l'étude de la morphologie des flots protostellaires, la sélection du schéma numérique s'est uniquement basée sur la morphologie de la cavité poussée par le jet à la base du flot. Nous cherchons un schéma qui minimise les instabilités qui apparaîtraient dans les chocs ainsi que les oscillations numériques identifiables dans cette région. Nous cherchons aussi le schéma et limiteur reproduisant le mieux la forme parabolique (voire conique) des couches de flots observés en CO (e.g. HH46/47, voir Zhang et al. 2016, 2019), en minimisant l'ouverture de la couche à la base du flot.

### 6.1.2. Résultats

En tenant compte de nos critères ainsi que du fonctionnement des schémas numériques susnommés, les solveurs HLL et TVDLF sont ceux qui maximisent l'ouverture de la base des flots, tout en leur donnant une forme davantage cylindrique que parabolique. En effet, nous avons vu dans le Chapitre III.2.2 que ce sont ces deux schémas qui font les approximations les plus fortes des ondes solutions dans le problème de Riemman (Fig. 21). TVDLF sélectionne l'onde la plus rapide, justifiant que les chocs résolus décrivent les bases de couches de flots les plus ouvertes parmi les schémas testés. D'autre part, HLL ne prend en compte que l'onde de choc (la plus rapide) et de raréfaction (la plus lente) de l'éventail de Riemann. En contraste, les solveurs HLLC et HLLCD sont moins diffusifs et permettent de réduire l'ouverture des couches de flots tout en maintenant une forme de cavité plus proche de celle observée. De plus, ces deux solveurs sont les plus exacts des solveurs.

Au niveau des limiteurs, nous pouvons d'abord déjà remarquer que le limiteur de Koren produit systématiquement des oscillations numériques identifiables à la base de la couche principale, indépendamment du schéma numérique. En effet, ce limiteur est d'un ordre de précision spatiale plus élevé que MINMOD, introduisant donc moins de dissipation numérique. Il permet de capturer plus précisément les ondes de chocs, mais aussi toutes perturbations instables (e.g. instabilités de pression dites pair-impair) dans l'onde de choc. Contrairement à MINMOD, ces perturbations sont moins bien dissipées et peuvent alors aboutir à des oscillations numériques. Koren est ainsi davantage adapté aux problèmes d'écoulements continus/lisses que pour des chocs forts (Marcowith et al. 2020) comme ceux à la base de notre couche. La PPM introduit davantage de dissipation numérique et empêche la formation d'oscillations, mais la forme de couche simulée s'éloigne encore de ce qui est observé. Par conséquent, le limiteur de pente le plus adapté à notre problème est MINMOD. Comme le souligne Marcowith et al. (2020), c'est en effet le limiteur de pente le plus
stable en présence discontinuités physiques fortes. De plus, c'est le limiteur de pente le plus rapide des trois.

En conclusion, nous utiliserons HLLC, du moins sur de courts temps d'intégration (< 1000 ans, Chapitre V), et employerons la variante hybride HLLCD pour des temps d'intégration plus longs (> 1000 ans, Chapitre VI), afin de tirer tout de même avantage de la précision de HLLC (là où il est appliqué) et de l'économie en temps de calcul apporté par TVDLF sur ces longues simulations (temps de calculs mesurés à ~2 – 4 fois plus rapide qu'avec HLLC, durant les tests). Pour les limiteurs de pente, nous utiliserons exclusivement MINMOD, réduisant davantage le temps de calcul d'un facteur ~1.5 par rapport aux autres limiteurs.

### 6.1.3. Schéma d'intégration en temps

Dans toutes nos simulations, nous utilisons un schéma d'intégration temporel avec une précision d'ordre deux en temps. Parmi tous les choix possibles, nous avons choisi un schéma de type prédicteur-correcteur en deux étapes (two-step predictor-corrector scheme). Dans nos équations de conservation résolues, et qui s'écrivent sous la forme semi-discrète :

$$ \frac{du}{dt} = \frac{1}{\Delta x} (F_{i-1/2} - F_{i+1/2}) = g(t, u) \qquad (IV.1) $$

et où $u(t) = u^n$ est une variable d'état résolue évaluée au pas de temps $t$ (densité, vitesses, énergie, traceurs), le schéma prédicteur-correcteur en deux étapes évalue l'état $u(t + \Delta t) = u^{n+1}$ au pas de temps $t + \Delta t$ comme :

$$ u^{n+1/2} = u^n + \frac{\Delta t}{2} g(t, u^n) \qquad (IV.2) $$

$$ u^{n+1} = u^n + \Delta t g \left( t + \frac{\Delta t}{2}, u^{n+1/2} \right) \qquad (IV.3) $$

## 6.2. Résolution spatiale : influence du niveau de raffinement par AMR

La Sect. III.1.1 présente aussi la fonctionnalité d'AMR du code. Lorsque la taille d'une cellule est divisée par 2 à chaque niveau de raffinement supplémentaire, la mémoire utilisée pour stocker la grille raffinée et les champs de variables associés augmente, limitant dans la pratique l'usage de niveaux d'AMR trop élevés et trop coûteux en mémoire. En outre, sur les machines que j'ai employées, la taille des cellules raffinées typiquement au-delà d'un niveau d'AMR de l'ordre de 20 ($\le 1/2^{19} = 1/524288$ fois la taille d'une cellule de niveau 1 d'une taille de 14 au) peut être affectée par les erreurs sur les opérations de virgule flottante, et ensuite affecter les opérations numériques (calculs de flux numériques, etc...) faites sur ces cellules.

Dans cette section, je détermine la valeur optimale de niveau maximal de raffinement de notre domaine, à partir d'une résolution de 83 au (au niveau 1) adoptée dans les simulations avec AMR des Chapitres V et VI. Ainsi, la Fig. 31 compare des flots avec des densités $n_H = 10^5 \text{cm}^{-3}$ pour l'enveloppe et $n_H = 10^7 \text{cm}^{-3}$ pour le jet, montrées à $t = 204$ ans en fonction du niveau maximal de raffinement. Le jet est ici 100 fois plus dense que l'enveloppe. Ici aussi : $P = 50$ ans, $v_0 = 100$ km/s, et $\Delta v = 30$ km/s.
<figure><img src="image_46.png" /></figure>

FIGURE 31 – Comparaison des cartes de densité à $t = 204$ ans, pour une simulation de jet identique lancée avec différentes valeurs de niveau maximal de raffinement (entre 1 et 7), avec une résolution spatiale entre 84 au et 1 au (voir Table 3).

### 6.2.1. Critères de sélection

Pour choisir le niveau de raffinement maximal pour nos simulations et optimiser le rapport gain de précision/temps de calcul, nous imposons comme critère seuil que la largeur totale de flot $W_{800}$ à $z = 800$ au varie de moins de 10%, après avoir itéré d'un niveau supplémentaire de raffinement maximal.

Pour mesurer la position de la largeur du flot, j'effectue une coupe de densité à $z = 800$ au et en partant de l'extrémité $R = R_{max}$ de la grille. La position recherchée est celle où le profil de densité atteint son premier maximum, correspondant au choc en avant délimitant la couche principale (Chapitre V).

### 6.2.2. Résultats

De cette façon, nous sommes descendus jusqu'au niveau 7 d'AMR. Le critère seuil sur $W_{800}$ est atteint à partir du niveau 5, même si cette quantité continue de diminuer de 9% en passant d'un niveau maximal de 5 à un niveau 6, et de 7% en passant du niveau 6 à 7.

Cependant, en parallèle de ce gain de précision, le temps de calcul croît exponentiellement avec le niveau d'AMR maximal autorisé (chaque niveau supplémentaire de raffinement divisant une même cellule 2D en 4 plus petites cellules <fnref n="1" />). La Table 3 compile la précision maximale spatiale obtenue et

<footnote n="1">1. faisant à la fois augmenter rapidement le nombre de cellules dans le domaine (et donc le nombre d'opérations numériques par pas de temps), mais aussi diminuer le pas de temps $\Delta t$ imposé par la condition de convergence (pour résoudre les équations aux dérivées partielles utilisées) de Courant-Friedrichs-Lewy (CFL) sur le nombre de Courant</footnote>
<table><caption>TABLE 3 – Résolution et temps de calcul en fonction du niveau maximal de raffinement</caption><thead><tr><td>Niveau max. de raffinement</td><td>Résolution max. (au)</td><td>Temps de simulation</td></tr></thead><tbody><tr><td>1</td><td>83</td><td>~5 min</td></tr><tr><td>2</td><td>42</td><td>~7 min</td></tr><tr><td>3</td><td>21</td><td>~30 min</td></tr><tr><td>4</td><td>11</td><td>~4 h</td></tr><tr><td>5</td><td>5</td><td>~15 h</td></tr><tr><td>6</td><td>3</td><td>~1 jrs</td></tr><tr><td>7</td><td>1</td><td>~3 jrs</td></tr></tbody></table>

le temps de calcul à chaque niveau maximal. Le temps de calcul pour les simulations dépassant le niveau 6 atteint 24 h, pour un temps d'intégration n'allant que jusqu'à seulement 604 ans. De plus, chaque fichier sorti dépasse 10 gigaoctet, rendant leur post-traitement irréalisable avec les outils à ma disposition (ParaView, Python). Or, nous cherchons dans les simulations des Chapitres V et VI à atteindre un âge de 10 000 ans.

Par conséquent, augmentant le niveau de raffinement au-delà d'un niveau 5 est vite apparu comme trop coûteux en temps de calcul, comparé à la diminution relative de largeur de flot $W_{800}$ induite à la base (< 10%).

## 6.3. Paramètres de parallélisation avec MPI

Pour finir, nous testons ici l'efficacité apportée au temps d'exécution du code par la parallélisation sous la norme MPI (voir Keppens et al. 2012; Keppens et al. 2021 pour l'implémentation de celle-ci dans le code), dans le cadre spécifique de notre problème. Nous voulons trouver un nombre optimal de cœurs de calculs à utiliser sur nos machines (en l'occurrence, ici Tycho).

Nos simulations étant réalisées sur la grappe de calcul interne Tycho de l'Observatoire de Paris, nous avons décidé d'y employer la version Intel® de la librairie OpenMPI (²). Nous pouvons ainsi tirer parti des performances des compilateurs des processeurs Intel® de Tycho.

### 6.3.1. Test d'extensibilité forte et loi d'Amdhal

Dans le cadre d'un test dit d'extensibilité forte : nous testons le gain d'efficacité lorsque le nombre de cœurs de calcul sur lesquels la charge totale est répartie augmente, tout en gardant une taille de domaine de simulation constante (en dimensions et nombre de cellules). La charge de travail par cœur de calcul diminue alors. Le test d'extensibilité forte est employé pour de longues simulations comme les nôtres, afin de trouver la configuration (nombre de cœurs) optimale sur le matériel utilisé, pour un coût modéré en ressources. En effet, je n'aurai pas toujours accès durant les périodes de simulation de la thèse au nombre maximum de cœurs.

L'accélération informatique $S_p$ (en anglais speedup) mesure ce gain de performance en fonction du nombre de cœurs $N_p$:

$$CFL_{k,l} = \Delta t \sum_{D=z,R} \frac{v_{D,k,l}}{\Delta D}$$

<footnote n="2">2. voir https://docs.hpc.shef.ac.uk/en/latest/sharc/software/parallel/openmpi-intel.html</footnote>
$$S_p = t(1)/t(N_p) \qquad (IV.4)$$

où $t(1)$ est le temps de calcul d'un programme en séquentiel (sur un seul cœur), et $t(N_p)$ est le temps d'exécution du même programme sur $N_p$ cœurs. Je cherche ainsi à optimiser $S_p$ et $N_p$. Cette optimisation dépend aussi des performances des processeurs servant au calcul en parallèle, mais aussi de par la fraction du programme MPI-AMRVAC qui peut être passée sur les parties parallélisées. En mettant d'autres facteurs (comme les fluctuations de performances des processeurs), cette fraction ajuste théoriquement la courbe d'accélération $S_p = S_p(N_p)$ suivant la loi d'Amdahl (Amdahl 1967):

$$S_p = 1/(s + p/N_p) \qquad (IV.5)$$

où $s$ est la proportion de temps d'exécution passée sur les parties séquentielles; $p$ est la proportion du même temps passée sur les parties parallélisées.

Dans le cas idéal, l'accélération informatique d'un problème devrait évoluer linéairement, et accélérer un programme d'un facteur $N_p$ sur $N_p$ cœurs de calculs. Dans ce cas idéal, le code passe l'intégralité de son temps d'exécution sur des parties parallélisées. Cependant, $s > 0$ dans la pratique et il arrive un seuil de $N_p$ où cette loi de proportionnalité ne tient plus. Afin d'améliorer l'accélération $S_p = S_p(N_p)$ sans modifier le programme, il faudrait alors une machine de calcul plus puissante, ce qui n'est pas toujours faisable ni envisageable.

### 6.3.2. Méthode

Sur différents nombres de cœurs, nous lançons la même simulation de référence d'un jet conique pulsé (Sects. IV.4.2 et IV.4.3) lancé à travers un cœur singulier aplati (e.g. Chapitre V) jusqu'à $t = 1000$ ans, tout en gardant une taille de domaine constante. Pour simplifier le calcul, nous n'activons pas l'AMR au sein du code et décidons de lancer la simulation sur un domaine d'extension $(R, z) \in [0 \text{ cm}, 3 \times 10^{17} \text{ cm}] \times [0 \text{ cm}, 3 \times 10^{17} \text{ cm}]$, avec $n_R \times n_z = 240 \times 240$ cellules, comme dans le Chapitre V au niveau 1 d'AMR. La simulation est lancée sur le même nœud de calcul de Tycho (tycho70), pour minimiser tout biais lié à une différence de caractéristiques entre les nœuds de Tycho. Mis en service sur Tycho en Mars 2021, le nœud tycho70 dispose d'un processeur Intel Xeon Gold 5220R cadencé à 2.20GHz, pour un total de 48 cœurs de calcul, une mémoire vive de 192 Go, et un stockage de 1.5 To. Pour chaque $N_p$ donné, le temps de calcul est aussi moyenné sur 8 simulations identiques, pour tenir compte des fluctuations statistiques dans les opérations de virgule flottante, ainsi que dans les communications entre les cœurs du nœud.

### 6.3.3. Résultats

La Fig. 32 montre l'accélération informatique $S_p$ mesurée en fonction du nombre de cœurs de calculs $N_p$, décomposé selon trois temps d'exécution complémentaires, que le code peut mesurer et afficher à la fin de chaque simulation : le temps total ($t_{tot}$) de la simulation, le temps passé sur la partie lecture/écriture de fichiers ($t_{IO}$), et celui passé en dehors de ces lectures/écritures ($t_r$, temps dans la boucle d'intégration).

Premièrement, la parallélisation de notre problème dans MPI-AMRVAC permet d'accélérer le temps de calcul (total et hors lecture/écriture) en augmentant le nombre de cœurs. Par rapport à une exécution séquentielle, le gain atteint un facteur $\gtrsim 15-30$ avec 48 cœurs.
<figure><img src="image_47.png" /><figcaption>FIGURE 32 – Accélération informatique du temps de calcul pour notre programme avec MPI-AMRVAC, dans le cas d’un test d’extensibilité forte. Les courbes noire, bleue et verte sont associés respectivement au temps calcul et à l’accélération informatique totaux ($t_{tot}$) de la simulation, passés sur la partie lecture/écriture de fichiers ($t_{IO}$), et passés en dehors de ces lectures/écritures ($t_r$).</figcaption></figure>

Deuxièmement, je confirme que l'accélération informatique du programme s'ajuste sur la loi d'Amdahl Eq. (IV.5), de paramètres $s \simeq 5.1\%$ et $p = \simeq 94.9\%$ pour $t_{tot}$, ainsi que $s \simeq 1\%$ et $p = \simeq 99\%$ pour $t_r$. Ces valeurs donnent une idée des proportions séquentielles et parallélisées du programme avec MPI-AMRVAC. L'adéquation avec la loi théorique confirme que la perte de performance par rapport au cas idéal semble effectivement liée à une proportion non-nulle de temps passé sur les parties séquentielles du code.

Troisièmement, nous confirmons qu'entre ~32 et ~48 cœurs, l'accélération informatique gagnée commence à saturer vers $S_p = \frac{1}{s+p/48} \sim 14$ pour $t_{tot}$. Par conséquent, au-delà de ~24 cœurs, le gain de temps de calcul total lié à la parallélisation sous MPI devient négligeable dans le cas de notre problème. Lorsque les ressources de la gruppe de calcul Tycho utilisées par d'autres utilisateurs à l'Observatoire de Paris se font moins disponibles, j'ai ainsi pu économiser le nombre de cœurs demandés en descendant jusqu'à 24 cœurs.
<figure><img src="image_48.png" /></figure>
# CHAPITRE

## V

## FLOTS POUSSÉS PAR UN JET VARIABLE: INFLUENCE DE L'ENVIRONNEMENT ET DE L'ÂGE

### Sommaire

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^4$ 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
**Résumé** : Pour la première fois, à l'aide de simulations numériques, nous étudions l'influence que la stratification de l'enveloppe environnante peut avoir sur la cinématique et la morphologie des flots résultants de l'entraînement par un jet variable (Sect. V.1).

Nous comparerons aussi la forme des flots poussés par un jet pulsé cylindrique avec la forme des flots poussés par un vent grand angle, lancé avec le même flux de masse à travers un cœur singulier aplati identique (Sect. V.1.4).

Puis, nous étudierons l'influence des paramètres libres du problème sur la forme et la cinématique du flot protostellaire produit par un jet pulsé conique, à travers un cœur singulier aplati (Sect. V.2).

Enfin, nous présenterons les premières simulations de l'évolution du flot à temps long $\simeq 10^4$ ans, dans un cas de référence. Nous montrerons que cela résoud les trois critiques formulées dans le passé à l'encontre de ce type de modèles (Sect. V.3).

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

## 1.1. Paramètres du jet et du domaine

Dans les simulations des Sects. V.1 et V.1.4, nous adoptons la même géométrie et les mêmes conditions initiales dans la valve d'éjection que dans la simulation de jet pulsé de Lee et al. (2001), à savoir un jet pulsé cylindrique avec :

- un rayon $R_j \simeq 167$ au,

- une température $T_j = 270$ K,

- un profil sinusoïdal de variabilité temporelle :

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

où $v_0 = 120$ km/s est la vitesse médiane du jet, $\Delta V = 60$ km/s est la semi-amplitude de variation de la vitesse du jet au cours du temps, P est la période de variabilité du jet,

- un profil de variabilité de la densité $\rho_j(t)$ inverse de celle de $v_j(t)$ :

$$\rho_j(t) = \frac{\rho_{j0}v_j(t)}{v_0} = \frac{\rho_{j0}}{1 + \frac{\Delta V}{v_0} \sin\left(\frac{2\pi t}{P}\right)}, \qquad (V.2)$$

de façon à induire un flux de masse hémisphérique stationnaire

$$\dot{M} = \rho_0 v_0 \pi R_j^2 \qquad (V.3)$$

Nous adoptons aussi :

- la même taille de domaine $(R, z) = 4679$ au $\times 9358$ au et une résolution spatiale uniforme de 14 au,

- la même loi de refroidissement radiatif : celle de Dalgarno & McCray (1972) en dessous de $T_c = 10^4$ K et celle de MacDonald & Bailey (1981) au-dessus de $10^4$ K,

- une température identique pour le cœur environnant ($T_a = 100$ K).

Cette configuration est résumée dans la Table 4, avec les valeurs spécifiques de chaque modèle.
<table><caption>TABLE 4 – Paramètres des simulations des jets pulsés cylindriques de la Sect. V.1</caption><thead><tr><th colspan="2">Paramètres communs (de Lee et al. 2001)</th></tr></thead><tbody><tr><td>flux de masse hémisphérique du jet</td><td>M = 6.0 × 10<sup>-8</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>demi-ouverture angulaire du jet</td><td>θ<sub>j</sub> = 0°</td></tr><tr><td>rayon du jet</td><td>R<sub>j</sub> = 2.5 × 10<sup>15</sup> cm = 167 au</td></tr><tr><td>vitesse moyenne du jet</td><td>v<sub>0</sub> = 120 km/s</td></tr><tr><td>variation en vitesse du jet</td><td>v<sub>j</sub>(t) = v<sub>0</sub> + ΔV sin (2πt/P)</td></tr><tr><td>demi-amplitude</td><td>ΔV = 60 km/s</td></tr><tr><td>variation en densité du jet</td><td>ρ<sub>j</sub>(t) = ρ<sub>j0</sub> [v<sub>0</sub>/v<sub>j</sub>(t)] (flux de masse constant)</td></tr><tr><td>densité initiale du jet</td><td>ρ<sub>j0</sub> = 1.6 × 10<sup>-20</sup> g cm<sup>-3</sup></td></tr><tr><td>température du jet</td><td>T<sub>j</sub> = 270 K</td></tr><tr><td>temperature de l'enveloppe</td><td>T<sub>a</sub> = 30 K</td></tr><tr><td>fonction de refroidissement radiatif</td><td>Λ (270 K ≤ T < 10<sup>4</sup> K) à partir de Dalgarno &amp; McCray (1972)<br>Λ (T ≥ 10<sup>4</sup> K) à partir de MacDonald &amp; Bailey (1981)<br>(R, z) = (7.0 × 10<sup>16</sup> cm, 1.4 × 10<sup>17</sup> cm) = (4679 au, 9358 au)</td></tr><tr><td>domaine de simulation</td><td>n<sub>R</sub> × n<sub>z</sub> = 336 × 672</td></tr><tr><td>nombre de cellules</td><td>ΔR = Δz = 2.08 × 10<sup>14</sup> cm = 14 au</td></tr><tr><td>Résolution spatiale</td><td></td></tr></tbody></table><table><thead><tr><th colspan="4">Paramètres spécifiques de chaque modèle</th></tr><tr><th>Nom du modèle</th><th>PCJ-U<sup>a</sup><br>(Fig. 33a)</th><th>PCJ-Z<br>(Fig. 33b)</th><th>PCJ-R<br>(Fig. 33c)</th><th>PCJ-RW<sup>b</sup><br>(Fig. 38)</th></tr></thead><tbody><tr><td>Profile de densité de l'enveloppe</td><td>Uniforme</td><td>stratifié en z</td><td>Cœur singulier applati</td><td></td></tr><tr><td>ρ<sub>a</sub>(r)</td><td>ρ<sub>a0</sub></td><td>ρ<sub>a0</sub> (1 + z/z<sub>c</sub>)<sup>-2</sup><br>z<sub>c</sub> = 10<sup>16</sup> cm</td><td>ρ<sub>a0</sub> sin<sup>2</sup> θ r<sub>0</sub><sup>2</sup> (R<sup>2</sup> + z<sup>2</sup>)<sup>-1</sup><br>r<sub>0</sub> = 2.5 × 10<sup>15</sup> cm</td><td></td></tr><tr><td>Densité du cœur en z = 0, R = 2.5 × 10<sup>15</sup> cm</td><td>ρ<sub>a0</sub> = 1.6 × 10<sup>-20</sup> g cm<sup>-3</sup></td><td></td><td>1.6 × 10<sup>-18</sup> g cm<sup>-3</sup></td><td></td></tr><tr><td>contraste de densité cœur/jet en z = 0</td><td>η = ρ<sub>j0</sub>/ρ<sub>a0</sub> = 1</td><td></td><td>0.01</td><td></td></tr><tr><td>période de variabilité du jet</td><td>P = 310 ans</td><td></td><td>115 ans</td><td></td></tr><tr><td>Age de simulation présentée</td><td>t = 610 ans</td><td></td><td>296 ans</td><td></td></tr></tbody></table>

<footnote><sup>a</sup> paramètres identiques à ceux de la simulation de jet pulsé de Lee et al. (2001) (cf. leur Fig. 7)</footnote>

<footnote><sup>b</sup> paramètres identiques à ceux de la simulation de vent grand angle de Lee et al. (2001) (cf. leur Fig. 12), mais où ce dernier est remplacé par un jet cylindrique ayant le même flux de masse, rayon d'injection et variabilité en vitesse.</footnote>
## 1.2. Choix des modèles de stratification du cœur

Dans le Chapitre IV, nous avons vu qu'il est possible d'initialiser la densité de l'enveloppe suivant différents choix possibles de champs préalablement implémentés. Dans nos simulations, nous avons pris trois champs où nous avons choisi de fixer la même valeur de densité $\rho_{a0} = 1.6 \times 10^{-20} \text{ g cm}^{-3}$ en $(R, z)=(R_j, 0)$.

En premier lieu, dans le but d'une comparaison aux résultats du modèle de jet de Lee et al. (2001), nous simulons d'abord le cas classique d'un jet cylindrique lancé à travers une enveloppe environnante statique uniforme en densité et en pression, avec $\bar{v}(t=0) = \bar{0}$.

Le profil de densité de l'enveloppe est ainsi exprimée comme :

$$ \rho_a(R, z) = \rho_{a0}, \quad \text{avec} \quad \rho_{a0} = 1.6 \times 10^{-20} \text{ g cm}^{-3}, \qquad (V.4) $$

Par la suite, nous désignons ce modèle simulé par la dénomination PCJ-U (U=uniforme), qui est en fait la même simulation que celle de la Fig. 7 de Lee et al. (2001).

En second lieu, nous nous intéressons au cas d'un environnement stratifié suivant la direction z, tel que prescrit par Cabrit et al. (1997):

$$ \rho_a(R, z) = \rho_{a0} \left(1 + \frac{z}{z_c}\right)^{-2} \quad \text{avec} \quad z_c = 700 \text{ au, } \rho_{a0} = 1.6 \times 10^{-20} \text{ g cm}^{-3}, \qquad (V.5) $$

où la valeur de l'échelle caractéristique d'aplatissement $z_c$ est justifiée par des observations de cœurs protostellaires montrant une stratification abrupte de densité au moins à partir de 1000 au et jusqu'à 20 000 au (see e.g., Motte & André 2001; van Kempen et al. 2009).

Par la suite, nous désignons ce modèle simulé par la dénomination PCJ-Z.

En dernier lieu, nous adoptons le même champ stratifié dit de "coeur singulier aplati" (flattened singular core en anglais) que dans les simulations de vent grand angle de Lee et al. (2001), à savoir :

$$ \rho_a(r, \theta) = \rho_{a0} \sin^2 \theta \left(\frac{r_0}{r}\right)^2 \quad \text{avec} \quad r_0 = 2.5 \times 10^{15} \text{cm}, \rho_{a0} = 1.6 \times 10^{-20} \text{ g cm}^{-3}, \qquad (V.6) $$

où $r = (R^2 + z^2)^{1/2}$ désigne le rayon sphérique et $\theta$ est l'angle polaire par rapport à l'axe du jet.

Par la suite, nous désignons ce modèle par la dénomination PCJ-R.

La décroissance en $1/r^2$ de la loi de densité Eq. (V.6) est appropriée pour un cœur isotherme singulier auto-similaire. La dépendance en $\sin^2 \theta$ de la loi de densité est quant à elle une approximation de la solution d'équilibre magnétostatique de Li & Shu (1996) pour $n = 2$ (où ce paramètre mesure le degré de support magnétique and d'aplatissement équatorial du cœur). Comme nous l'avons abordé dans le Chapitre II, ce champ est favorisé pour expliquer la distribution masse-vitesse et la forme des flots observés. Cependant, à notre connaissance, le modèle d'un cœur singulier aplati n'a jusqu'ici jamais été utilisé dans les simulations de flots induits par des jets.

Comme la densité décroît de plus en plus vite avec la distance à l'origine en passant de PCJ-U à PCJ-Z, puis à PCJ-R, nous dirons par la suite que "le degré de stratification de l'enveloppe augmente".
<figure><img src="image_49.png" /></figure>

FIGURE 33 – Cartes de densité $n_H = n(H + H^+) + 2n(H_2)$, à $t = 296$ ans et $t = 610$ ans, de nos trois simulations avec le même jet cylindrique pulsé que dans Lee et al. (2001), pour trois stratifications différentes de la densité de l'enveloppe environnante : à gauche (a,d) : enveloppe uniforme $\rho(R, z) = \rho_{a0}$ (seul cas présent dans Lee et al. 2001); au milieu (b,e) : enveloppe stratifiée suivant $\rho(R, z) = \rho_{a0}/(1 + z/z_c)^2$ avec $z_c = 10^{16}$ cm; à droite (c,f) : cœur singulier aplati $\rho(r, \theta) = \rho_{a0} \sin \theta (r_0/r)^2$, où r est le rayon sphérique et $\rho_0 = 2.5 \times 10^{15}$ cm. Tous les profils ont la même densité en $z = 0$, $R = r_0 : \rho_{a0} = 1.6 \times 10^{-20} \text{ g cm}^{-3}$, valeur qui coïncide avec celle du jet à $t = 0$. Nous pouvons remarquer que le choc d'étrave le plus externe ainsi que les chocs d'étraves imbriqués deviennent d'autant plus large que le jet se propage à travers une enveloppe de plus en plus stratifiée. L'ensemble des paramètres de ces 3 simulations sont résumés dans la Table 4.

## 1.3. Résultats

### 1.3.1. Coupes de densité à $t = 300 - 600$ ans

La figure 33 présente les cartes de densité issues de nos simulations pour les trois différentes enveloppes que nous avons pris (Eqs. (V.4) à V.6). La comparaison est ici faite aux mêmes âges, à temps courts : $t = 296$ ans et 610 ans.
Les Fig. 33a et Fig. 33d nous montrent le cas du cœur environnant uniforme (PCJ-U,Eq. (V.4)).

La cavité creusée par le choc d'étrave à l'avant du jet s'élargit au fur et à mesure du temps et que le jet progresse dans le cœur environnant.

Nous pouvons vérifier que le jet lancé dans cet environnement y avance à la même vitesse que dans la simulation du même jet en Fig. 7 de Lee et al. (2001). À $t = 610$ ans, le choc d'étrave principal a la même longueur $z \simeq 1.3 \times 10^7$ cm, ainsi que les mêmes (demi-)largeurs maximales $z \simeq 1 \times 10^{16}$ cm. Ce résultat est cohérent avec le fait que le modèle uniforme PCJ-U ait exactement la même configuration que dans le modèle de jet pulsé présenté sur la Fig. 7 de Lee et al. (2001).

Comme déjà noté par Lee et al. (2001) et Ostriker et al. (2001), nous confirmons qu'un jet cylindrique lancé à travers un milieu uniforme génère une cavité étroite et de forme grossièrement cylindrique, contrairement aux observations de flots. Ce résultat constitue la principale objection d'Ostriker et al. (2001) pour rendre compte des formes de cavités observées à l'aide du modèle d'entraînement par un jet.

Les Fig. 33b et Fig. 33e. nous montrent le cas de PCJ-Z, où l'enveloppe est stratifiée selon $z$ (Eqs. (V.5)). En comparaison avec le cas uniforme, la pression dynamique exercée sur le choc d'étrave à l'avant du jet diminue et lui permet de se propager plus rapidement, élargissant plus rapidement la cavité. En effet, à $t = 610$ ans, la largeur totale du choc d'étrave à $z = 800$ au augmente, passant de $W_{800} = 0.95 \times 10^{16}$ cm pour l'enveloppe uniforme (PCJ-U) à $1.16 \times 10^{16}$ cm pour l'enveloppe stratifiée en $1/z^2$ (PCJ-Z). Au sommet du domaine ($z = 9358$ au), cette augmentation atteint même un facteur 4 en passant de PCJ-U à PCJ-Z. Ainsi, une stratification en densité abrupte en $1/z^2$ permet d'observer un élargissement général de la largeur du choc d'étrave à toutes les altitudes jusqu'en haut du domaine. Ce résultat diffère de la conclusion de la Section 3.6 de Lee et al. (2001), où la stratification en densité était décrite comme n'ayant pas d'effet notable. Cependant, ils avaient considéré un $z_c = 8350$ au douze fois plus long que dans notre modèle PCJ-Z, ne faisant chuter la densité d'enveloppe $\rho = \rho_{j0}/(1 + (z/z_c)^2)$ que d'un facteur deux sur l'ensemble du domaine de calcul. En outre, une échelle de 800 au est en meilleur accord avec les profils de densité déterminés dans les enveloppes protostellaires qu'avec la valeur de 8350 au alors utilisée par Lee et al. (2001) (e.g. van Kempen et al. 2009).

Alors que le choc d'étrave à l'avant du jet s'élargit plus rapidement, la matière qui se retrouve à l'arrière du choc, dans la cavité, est répartie dans un plus grand volume. En conséquence, la densité du gaz poussé par les chocs d'étraves imbriqués (et la pression dynamique exercée sur ces derniers) est plus faible (au même âge) que dans le cas d'une enveloppe uniforme : les chocs d'étraves imbriqués s'élargissent alors plus rapidement en passant de PCJ-U à PCJ-Z.

Un point important est cependant à remarquer. Même si les couches imbriquées deviennent plus larges suivant $R$, les zones d'interactions internes au jet (internal working surface en anglais) se déplacent le long de ce jet à la même vitesse, indépendamment de la stratification en densité du cœur environnant, et indépendamment de l'âge. Comme le montrent les relations de saut de Rankine-Hugoniot ainsi que les simulations de Raga et al. (1990), cela est dû au fait que la vitesse de propagation d'une zone d'interaction ne dépend que des conditions en vitesse et en densité du jet en amont et en aval de la zone d'interaction.

Or ces conditions sont uniquement et entièrement déterminées par les propriétés de variabilité du jet (i.e. $v_j(t)$ et $\rho_j(t)$), qui restent ici identiques quel que soit nos modèles de stratification en densité du cœur environnant.

Dans le cas d'un cœur singulier aplati (PCJ-R, Fig. 33c et Fig. 33f), la densité chute encore plus rapidement avec la distance à l'origine que pour une loi de décroissance en $1/z^2$, faisant
davantage chuter la pression dynamique exercée sur tous les chocs d'étrave, dont l'expansion se retrouve d'autant augmentée. Par exemple, à $t = 610$ ans, la largeur totale de choc d'étrave externe à $z = 800$ au augmente, allant de $W_{800} = 1.16 \times 10^{16}$ cm pour un cœur stratifié en $1/z^2$ (PCJ-Z) à $1.75 \times 10^{16}$ cm pour le cœur singulier aplati (PCJ-R). Au sommet du domaine de simulation, cette augmentation atteint même un facteur 2 entre PCJ-Z et PCJ-R, et un facteur 8 entre PCJ-U et PCJ-R.

### 1.3.2. Champs de vitesse

Les figures Fig. 34 et 35 montrent respectivement les cartes de vitesses radiale ($v_R$) et axiale ($v_z$) de nos trois simulations précédentes, aux mêmes âges : $t = 296$ ans et 610 ans.

Dans chacune des trois simulations, la vitesse radiale du gaz augmente à la sortie des zones d'interaction internes au jet pour atteindre son maximum au niveau des ailes des chocs d'étraves. La vitesse axiale, qui atteint un maximum de ~180 km/s le long de l'axe du jet, diminue en s'éloignant du faisceau du jet, avant de chuter rapidement de 50 à 0 km/s à partir d'une certaine distance radiale $R$.

D'un côté, dans le cas d'une enveloppe uniforme (panneaux a et d), les vitesses radiales atteintes dans les cavités creusées avec un jet ne dépassent pas ~10 km/s. Dans le cas stratifié suivant $z$ (b et e), la limite maximale est de $v_R$~30 km/s, tandis qu'elle atteint $v_R$~50 km/s pour le cœur aplati singulier (c et f). La Fig. 34 explique ainsi d'un point de vue cinématique le fait que les cavités creusées par chaque d'étrave s'ouvrent d'autant plus vite que le degré de stratification de l'enveloppe augmente. À $(R, z)$ donné, $v_z$ et $v_R$ augmentent en passant du cas uniforme au cas stratifié suivant $z$, suivi du cœur aplati singulier.

D'un autre côté, les cartes de vitesses axiales Fig. 35 nous montrent en dessous de ~50 km/s l'existence d'une structure de fort gradient de vitesse similaire à celle d'un cisaillement, et qui se développe le long des bords du choc d'étrave, dans les trois stratifications d'enveloppe. L'épaisseur de cette couche augmente avec le degré de stratification. Au sommet du domaine, elle est deux fois plus épaisse dans PCJ-Z que dans le cas uniforme PCJ-U, et jusqu'à 4 fois plus épaisse dans PCJ-R que dans PCJ-Z.

### 1.3.3. Diagrammes position-vitesse

La Fig. 36 présente les diagrammes position-vitesse (PV) longitudinaux dans PCJ-U, PCJ-Z et PCJ-R, calculés sur une épaisseur de 300 AU le long de l'axe du flot et avec une inclinaison de $i = 60°$ du jet par rapport à la ligne de visée (voir définitions et calculs du PV dans la Sect. VI.1). La méthode, ainsi que l'outil Python que j'ai programmé pour établir les diagrammes position-vitesse et les distributions masse-vitesse de mes simulations, sont décrites au début du Chapitre VI. Ici, nous avons représenté seulement la fraction de matière originaire de l'enveloppe environnante, $f_a = \theta_{tra}/(\theta_{tra} + \theta_{trj})$, à l'aide des traceurs $\theta_{tra}$, $\theta_{trj}$ (voir Sect. III.1.2.1). Ainsi, nous avons soustrait la contribution de la matière venant du jet, ne laissant ainsi idéalement à vitesses non-nulles que :

- la contribution de la matière de l'enveloppe balayée et comprimée au niveau des chocs d'étrave;

- la matière venant de l'enveloppe qui a re-rempli partiellement les cavités successives où elle s'est mélangée avec le matériau venant du jet (Raga & Cabrit 1993).

Tous les diagrammes PV de la Fig. 36 mettent en lumière les signatures caractéristiques de la cinématique des chocs d'étrave induits par des jets, à savoir : des structures convexes en forme
<figure><img src="image_50.png" /><figcaption>FIGURE 34 – Cartes de vitesse radiale $v_R$ des trois mêmes simulations que dans la Fig. 33. Nous pouvons remarquer que la vitesse atteinte dans les chocs d'étraves induits par le jet augmente d'autant plus que le jet se propage à travers une enveloppe de plus en plus stratifiée.</figcaption></figure>
<figure><img src="image_51.png" /><figcaption>FIGURE 35 – Cartes de vitesse axiale $v_Z$ des trois mêmes simulations que dans la Fig. 33. Ici, j’ai restreint l’affichage de la carte de couleurs en $v_z$ à 0–50 km/s, afin de garder le même intervalle de vitesse qu’atteint dans la Fig. 34, et qu’il s’agit de l’intervalle de vitesses typique observé dans les observations ALMA. Nous pouvons remarquer la décroissance de la vitesse axiale avec la distance au jet, formant systématiquement une structure de gradient de vitesse axiale formée au bord du choc d’étrave externe, et de plus en plus épaisse lorsque le degré de stratification augmente.</figcaption></figure>
<figure><img src="image_52.png" /></figure>

FIGURE 36 – Diagramme PV longitudinal le long de l'axe du flot et sommé sur une épaisseur $\Delta R = 300$ au, pour les modèles PCJ-U (jet cylindrique pulsé lancé à travers un cœur uniforme, gauche), PCJ-Z (cœur stratifié suivant $z$, au centre) et PCJ-R (cœur singulier aplati, à droite). On a supposé un angle d'inclinaison "médian" $i = 60°$ par rapport à la ligne de visée, à deux âges $t = 296$ et $t = 610$ ans. En suivant les traceurs et la fraction $f_a$ originaire de l'enveloppe, la masse représentée ne tient compte que de la fraction de gaz de l'enveloppe entrainée par les chocs d'étraves.

d'éperon (spur-like en anglais), telles que déjà prédites par les modèles analytiques (Masson & Chernin 1993; Cabrit et al. 1997; Ostriker et al. 2001) et vues dans les simulations avec enveloppe uniforme de Lee et al. (2001). Ces signatures (qui à $i = 60°$ sont présentes à la fois dans des vitesses observées négatives et positives) ne disparaissent donc pas aux âges jeunes $\le 600$ ans, même lorsque le cœur environnant est abruptement stratifié suivant les profils de densité que nous avons sélectionnés. Cependant, nous voyons que le degré de stratification du cœur environnant a un impact sur leur extension en vitesse. Les structures en forme d'éperon s'étendent ainsi à $t = 610$ ans entre -5 et 80 km/s pour PCJ-U, entre -10 km/s et $\sim 100$ km/s pour PCJ-Z, et entre -20 et 110 km/s pour PCJ-R. Plus le degré de stratification augmente, plus la pression du milieu ambiant
<figure><img src="image_53.png" /></figure>

FIGURE 37 – Distributions masse-vitesse calculées pour les simulations de jet pulsé cylindrique de la Fig. 33 lancé à travers un cœur environnant uniforme (bleu, PCJ-U), stratifié suivant *z* (orange, PCJ-Z), singulier aplati (vert, PCJ-R). Le jet est ici incliné de *i* = 60° par rapport à la ligne de visée. Les diagrammes sont évalués à *t* = 296 ans (haut) et *t* = 610 ans (bas); et les boîtes (*bins* en anglais) pour l'échantillonnage des vitesses observées *v*<sub>obs</sub> sont larges de 1 km/s en étant centrées sur ..., -1 km/s, 0 km/s, 1 km/s,... Nous isolons ici uniquement la contribution de la matière ambiante.

freinant l'expansion des chocs d'étrave est faible, et plus la matière emportée peut atteindre des vitesses maximales élevées.

À des niveaux plus bas, nous voyons dans les PV des signatures $v_{obs} \propto r$ (similaires à des *Hubble wedges* triangulaires) pouvant atteindre de hautes vitesses projetées $\ge 80$ km/s. Elles tracent le bord du jet où le gaz est mélangé avec une fraction de matière venant de l'enveloppe.

### 1.3.4. Distributions masse-vitesse

La Fig. 37 nous présente les diagrammes masse-vitesse (MV) pour les simulations PCJ-U, PCJ-Z et PCJ-R, avec un jet incliné d'un angle $i = 60^\circ$ par rapport à la ligne de visée, aux mêmes deux
âges que jusqu'ici. Nous isolons ici aussi la contribution de la matière originaire de l'enveloppe environnante (en supposant que c'est la seule qui serait observable en CO, si le jet est majoritairement atomique). Enfin, la figure montre séparément les composantes décalées vers le rouge et vers le bleu de la distribution. Les composantes en masse situées dans l'intervalle $v_{obs} \in [-0.5 \text{ km/s}; 0.5 \text{ km/s}]$, et majoritairement associées à l'enveloppe environnante encore statique, sont exclues.

Nous pouvons voir que comme démontré par Ostriker et al. (2001), comme dans les résultats du modèle analytique de Zhang & Zheng (1997), et comme dans les simulations de jets de Downes & Cabrit (2003); Moraghan et al. (2006), la distribution obtenue pour le lobe décalé vers le bleu se scinde entre deux parties à faibles ($<20 \text{ km/s}$) et hautes ($>20 \text{ km/s}$) vitesses, où la partie à faibles vitesses semble suivre une loi de puissance $m \propto v^{-\gamma}$ d'indice $\gamma$.

Nous pouvons observer qu'à ces âges courts, l'indice $\gamma$ de cette loi de puissance augmente avec le degré de stratification de l'enveloppe, en passant de PCJ-U ($\gamma(t = 296 \text{ ans}) = 1.05$ et $\gamma(t = 610 \text{ ans}) = 1.67$) à PCJ-Z ($\gamma(t = 296 \text{ ans}) = 1.21$ et $\gamma(t = 610 \text{ ans}) = 1.66$), puis à PCJ-R ($\gamma(t = 296 \text{ ans}) = 1.81$ et $\gamma(t = 610 \text{ ans}) = 1.99$). À titre de comparaison, Moraghan et al. (2006) ont réalisé une simulation de jet stationnaire et atomique, lancé à $v_j = 100 \text{ km/s}$ à travers une enveloppe atomique de densité $\rho(R, z) = \rho_a \propto [1 + z/z_c]^{-\kappa}$ ($\kappa > 0$), mais qui ne chute que d'un facteur 2 entre $z = 0$ et 16 000 au. Ils trouvent à $t = 1096 \text{ ans}$ des valeurs de $\gamma = 0.7 - 2.7$ (selon l'intervalle de $v_{obs}$ considéré) globalement plus faibles que pour PCJ-R à $t \le 610 \text{ ans}$.

En outre, dans nos simulations, la loi de puissance à $v_{obs} < 20 \text{ km/s}$ pour un même degré de stratification devient plus pentue avec l'âge. Cela est dû à la décélération progressive des ailes de chocs d'étrave lors de leur expansion.

Il est à noter également que comme vu dans la Sect. V.1.3.1, le gaz dans les cavités des chocs imbriqués devient de moins en moins dense au fur et à mesure que le degré de stratification du cœur augmente, ce qui se traduit dans les diagrammes MV par des signatures de masses par unité de vitesse moins élevées en passant de PCJ-U à PCJ-Z, puis à PCJ-R.

Dans l'hypothèse où seule la matière de l'enveloppe poussée est observable en CO et que le jet demeure majoritairement atomique, nous pouvons intégrer la distribution MV de la Fig. 37 afin d'obtenir la masse totale observable du flot $M = \int \frac{dM}{dv} dv$.

Nous pouvons ainsi voir qu'à $t=296 \text{ ans}$ et $t=610 \text{ ans}$, plus le degré de stratification de l'enveloppe augmente et la densité ambiante décroit, plus la masse poussée par les chocs d'étraves diminue : dans notre cas, elle descend d'un ordre de grandeur entre chaque modèle de stratification (à densité identique en $z = 0$ et $R = R_j$).

## 1.4. Comparaison au flot poussé par un vent grand angle

Nous avons ensuite effectué une quatrième simulation avec l'objectif de confronter, dans les mêmes conditions, le flot entrainé par un jet avec le flot poussé par un vent grand angle simulé par Lee et al. (2001).

Pour la comparaison entre notre simulation de jet et la simulation de vent de Lee et al. (2001), nous devons adapter deux paramètres de la simulation PCJ-R précédente, de sorte à avoir dans ce modèle-ci (nommé PCJ-RW):

- le même flux de masse, rayon d'injection, température, vitesse d'injection $v_j(t)$ et densité injectée $\rho_j(t)$ que le vent grand angle de Lee et al. (2001): Cela implique de prendre une période de variabilité du jet $P = 115 \text{ ans}$, ~3 fois plus courte que PCJ-R;

- un environnement stratifié identique à celui de Lee et al. (2001), c'est à dire un cœur aplati singulier comme dans PCJ-R, mais avec une densité de base 100 fois plus importante : $\rho_{a0} =$
<figure><img src="image_54.png" /></figure>

FIGURE 38 – Carte de densité de la simulation de jet pulsé cylindrique PCJ-RW à $t = 296$ ans, avec le même flux de masse injecté, variabilité (en vitesse et en densité) et distribution de densité environnante que dans la simulation de vent grand angle de la Fig. 12 de Lee et al. (2001). À cet âge court, nous pouvons remarquer que la taille de la couche externe qui en résulte dans PCJ-RW est similaire à celle pour le vent grand angle (montrée sur la droite).
$$1.6 \times 10^{-18} \text{ g cm}^{-3}.$$

Comme pour la Fig. 33, nous obtenons la carte de densité Fig. 38 pour PCJ-RW à t = 296 ans, qu'il est possible de comparer directement avec la Fig. 12 de Lee et al. (2001).

Nous mesurons les valeurs maximales de largeur de la couche externe induite par le jet/vent ainsi que les ouvertures angulaires à la base. Ainsi, nous mesurons des largeurs maximales de couche similaires valant $W_{max} = 2.6 \times 10^{16} \text{ cm}$ ($\sim 10R_j$) pour chaque modèle; il faut noter cependant que ces largeurs maximales sont atteintes à des hauteurs différentes. En parallèle, nous mesurons des angles d'ouverture à la base, définis comme $\alpha_{800} = 2 \arctan(W_{800}/[2 \times 800 \text{ au}])$, où $W_{800}$ est la largeur totale à hauteur $z = 800 \text{ au}$. Les angles d'ouverture obtenus valent $\alpha_{800} = 70^{\circ}$ pour PCJ-RW et $\alpha_{800} = 80^{\circ}$ pour le vent grand angle.

En somme, nous pouvons conclure qu'à cet âge court, et pour la même stratification de cœur singulier aplati "standard", favorisée par Li & Shu (1996); Lee et al. (2001); Shang et al. (2006); Wang et al. (2015b); Shang et al. (2020), le choc d'étrave externe induit par un jet pulsé s'élargit autant que la couche poussée par le modèle de vent grand angle, et forme un angle d'ouverture à la base qui peut aussi atteindre des valeurs élevées $\approx 70^{\circ}$.

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

Dans cette section, nous nous intéressons à l'influence de plusieurs paramètres libres (attendus comme étant variables parmi les sources observées) sur la morphologie et la cinématique des couches induites par un jet variable.

Cependant, contrairement aux sections précédentes et aux travaux réalisés par Lee et al. (2001), et dans un souci de plus grand réalisme, nous choisissons ici de simuler un jet conique (et non plus cylindrique). En effet, les observations de jet suggèrent des demi-ouvertures angulaires $\theta_j$ de quelques degrés aux échelles de nos simulations. Par exemple, nous avons $\theta_j = 2^{\circ}$ à une distance de 800 au dans le jet atomique RW Aur (Dougados et al. 2000); ou encore $\theta_j = 5^{\circ}$ entre 1000 et 12 000 au dans le jet CO de IRAS04166+2706 (Santiago-García et al. 2009). De ce fait, afin d'approcher le plus fidèlement l'intervalle de valeurs observées pour $\theta_j$, nous avons exploré un intervalle d'ouvertures angulaires similaires.

Pour l'enveloppe, nous gardons la configuration de cœur singulier aplati de l'Eq. (V.6), adoptée dans les modèles "standard" de flot poussés par un vent grand angle.

### 2.1. Configuration numérique

#### 2.1.1. Dispersion angulaire

De manière similaire aux travaux de Völker et al. (1999), nous introduisons une dispersion angulaire dans la valve d'entrée du jet, en maintenant un vecteur vitesse qui diverge radialement d'un point virtuel

$$(R, z) = (0, -z_0), \quad \text{with } z_0 = R_j / \tan \theta_j, \qquad (V.1)$$

où $R_j$ est le rayon du jet à une altitude $z = 0$ et $\theta_j$ est la demi-ouverture angulaire du jet.

Afin de conserver le flux de masse tout au long de la valve d'entrée du jet, le champ de densité est varié à la fois inversement au carré de la distance au point virtuel susnommé, mais aussi inversement en temps avec le module de la vitesse, c'est-à-dire tel que
$$ \rho_j(t) = \rho_{j0} \left( \frac{R_j^2 + z_0^2}{R^2 + (z + z_0)^2} \right) \frac{v_0}{v_j(t)}, \qquad (V.2) $$

où $\rho_{j0}$ est la "densité de base" du jet, c'est-à-dire la valeur de la densité dans la valve d'entrée du jet évaluée à $z = 0$, $R = R_j$, et $t = 0$. Le flux de masse constant du jet est alors donné par :

$$ \dot{M} = 2\pi (1 - \cos \theta_j) (R_j^2 + z_0^2) \rho_{j0} v_0 \qquad (V.3) $$

$$ = 2\pi (1 - \cos \theta_j) (R_j / \sin \theta_j)^2 \rho_{j0} v_0 \qquad (V.4) $$

### 2.1.2. Refroidissement radiatif et grille adaptative

À partir de cette section V.2 du chapitre V, nous employons la fonction de refroidissement radiatif $\Lambda$ de Schure et al. (2009), à hautes températures ($T \ge 10^4$ K). Cette fonction de refroidissement est une mise à jour de la fonction de MacDonald & Bailey (1981) utilisée par Lee et al. (2001), tout comme dans nos simulations de la Sect. V.1.

De plus, par rapport aux simulations de Lee et al. (2001) et des sections précédentes, où la résolution était fixée à $\Delta R = \Delta z = 14$ au, nous améliorons la résolution ainsi que la précision des simulations en activant la grille adaptative (AMR), en allant jusqu'au niveau 5 de raffinement. Cela nous donne une résolution spatiale allant de 5.2 au (résolution maximale) à 84 au (résolution minimale). Le domaine de simulation est aussi étendu jusqu'à $1.9 \times 10^{17}$ cm en $z$.

## 2.2. Grilles de modèles

Dans un premier temps, un ensemble de paramètres typiques sont définis pour construire un modèle de référence. De façon cohérente, nous gardons la même loi de vitesse d'éjection et la même distribution en densité environnante (cœur singulier aplati) que dans notre modèle PCJ-RW dans la Table 4. Si ce n'est que dans le cas ici de notre jet conique, nous adoptons une valeur de référence $\theta_j = 3^\circ$ pour la demi-ouverture angulaire du jet, ainsi qu'un rayon de référence du jet valant $R_j = 50$ au en $z = 0$.

En ce qui concerne la densité de base du jet $\rho_{j0}$, nous adoptons une valeur à "haute densité", c'est-dire telle que le flux de masse hémisphérique résultant pour le jet vaut $\dot{M} = 6 \times 10^{-6} M_\odot \text{yr}^{-1}$, pour nos valeurs de références de $\theta_j$ et $R_j$. Ce flux de masse est typique des protoétoiles de Classe 0 à forte activité.

Dans un second temps, nous lançons sept modèles où un seul paramètre à la fois est modifié par rapport au modèle de référence (dénommé $H_{REF}$ dans la suite). Ces modèles "modifiés" sont ainsi :

1. $H_{DENSA}$: avec une densité environnante de base $\rho_{a,0} = 1.6 \times 10^{-20} \text{ g cm}^{-3}$ au lieu de $1.6 \times 10^{-18} \text{ g cm}^{-3}$, soit 100 fois plus faible;

2. $H_{THETA}$: avec une demi-ouverture angulaire du jet $\theta_j = 7^\circ$ au lieu de $3^\circ$;

3. $H_{DENSJ}$: avec une densité de base du jet $\rho_{j0} = 100$ fois plus grande que dans le modèle de référence;

4. $H_{VARAMP}$: avec une demi-amplitude de variabilité $\Delta V = 90 \text{ km/s}$ au lieu de 60 km/s;

5. $H_{PER}$: avec une période de variabilité $P = 300$ ans au lieu de 115 ans;
<table><caption>TABLE 5 – Paramètres des simulations de jet pulsé conique à hautes densités, avec les ouvertures angulaires résultantes et les largeurs totales associées (Fig. 40)</caption><thead><tr><th colspan="2">Paramètres fixes</th></tr></thead><tbody><tr><td>vitesse moyenne du jet</td><td>v<sub>0</sub> = 120 km/s</td></tr><tr><td>variation en densité du jet</td><td>ρ<sub>j</sub>(t) = ρ<sub>j0</sub>[v<sub>0</sub>/v<sub>j</sub>(t)] × (R<sub>j</sub><sup>2</sup> + z<sub>0</sub><sup>2</sup>)(R<sup>2</sup> + [z + z<sub>0</sub>]<sup>2</sup>)<sup>-1</sup>, avec z<sub>0</sub> = R<sub>j</sub>/tan θ<sub>j</sub><br>(flux de masse constant avec l'altitude)</td></tr><tr><td>flux de masse du jet</td><td>Ṁ = 6 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup></td></tr><tr><td>profil de densité du cœur</td><td>Cœur singulier appli ρ<sub>a</sub>(r) = ρ<sub>a0</sub> sin<sup>2</sup> θ r<sub>0</sub><sup>2</sup> (r<sup>2</sup>)<sup>-1</sup>, avec r<sub>0</sub> = 2.5 × 10<sup>15</sup> cm</td></tr><tr><td>température du jet</td><td>T<sub>j</sub> = 100 K</td></tr><tr><td>température du cœur environnant</td><td>T<sub>a</sub> = 100 K</td></tr><tr><td>fonction de réfroidissement</td><td>Λ (100 K ≤ T < 10<sup>4</sup> K) tiré de Dalgarno & McCray (1972)</td></tr><tr><td>radiatif</td><td>Λ (T ≥ 10<sup>4</sup> K) tiré de Schure et al. (2009)<br>(R, z) = (7.0 × 10<sup>16</sup> cm, 1.9 × 10<sup>17</sup> cm) = (4679 au, 12 700 au)</td></tr><tr><td>domaine de simulation</td><td>n<sub>R</sub> × n<sub>z</sub> = 56 × 152 pour la grille complète au niveau 1 d'AMR</td></tr><tr><td>nombre de cellules</td><td>7.8 × 10<sup>13</sup> cm = 5.2 au au niveau 5 d'AMR</td></tr><tr><td>résolution maximale</td><td></td></tr><tr><td>Âge du rélevé</td><td>700 ans</td></tr><tr><th>Paramètre</th><th>Modèle de référence H<sub>REF</sub></th><th>Paramètre modifié <sup>a</sup></th><th>Nom du modèle modifié</th><th>α<sub>800</sub><sup>b</sup> (°)</th><th>W<sub>800</sub><sup>c</sup> (10<sup>16</sup> cm)</th><th>W<sub>12 700</sub><sup>d</sup> (10<sup>16</sup> cm)</th></tr><tr><td>Densité de base du cœur ρ<sub>a0</sub></td><td>1.6 × 10<sup>-18</sup> g cm<sup>-3</sup></td><td>1.6 × 10<sup>-20</sup> g cm<sup>-3</sup></td><td>H_DENSA</td><td>112</td><td>3.5</td><td>13.8</td></tr><tr><td>Demi-ouverture angulaire du jet θ<sub>j</sub></td><td>3°</td><td>7°</td><td>H_THETA</td><td>95</td><td>2.6</td><td>10.5</td></tr><tr><td>Densité de base initiale du jet <sup>e</sup> ρ<sub>j0</sub></td><td>1.8 × 10<sup>-17</sup> g cm<sup>-3</sup></td><td>1.8 × 10<sup>-15</sup> g cm<sup>-3</sup></td><td>H_DENSJ</td><td>86</td><td>2.2</td><td>8.4</td></tr><tr><td>Demi-amplitude ΔV</td><td>60 km/s</td><td>90 km/s</td><td>H_VARAMP</td><td>91</td><td>2.4</td><td>8.0</td></tr><tr><td>Période de variabilité du jet P</td><td>115 ans</td><td>300 ans</td><td>H_PER</td><td>90</td><td>2.4</td><td>7.4</td></tr><tr><td>Profil de variabilité du jet h(t)</td><td>1 + ΔV/v<sub>0</sub> sin 2πt/P</td><td>1 + ΔV/v<sub>0</sub> {1 - 2 · mod(t/P, 1)}</td><td>H_SAWT</td><td>86</td><td>2.2</td><td>7.2</td></tr><tr><td>Rayon du jet R<sub>j</sub></td><td>7.5 × 10<sup>14</sup> cm</td><td>3.0 × 10<sup>14</sup> cm</td><td>H_RAD</td><td>84</td><td>2.2</td><td>6.9</td></tr></tbody></table>

a. Dans chaque modèle modifié, seulement un paramètre à la fois est changé par rapport au modèle de référence.

b. Ouverture angulaire totale de la couche à z = 800 au. Nous obtenons α<sub>800</sub> = 88° pour le modèle de référence.

c. Largeur totale de couche à z = 800 au. Nous obtenons W<sub>800</sub> = 2.3 × 10<sup>16</sup> cm pour le modèle de référence.

d. Largeur totale de couche à z = 12 700 au (haut du domaine). Nous obtenons W<sub>12 700</sub> = 7.6 × 10<sup>16</sup> cm pour le modèle de référence.

e. Cas du jet haute densité, c'est-à-dire ayant un flux de masse hémisphérique Ṁ = 6 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup> pour les valeurs de références données pour θ<sub>j</sub> et R<sub>j</sub>.
<table><caption>TABLE 6 – Paramètres des simulations de jet conique à "densité moyenne", avec les largeurs totales et angles d'ouverture résultants</caption><caption>Paramètres différents de ceux de la Table 5</caption><thead><tr><th>flux de masse du jet</th><th colspan="2">M = 6 × 10<sup>-8</sup> M<sub>☉</sub> yr<sup>-1</sup></th><th colspan="4">(1 lobe)</th></tr><tr><th>Paramètre</th><th>Modèle de référence M_REF</th><th>Paramètre modifié<sup>a</sup></th><th>Nom du modèle modifié</th><th>α<sub>800</sub><sup>b</sup> (°)</th><th>W<sub>800</sub><sup>c</sup> (10<sup>16</sup> cm)</th><th>W<sub>12700</sub><sup>d</sup> (10<sup>16</sup> cm)</th></tr></thead><tbody><tr><td>Densité de base du cœur ρ<sub>a0</sub></td><td>1.6 × 10<sup>-18</sup> g cm<sup>-3</sup></td><td>1.6 × 10<sup>-20</sup> g cm<sup>-3</sup></td><td>M_DENSA</td><td>112</td><td>3.5</td><td>12.4</td></tr><tr><td>Demi-ouverture angulaire du jet θ<sub>j</sub></td><td>3°</td><td>7°</td><td>M_THETA</td><td>81</td><td>2.0</td><td>8.4</td></tr><tr><td>Densité de base initiale du jet ρ<sub>j0</sub><sup>e</sup></td><td>1.8 × 10<sup>-19</sup> g cm<sup>-3</sup></td><td>1.8 × 10<sup>-17</sup> g cm<sup>-3</sup></td><td>M_DENSJ<sup>f</sup></td><td>88</td><td>2.3</td><td>7.6</td></tr><tr><td>Demi-amplitude ΔV</td><td>60 km/s</td><td>90 km/s</td><td>M_VAMP</td><td>77</td><td>1.9</td><td>8.8</td></tr><tr><td>Période de variabilité du jet P</td><td>115 ans</td><td>300 ans</td><td>M_PER</td><td>74</td><td>1.8</td><td>7.4</td></tr><tr><td>Profil de variabilité du jet h(t)</td><td>1 + ΔV/v<sub>0</sub> sin 2πt/P</td><td>1 + ΔV/v<sub>0</sub> {1 - 2 · mod(t/P, 1)}</td><td>M_SAWT</td><td>74</td><td>1.8</td><td>7.0</td></tr><tr><td>Rayon du jet R<sub>j</sub></td><td>7.5 × 10<sup>14</sup> cm</td><td>3.0 × 10<sup>14</sup> cm</td><td>M_RAD</td><td>65</td><td>1.5</td><td>6.4</td></tr></tbody></table>

a. Dans chaque modèle modifié, seulement un paramètre à la fois est changé par rapport au modèle de référence.

b. Ouverture angulaire totale de la couche à z = 800 au. Nous obtenons α<sub>800</sub> = 72° pour le modèle de référence.

c. Largeur totale de couche à z = 800 au. Nous obtenons W<sub>800</sub> = 1.7 × 10<sup>16</sup> cm pour le modèle de référence.

d. Largeur totale de couche à z = 12 700 au (haut du domaine). Nous obtenons W<sub>12700</sub> = 8.0 × 10<sup>16</sup> cm pour le modèle de référence.

e. Cas du jet moyenne densité, c'est-à-dire ayant un flux de masse hémisphérique Ṁ = 6 × 10<sup>-8</sup> M<sub>☉</sub> yr<sup>-1</sup> pour les valeurs de références données pour θ<sub>j</sub> et R<sub>j</sub>.

f. Ce modèle équivaut à H_REF.
6. H_SAWT: avec un profil de variabilité en forme de dent de scie (décroissant avec le temps) au lieu d'un profil sinusoïdal (cf. Fig. 39);

7. Hiem: avec un rayon de jet à la base R_j = 20 au au lieu de 50 au.

<figure><img src="image_55.png" /><figcaption>FIGURE 39 – Comparaison entre le profil de variabilité temporelle sinusoïdal (H_REF) et en dent de scie (H_SAWTH, tel qu'adopté par Suttner et al. 1997), tous deux utilisés dans la Sect. V.2.</figcaption></figure>

Nous pouvons noter que les modèles modifiés ont la même valeur de flux de masse $\dot{M}$ du jet que dans le modèle de référence H_REF, à l'exception de H_THETA, H_DENSJ, et Hiem, puisque $\dot{M}$ varie avec $\theta_j$, $\rho_{j0}$, et $R_j$, d'après l'Eq. (V.3). Tous les paramètres d'entrée pour le modèle de référence ainsi que les modèles modifiés dans ce cas de jet "haute densité" sont résumés dans la Table 5.

Afin de vérifier la robustesse de nos conclusions, nous réalisons aussi des simulations pour une seconde grille de modèles avec les mêmes changements de paramètres, mais une densité de référence du jet 100 fois plus faible que précédemment.

Le modèle de référence (M_REF) de cette seconde série a alors un flux de masse $\dot{M} = 6 \times 10^{-8} M_{\odot} \text{yr}^{-1}$ qui est 100 fois plus faible que pour le modèle de référence de la première série (H_REF), et typique des jets à des stades davantage évolués, dans les Classes I.

Les paramètres de cette seconde série de modèles (auquel nous référons comme étant le cas à "densité moyenne") sont dans la Table 6.

## 2.3. Résultats : effet des paramètres libres

La Fig. 40 montre les différences dans la forme et la cinématique des chocs d'étraves entre le modèle de référence et chaque modèle modifié, dans le cas du jet "haute densité". La Fig. 41 montre les résultats pour la seconde série de simulations (jet à "densité moyenne"). Ces figures confrontent
<figure><img src="image_56.png" /><figcaption>FIGURE 40 – Influence de sept paramètres libres sur la géométrie et la cinématique des couches induites par un jet pulsé conique dans un cœur singulier aplati, à un âge $t = 700$ ans. Chaque ligne correspond à un différent paramètre modifié (comme indiqué dans la marge de gauche) et compare les cartes de densité en noyaux d'hydrogènes $n_H$ (gauche), de vitesse radiale $v_R$ (milieu) et de vitesse axiale $v_z$ (droite) pour le modèle de référence (moitié gauche de chaque carte) et le modèle modifié (moitié droite de chaque carte). Les valeurs du paramètre modifié correspondant sont indiquées à la base du panneau de $n_H$, tandis que les noms des modèles sont indiqués à la base du panneau de $v_z$. L'échelle de couleur des vitesses est tronquée sur l'intervalle typiquement observé dans les flots CO : entre 0 et 50 km/s. Le modèle de référence dans cette figure a un flux de masse "élevé", c'est-à-dire égal à $6 \times 10^{-6} \text{ M}_\odot \text{yr}^{-1}$ (voir Table 5 pour la liste complète des paramètres des modèles). La plus forte influence sur la forme de la cavité est ici observée lorsque la densité du cœur est modifiée.</figcaption></figure>
<figure><img src="image_57.png" /><figcaption>FIGURE 40 – Suite.</figcaption></figure>
<figure><img src="image_58.png" /><figcaption>FIGURE 41 – Même principe que la Fig. 40, mais avec un modèle de référence ayant un flux de masse du jet 100 fois plus faible : $\dot{M} = 6 \times 10^{-8} M_{\odot} \text{yr}^{-1}$.</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>FIGURE 41 – Suite.</figcaption></figure>
les modèles au même âge $t = 700$ ans, via leurs cartes respectives de densité en nucléons $n_H$, de vitesse radiale $v_R$, et de vitesse axiale $v_z$.

### 2.3.1. Cartes de vitesses

Pour les jets "haute densité" et "densité moyenne", les cartes de $v_R$ montrent l'expansion latérale de la cavité induite par les chocs d'étrave successifs. En parallèle, les cartes de $v_z$ mettent à nouveau en lumière une structure de gradient de vitesse similaire à celle d'un cisaillement, et qui se développe le long des bords de la couche principale.

### 2.3.2. Forme de la cavité principale – méthode

J'utilise les cartes de densité des Figs. 40 et 41, afin de mesurer la largeur totale $W(z)$ de la couche principale comme une fonction de l'altitude $z$. Nous obtenons alors la Fig. 42, qui superpose les formes obtenues pour les cavités de chaque modèle, à la fois pour le jet "haute densité" et de "densité moyenne". De plus, à partir de la largeur totale du choc d'étrave $W_{800}$ mesurée dans chaque modèle à $z = 800$ au, nous calculons $\alpha_{800}$, l'ouverture angulaire totale résultante près de la base du flot, suivant la définition employée par Dutta et al. (2020):

$$ \alpha_{800} = 2 \arctan(W_{800} / [2 \times 800 \text{ au}]) \quad (V.5) $$

Les valeurs résultantes de $\alpha_{800}$, $W_{800}$, et $W_{12700}$ (largeur totale mesurée au sommet du domaine) sont dans la Table 5 pour le cas "haute densité" et dans la Table 6 pour le cas "moyenne densité".

### 2.3.3. Formes de la cavité principale – analyse

Les mêmes comportements sont obtenus pour les deux grilles de modèles: jets "haute densité" et de "densité moyenne". Voici les résultats des Figs. 40 et 41:

L'effet avec le maximum d'impact à la fois sur la morphologie et la cinématique de la couche principale est ici obtenue lorsque la densité de base de l'enveloppe est diminuée d'un facteur 100. En effet, les panneaux (a) de chaque figure montrent que le choc d'étrave le plus externe est jusqu'à deux fois plus large que dans le modèle de référence. Le champ de vitesses $v_R$ atteint de plus hautes valeurs dans les chocs d'étrave (externe et imbriqués). La couche de gaz dans laquelle un gradient de $v_z$ se forme le long du choc externe (par exemple entre 0 et 50 km/s) s'épaissit aussi;

Accroître l'ouverture angulaire du jet (panneaux b) produit un effet notable mais moins important, avec une augmentation maximale (sur tous les $z$) de ~30% de la largeur totale.

Augmenter d'un facteur 100 la densité de base du jet (panneaux c) produit aussi un effet modéré, avec un élargissement maximum de 40% dans le cas le plus affecté. (comparer aussi les courbes de mêmes couleurs entre le panneau gauche et droit de la Fig. 42; et les largeurs totales compilées entre les tables Table 5 et 6).

Pour finir, la morphologie et la cinématique du choc d'étrave externe sont affectés de manière encore plus modérée lorsque :

- la demi-amplitude $\Delta V$ ou la période de variabilité $P$ sont augmentées (panneaux d et e);

- un profil en dent de scie $h(t)$ est imposé dans la valve d'entrée du jet, à la place d'un profil sinusoïdal (panneaux f);

- le rayon du jet $R_j$ est diminué (panneaux g).
<figure><img src="image_60.png" /></figure>

FIGURE 42 – Largeur totale ($W(z)$) de la couche principale à $t = 700$ ans, en fonction de l'altitude $z$, pour les modèles de jet "haute densité" (haut) et de "densité moyenne" (bas) des Figs. 40 & 41. Les courbes colorées sont obtenues en variant un seul paramètre par rapport à la simulation de référence (en pointillés et croix noires) : la densité de base du cœur $\rho_{a0}$; la demi-ouverture angulaire du jet $\theta_j$; la densité de base initiale du jet $\rho_{j0}$; la demi-amplitude de variabilité $\Delta V$; la période $P$; le profil de variabilité (en dent de scie au lieu de sinusoïdal); ou le rayon du jet $R_j$. Les valeurs du paramètre que nous faisons varier dans chaque modèle sont indiquées dans la Table 5.

Néanmoins, nous pouvons souligner qu’en comparaison du modèle de référence avec un profil de variabilité sinusoïdal, adopter un profil en dent de scie (et décroissant avec le temps) produit, dans les deux séries de modèles, une distribution plus régulière de la densité et de la vitesse dans les chocs d’étraves et les cavités du flot, qui sont moins sujets aux instabilités (cf. panneau f des Figs. 40 et 41).

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

Nous avons ensuite étudié, pour la première fois, l'évolution à long terme des flots poussés par un jet pulsé conique à travers un cœur singulier aplati, identique à celui invoqué dans le modèle "standard" de flots poussés par un vent grand angle (Lee et al. 2001). Pour ce faire, nous avons réalisé pour la première fois une simulation de ce type jusqu'à des échelles de temps atteignant 10 000 ans.

### 3.1. Choix des paramètres de la simulation

Dans la Sect. V.2, nous avons pu observer qu’adopter un profil de variabilité du jet en dent de scie (monotone décroissant par morceaux avec le temps) permet de minimiser le développement d’instabilités par rapport à un profil de variabilité sinusoïdale, sans pour autant affecter la morphologie et la cinématique qui en résultent, de manière globale.
Par conséquent, nous avons choisi pour notre simulation à temps long d'adopter un profil en dent de scie pour la variation temporelle de la vitesse. L'avantage principal est que la réduction du nombre et de l'importance des instabilités que ce choix apporte nous permet de diminuer le nombre de cellules nécessitant le niveau le plus haut du raffinement AMR, et donc d'atteindre les échelles d'âges $\simeq 10^4$ ans recherchées au niveau de cette étude avec un coût en temps de calcul plus raisonnable.

En ce qui concerne les autres paramètres, nous avons par cohérence gardé les mêmes valeurs que notre modèle de "densité moyenne" M_SAWT, excepté pour la période de variabilité où nous avons adopté une valeur de 300 ans (au lieu de 115 ans), afin de se rapprocher des valeurs typiques estimées dans le flot de HH46/47 (Zhang et al. 2019).

## 3.2. Validation du calcul de la force de rappel

La Fig. 43 montre les cartes simulées à $t = 10^4$ ans de :

(a) densité $n_H$ en noyaux d'hydrogène,

(b) fraction $f_a$ de matériau originaire du cœur environnant,

(c) et (d) vitesses radiales $v_R$ et axiales $v_z$ dans le référentiel du cœur environnant.

La sous-figure du bas (Fig. 43ii) reprend la Figure 5 de l'article Rabenanahary et al. (2022), où la force de gravité effective $\vec{F}_p = \nabla (p(t=0))$ permettant de maintenir les gradients initiaux de pression du cœur environnant a été calculée de façon numérique, exacte à l'ordre deux (méthode 1 décrite dans la Sect. IV.5.1.1 du Chapitre 4).

En comparaison, la sous-figure du haut (Fig. 43i) montre une simulation où nous avons utilisé l'expression analytique exacte du gradient de pression initial, que nous n'avions pas encore établie au moment de la publication de Rabenanahary et al. (2022), et qui est décrite dans la Sect. IV.5.1.2 du Chapitre 4.

Nous pouvons ainsi vérifier que la méthode des gradients numériques, qui est la seule méthode possible dans le cas général sans solution analytique, produit des résultats identiques (au niveau de quelques pourcents) à l'expression exacte de $\vec{F}_p$.

Nous vérifions ainsi la validité des résultats sur la morphologie et la cinématique des flots obtenus par Rabenanahary et al. (2022) avec une évaluation numérique approchée de la force de rappel $\vec{F}_p$.

Dans la suite de ce manuscript, nous présenterons les figures obtenues avec la méthode exacte, qui améliorent encore la précision de notre modélisation par rapport à Rabenanahary et al. (2022).

## 3.3. Comparaison avec les résultats aux temps courts

Les cartes de la Fig. 43i nous montrent que la structure des flots simulés à $t = 10^4$ ans diffère beaucoup de celle observée à des stades d'évolution peu avancés (comparer aux cartes Fig. 33).

### Structure et géométrie des cavités du flot

1. La couche surdense, qui formait une seule structure à $t \sim 300 - 600$ ans, se retrouve séparée à $t \sim 10^4$ ans en deux couches distinctes :

— une couche "externe" qui trace le "choc en avant" (en anglais *forward shock*) se propageant à faible vitesse dans le cœur environnant.
<figure><img src="image_61.png" /></figure>

(i) Cartes obtenues avec l'expression analytique exacte de la force de rappel $F_p$ pour maintenir le gradient de pression.

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

(ii) Cartes obtenues avec un calcul numérique à l'ordre 2 du gradient de pression (tirées de Rabenanahary et al. 2022).

FIGURE 43 – Cartes à $t = 10\,000$ ans de la (a) densité de noyaux d'hydrogène $n_H = n(H + H^+) + 2n(H_2)$, (b) fraction de gaz originaire du cœur $f_a$, (c) vitesse transverse $v_R$, et (d) vitesse axiale dans notre simulation de jet conique pulsé à travers un cœur singulier aplati. Les paramètres sont ceux de la Sect. V.3.1. L'échelle de couleur des vitesses est tronquée sur l'intervalle typiquement observé dans les flots CO : entre 0 et 50 km/s. Le contour blanc sur la moitié gauche de chaque panneau trace l'iso-contour $f_a = 50\%$.
— une couche "principale" plus interne qui trace le "choc en retour" (en anglais *reverse shock*), et dans laquelle les ailes des chocs d'étrave successifs produits par le jet variable s'empilent au cours du temps.

2. Les couches externe et principale ont toutes deux une structure "parabolique" avec un angle d'ouverture plus large à la base (eg. 86° à z = 800 au vs. 41° à z = 20 000 au pour la cavité principale).

3. Entre ces deux couches se trouve un cocon de gaz ambiant accéléré à très faible vitesse (par rapport à l'intérieur de la couche principale).

4. Nous pouvons observer que la couche principale correspond aussi approximativement à l'interface de contact délimitée par l'isocontour $f_a = 50\%$, c'est-à-dire où le gaz du cœur et du jet se mélangent à proportions égales (cf. contours blancs dans la moitié gauche des panneaux de la Fig. 43i).

La présence de deux couches de vitesses et d'angles d'ouvertures différents à la base du flot est semblable de la structure en forme "d'araignée" observée par ALMA dans B5-IRS1 (Zapata et al. 2014).

## Cartes de vitesses

Les Fig. 43(c) et (d) montrent les cartes de vitesses radiale et axiale, respectivement. Nous pouvons y voir que les chocs d'étrave générés par la variabilité du jet atteignent des vitesses plus élevées et forment des ailes plus larges qu'aux courts âges, car ils rencontrent moins de résistance. En effet, la cavité dans laquelle ils se propagent est de moins en moins dense au cours du temps, du fait de son expansion globale.

Cependant, les schémas caractéristiques de la cinématique restent similaires entre âges avancés et courts, à savoir :

1. des ailes de chocs d'étrave où le gaz atteint sa vitesse radiale maximale, avant de diminuer vers les bords de la cavité

2. une structure de cisaillement en vitesse axiale le long de la couche principale

Nous verrons dans le Chapitre VI que c'est ce schéma cinématique caractéristique qui va permettre de former une structure en forme de "cloche" dans les diagrammes position-vitesse transverses, très similaire à celle observée avec ALMA, par exemple dans HH46/47.

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

### 3.4.1. Décélération à long terme

Dans le Chapitre II, nous avons montré que les chocs d'étrave produits par des jets sont censés décélérer au cours du temps. Cependant, contrairement au modèle de flots poussés par un vent radial, aucune estimation analytique de cette décélération n'est disponible dans le cas où le cœur environnant n'est pas uniforme. Nos simulations numériques à long terme permettent, pour la première fois, de quantifier la décélération des cavités produites par un jet dans un cœur stratifié.

La Figure 44 montre l'évolution dans le temps des largeurs totales des deux couches (externe et principale) observées dans la Fig. 43i(a), à deux altitudes différentes : à la base du flot ($z = 800$ au) et au sommet du domaine de simulation ($z = 20\,000$ au), une fois que ce dernier est atteint par l'apex de la couche considérée. Ci-dessous, ces largeurs sont nommées $W_{800}$ et $W_{20000}$, respectivement.

La Figure 44 nous montre que la séparation des deux couches a lieu autour de $t = 1000$ ans pour notre choix de paramètres. Nous pouvons aussi voir que chacune de ces deux couches décélère
<figure><img src="image_63.png" /></figure>

FIGURE 44 – Evolution dans le temps de la largeur totale de la couche principale (interface jet-cœur, croix) et externe (choc en avant, cercles) pour la simulation de la Fig. 43i. Les largeurs sont mesurées aux altitudes $z = 800$ au (noir) et $z = 20000$ au (rouge). Les prédictions analytiques issues du modèle équivalent de "vent X grand angle modifié" de Lee et al. (2001) (voir Sect. V.3.4.2) sont indiquées en traits noirs/rouges pointillés avec des pentes $t^{1/2}$ (voir texte). Est aussi indiqué en gris clair l'intervalle de largeurs observées à une hauteur projetée de 800 au dans 22 flots par Dutta et al. (2020) (deuxième et troisième quartiles en gris foncé, médiane en vert). La largeur totale à $z = 20000$ au du flot HH46/47 (observé en CO par Zhang et al. 2016) est représentée en bleu.

à des âges avancés. Cette décélération est plus importante (et intervient plus tôt, et à un rythme plus rapide) à faible altitude, où la densité ambiante est plus forte. Ainsi, à $z = 800$ au, la couche
principale cesse de grossir après ~6000 ans pour atteindre une largeur totale finale autour de $W_{800} \simeq 1500$ au. En haut du domaine, par contre, à $z = 20 000$ au, elle n'atteint sa valeur finale qu'à $t \ge 10^4$ ans.

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

Les lignes pointillées tracées sur la Fig. 44 montrent les largeurs à $z = 800$ au et 20000 au prédites au cours de temps pour la cavité poussée par le modèle de vent grand angle de Lee et al. (2001), en supposant le même flux de masse injecté et la même vitesse moyenne que le jet que nous simulons ici, ainsi que le même cœur singulier aplati environnant.

Dans ce modèle dit de "vent X grand angle modifié", la vitesse du vent décroît vers l'équateur comme

$$V_w = V_{w0} \cos \theta. \qquad (V.1)$$

Cette hypothèse diffère un peu du véritable "vent X" (qui a une vitesse identique à toutes les latitudes) mais comme nous l'avons montré dans le Chapitre II, elle a l'avantage d'offrir une solution analytique simple pour la forme de la base de la cavité : une parabole, paramétrisée par $z = R^2/R_0$. Celle-ci est en expansion auto-similaire au fil du temps selon $R_0 = V_0 t$, où $V_0 = V_{w0} \sqrt{\eta}$ et $\eta$ est le constraste (fixe) de densité vent-à-cœur au niveau de l'équateur. La largeur de la parabole au temps $t$ et à la hauteur $z$ est alors simplement donnée par :

$$W_z = 2R = 2\sqrt{zV_{0t}} \qquad (V.2)$$

Pour le modèle de vent équivalent au jet considéré ici, on a :

— $V_{w0} = V_{j0} = 120$ km/s,

— $\eta = 2.2 \times 10^{-3}$ (Lee et al. 2001),

donc $V_0 = 5.6$ km/s. Cette valeur de $V_0$ est tout à fait typique de celles obtenues par ajustement du même modèle à plusieurs flots observés, $V_0 \simeq 1.5 - 9$ km/s (Lee et al. 2000; Zhang et al. 2019) (en notant que Lee et al. 2000 utilisent des notations différentes : $C = 1/R_0$ et $v_0 = 1/t$).

Les prédictions de l'Eq. (V.2) sont tracées en lignes pointillées dans la Figure 44. Elles confirment un résultat déjà noté dans la Sect. V.1.4 : aux âges courts ($\le 500$ ans), la couche produite dans un cœur singulier aplati par un vent X modifié à large ouverture angulaire a une largeur très similaire à celle créée par un jet variable de même vitesse $V_{j0} = V_{w0}$ et de même flux de masse.

Par contre, elles montrent qu'aux âges longs $t \simeq 10^3 - 10^4$ ans, l'expansion radiale auto-similaire et sans décélération de la couche poussée par le vent à large ouverture lui fait dépasser rapidement les largeurs obtenues dans notre simulation de jet où la couche interne principale du flot décélère fortement.

Tout récemment, des simulations numériques de couches poussées par un vent X non modifié (ie. sans variation de $V_w$ avec la latitude) et incluant le champ magnétique dans le vent et le cœur environnant, ont été conduites jusqu'à des âges de 500-1000 ans (Shang et al. 2020). Ces simulations confirment que la couche grossit à vitesse constante et de manière auto-similaire, à savoir avec une largeur $R_0 \propto t$ et $V \propto R$, et prédisent des largeurs similaires au modèle analytique tracé en pointillés dans la Figure 44 qui offre donc une bonne approximation du comportement attendu.
### 3.4.3. Comparaison aux largeurs de flots observées

Comme nous l'avons expliqué dans le Chapitre II, le rapport longueur/largeur typique $q \simeq 5$ cité dans la littérature pour les flots moléculaires est très probablement sous-estimé. Les raisons sont la taille trop limitée des cartographies ciblées autour d'une protoétoile donnée (voir l'exemple du flot de B5-IRS1, Arce et al. 2010) et/ou l'absence de signatures moléculaires détectables dès lors que le flot sort du nuage parent (voir l'exemple du lobe bleu "tronqué" du flot CO de HH46/47 (Zhang et al. 2016)). Pour cette raison, nous considérons que la quantité la plus pertinente dans une comparaison entre simulations et observations visant à distinguer entre les différents modèles (vent X ou jet) n'est pas le rapport longueur/largeur du flot, mais sa largeur absolue, mesurée à une hauteur donnée.

Pour effectuer cette comparaison, nous indiquons tout d'abord dans la Figure 44 par des bandes grises la répartition des largeurs totales $W_{800}$ à une hauteur (projetée) de $z = 800$ au, dans un échantillon non-biaisé de 22 flots CO, principalement de classe 0. Ces flots sont observés à la même résolution avec ALMA dans le nuage d'Orion (Dutta et al. 2020). Les premiers et derniers quartiles sont représentés en gris clair, et les deuxièmes et troisièmes quartiles en gris foncé. Nous voyons que l'évolution des largeurs $W_{800}$ de la couche principale prédites par notre simulation est en excellent accord avec l'intervalle observé par Dutta et al. (2020).

La largeur asymptotique de 1500 au atteinte à $t \ge 10^4$ ans correspond à une ouverture angulaire finale à la base de la cavité de $\alpha_{800} = 2 \arctan(W_{800}/[2 \times 800 \text{ au}]) = 86^\circ$. Ce comportement est cohérent avec une étude observationnelle sur un autre échantillon, de sources d'âges variés, qui suggère que l'ouverture angulaire à la base des flots CO cesse de grossir après typiquement $t \simeq 8000$ ans, pour atteindre une valeur finale (non-corrigée de l'inclinaison) $\alpha_{\text{obs}} \simeq 90 - 100^\circ$ (Velusamy et al. 2014).

À une altitude intermédiaire de $z = 8000$ au, la couche principale dans notre simulation atteint une largeur de ~9000 au à un âge de $10^4$ ans (voir Fig. 43i). En comparaison, l'étude récente de Habel et al. (2021) de cavités en lumière diffusée de sources de Classe 0/1 dans le nuage d'Orion avec le Hubble Space Telescope (HST) rapporte des largeurs observées à la même altitude (projetée) dans l'intervalle 1100-8500 au dans 75% des cas (cf. les demi-ouvertures angulaires mesurées par Habel et al. 2021). Là encore, nous trouvons un bon accord avec nos prédictions.

Enfin, à une altitude de $z = 20000$ au = 0.1 pc (en haut de notre domaine de calcul), la couche principale de notre simulation atteint une largeur de $W_{20000} = 15000$ au à un âge de $10^4$ ans. Il n'existe pas (encore) d'étude observationnelle des largeurs de flots à cette distance de la source et qui soit non-biaisée et uniforme en résolution. Cependant, nous notons que notre résultat est en bon accord avec la largeur observée à la même hauteur (déprojetée) du flot CO de HH46/47, qui est issu d'une source de Classe 1 et donc en principe assez évolué (Zhang et al. 2016). Cette mesure est indiquée en bleu à $t = 2 - 3 \times 10^4$ ans dans la Fig. 44.

En conclusion, contrairement à la critique formulée par Ostriker et al. (2001) dans le cas d'un environnement uniforme, le choc d'étrave formé par un jet variable lancé à travers un cœur singulier aplati apparaît comme largement capable de reproduire les largeurs et angles d'ouverture de flots CO observées à la fois à petites et grandes échelles spatiales, pour des âges réalistes $\ge 10000$ ans.

Par contre, la Figure 44 montre que, après seulement 2000 ans, la couche produite par le modèle de vent X (en pointillés) à $z = 800$ dépasse rapidement les largeurs typiques observées dans les flots CO par Dutta et al. (2020). Et cet écart augmente au cours du temps en $W \propto t^{1/2}$ (cf. Eq. V.2). Rappelons ici que l'expansion radiale sans décélération prédite dans ce modèle est étroitement liée au choix d'un profil de densité en $1/r^2$ du cœur environnant, et à l'hypothèse d'un mélange instantané à l'interface vent-cœur (Shu et al. 1991). Un vent X radial pourrait, en fait, former
des couches stationnaires si le cœur environnant présente un gradient de densité moins abrupt (eg. $1/r^{1.5}$ dans une enveloppe en effondrement) et si le mélange au niveau de l'interface vent-cœur est inefficace (Smith 1986; Liang et al. 2020). Cependant, dans ce cas, la vitesse du gaz le long des flancs de la cavité serait uniforme (au lieu de proportionnelle à la distance) et ne reproduirait donc plus la "loi de Hubble" $V \propto R$ qui a motivé ce modèle à l'origine (Shu et al. 1991), et participé à son succès. Il s'agit donc tout de même d'un problème sérieux pour le scenario classique de formation des flots CO par un vent X stationnaire.

### 3.5. Superposition rouge-bleu et vitesse maximale

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

FIGURE 45 – Diagramme PV longitudinal calculé le long de l'axe du flot simulé dans la Fig. 43i à un âge de $t = 10\,000$ ans, pour différentes inclinaisons $i$ du flot par rapport à la ligne de visée : $30^{\circ} \le i \le 90^{\circ}$ (edge-on). Les deux lobes du flots sont représentés (le lobe rouge aux offsets positifs, et le lobe bleu aux offsets négatifs).

Nous allons maintenant examiner comment notre simulation répond aux deux autres critiques majeures formulées jusqu'ici contre les modèles de flots CO poussés par les chocs d'étrave d'un jet, c'est à dire,

- une superposition d'émission rouge-bleue excessive, à la plupart des inclinaisons (Lee et al. 2001),

- des vitesses prédites trop faibles aux largeurs de flots maximales observées ($V \simeq 0.03$ km/s pour $W \simeq 10\,000$ au Arce et al. 2007).

Ainsi, la Figure 45 présente les diagrammes position-vitesse (PV) longitudinaux que nous avons
calculés à partir des résultats de la simulation de la Fig. 43 à $t = 10\,000$ ans, pour des angles d'inclinaisons $i$ du lobe bleu par rapport à la ligne de visée variant de 30° à 90° (jet vu par la tranche, ou *edge-on*). L'échelle de couleurs reflète la quantité de matière en fonction de la vitesse projetée sur la ligne de visée ($v_{obs}$, en abscisse) et de l'altitude le long de l'axe du jet, projetée sur le plan du ciel ("offset", en ordonnée).

Cette figure démontre que la présence conjointe, dans un même lobe du flot, d'émission à des vitesses positives (décalées vers le rouge) et négatives (décalées vers le bleu), dite "superposition bleu-rouge", disparait rapidement en dessous d'un certain angle d'inclinaison $i < i_c \approx 80°$ par rapport à la ligne de visée, et n'est même plus observée du tout à $i \le 60°$.

Ce comportement diffère fortement des résultats obtenus sur des temps courts de quelques centaines d'années, où la superposition bleu-rouge dans les PV est encore très importante à $i = 60°$, quelle que soit la stratification adoptée pour le cœur ambiant, comme nous l'avons montré dans la Figure 36.

<figure><img src="image_65.png" /><figcaption>FIGURE 46 – Champ d'orientation des vecteurs vitesses pour la simulation à $t = 10\,000$ ans de la Fig. 43i, dans le référentiel au repos du cœur environnant. Pour plus de lisibilité, les flèches sont montrées pour des vitesses absolues au dessus de 0.1 km/s, et ont toutes la même longueur quel que soit le module de la vitesse. En arrière-plan couleur: carte de densité $n_H$.</figcaption></figure>

La disparition de la superposition bleu-rouge aux âges longs dans notre simulation peut se comprendre grâce à l'examen de la Fig. 46, où nous représentons la direction des vecteurs vitesse à l'intérieur du flot à $t = 10\,000$ ans, dans le référentiel au repos du cœur environnant. Au niveau de la couche externe (choc en avant) on observe que le gaz choqué est en expansion quasi perpendiculaire au front de choc. Cela produit donc en principe une superposition rouge-bleu, mais l'expansion de ce font de choc est tellement ralentie aux temps longs que l'émission résultante est projetée dans le bin $v_{obs} \simeq 0$ et est donc inobservable.

En ce qui concerne la couche principale, par contre, on observe que l'interaction et l'empilement de nombreux chocs d'étrave aux temps longs donne un champ de vitesses "conique", avec des vecteurs essentiellement orientés dans la direction opposée à la source, au lieu de perpendiculaire-
ment à la couche. Cette orientation directionnelle préférentielle, combinée à la forte diminution du module de vitesse vers les bords, réduit considérablement la superposition bleu-rouge obtenue dans la projection des vecteurs vitesses du gaz, sauf aux inclinaisons très proches de "edge-on".

On note enfin que, grâce à la variabilité du jet, plusieurs chocs d'étrave "jeunes" sont présents à $z \le 20000$ au dans la cavité principale, qui continuent à produire de l'émission jusqu'à des vitesses élevées dans les diagrammes PV, même aux temps longs et dans une cavité de grande largeur.

Ainsi, les deux autres critiques avancées par les études précédentes contre les flots CO poussés par un jet, (1) trop importante superposition rouge-bleu, et (2) vitesses trop faibles aux largeurs observées, n'apparaissent plus justifiées une fois que le modèle considéré prend en compte à la fois :

1. la variabilité du jet

2. une stratification réaliste en densité du cœur environnant

3. des âges de simulation $\simeq 10^4$ ans.
<figure><img src="image_66.png" /></figure>
# CHAPITRE

## VI

# COMPARAISON AUX OBSERVATIONS ALMA DE FLOTS TYPIQUES: HH46/47 ET CARMA-7

## Sommaire

<table><tr><td>1</td><td>Aniriniil : un code Python pour les cartes synthétiques</td><td>116</td></tr><tr><td>1.1</td><td>Cube de données</td><td>118</td></tr><tr><td>2</td><td>Comparaison avec le flot de HH46/47</td><td>120</td></tr><tr><td>2.1</td><td>Morphologie</td><td>120</td></tr><tr><td>2.2</td><td>Diagrammes PV transverses</td><td>121</td></tr><tr><td>2.2.1</td><td>Résultats</td><td>123</td></tr><tr><td>2.2.2</td><td>Comparaison au modèle de vent grand angle</td><td>123</td></tr><tr><td>2.3</td><td>Diagramme PV longitudinal</td><td>123</td></tr><tr><td>2.3.1</td><td>Résultats</td><td>123</td></tr><tr><td>2.3.2</td><td>Effet de la phase du jet</td><td>125</td></tr><tr><td>2.3.3</td><td>Comparaison au modèle de vent grand angle</td><td>125</td></tr><tr><td>2.4</td><td>Distribution masse-vitesse</td><td>125</td></tr><tr><td>3</td><td>Comparaison avec le flot de CARMA-7</td><td>128</td></tr><tr><td>3.1</td><td>Observations et paramètres adoptés</td><td>128</td></tr><tr><td>3.2</td><td>Diagrammes position-vitesse</td><td>129</td></tr><tr><td>3.2.1</td><td>Comparaison pour t = 500 ans</td><td>129</td></tr><tr><td>3.2.2</td><td>Facteurs favorisant un flot jeune</td><td>129</td></tr><tr><td>3.3</td><td>Quelques réserves et alternatives</td><td>132</td></tr><tr><td>4</td><td>Synthèse</td><td>133</td></tr></table>

Dans ce chapitre, nous présentons une comparaison qualitative des résultats de nos simulations numériques avec des flots protostellaires particuliers, à travers les mêmes diagnostics que dans le chapitre précédent.

Les prédictions des simulations de la Sect. V.3 en termes de diagrammes PV, de distribution MV, et de morphologie sont qualitativement comparées avec les cartographies de ces flots par l'Atacama Large Millimeter/sub-millimeter Array (ALMA).

Les principaux résultats clefs de cette comparaison ont été publiés dans l'article Rabenanahary et al. (2022). Ils sont complétés ici par des comparaisons supplémentaires qui renforcent nos con-
clusions et confirment l'intérêt de notre modèle pour reproduire la cinématique globale complexe des flots protostellaires telle que révélée par ALMA.

Notons que l'objectif principal de notre étude est de mettre en avant d'un point de vue qualitatif les schémas/comportements/structures comparables entre les simulations de notre modèle et ce que nous retrouvons dans les observations ALMA (et de façon à ce que les similarités trouvées ressortent de manière le plus visuelle/directe possible). Nous ne ferons donc pas un ajustement fin des paramètres de simulation (valeurs de densité du coeur/du jet, vitesse d'injection du jet, âge, inclinaison du flot/du coeur...) vis-à-vis des observations.

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

Avant de procéder à la suite, je présente ici l'outil Python que j'ai intégralement programmé pour obtenir des cartes synthétiques en post-traitement des sorties du code MPI-AMRVAC, à savoir pour obtenir : des diagrammes position-vitesse et la distribution masse-vitesse du gaz. Le principe est qu'à partir de la reconstruction en trois dimensions du flot incliné (par symétrie de révolution), je construis un cube de données 3D ($\delta, \alpha, v_{obj}$), où : ($\delta, \alpha$) sont la déclinaison et l'ascension droite correspondant aux coordonnées de la projection des positions des objets 3D sur le plan du ciel ; $v_{obj}$ est la projection de la vitesse d'un objet le long de la ligne de visée d'un observateur.

La Fig. 47 résume les coordonnées employées pour réaliser les opérations du code. L'axe du flot reconstruit en 3D y est noté $z$ et est incliné d'un angle $i$ par rapport à $\Gamma$, qui est compté positivement en pointant vers l'observateur. La vitesse observée $v_{obj}$ est projetée le long de $\Gamma$. La Fig. 48 montre une représentation du cube de données inféré par l'outil Python et ensuite utilisé pour en déduire des diagrammes position-vitesse (longitudinaux et transverses) ainsi que la distribution masse-vitesse du système simulé.

Le code source est disponible sur Github : https://github.com/mrabananahary/Aniriniil.git
La commande d'exécution du code pour un fichier paramètre <nameOfParFile>.par est : pyt\nhon3 geoCompiler.py <nameOfParFile>.par.

Les paramètres donnés au code sont :

- l'angle d'inclinaison $i$ du flot par rapport à la ligne de visée comptée positivement en direction de l'observateur,

- le nom \*.vtu du fichier de sortie de MPI-AMRVAC lu par Aniriniil,

- le nombre de cellules $n_{\phi}$ suivant la direction $\phi$ de révolution autour de $z$, autour duquel le système en trois dimensions est reconstruit,

- la résolution spectrale du cube de données 3D, suivant $v_{obj}$; ainsi que l'extension choisie [$v_{min}, v_{max}$] pour $v_{obj}$ dans le cube (toutes les vitesses en km/s),

- J'ai choisi de baser et de déduire automatiquement l'extension [$\delta_{min}, \delta_{max}$] × [$\alpha_{min}, \alpha_{max}$] en ($\delta, \alpha$) du cube de données à partir de celle du fichier \*.vtu de sortie de MPI-AMRVAC, et donc de l'extension du domaine simulé, afin que tout le système reconstruit en 3D soit inclus dans le cube 3D. La résolution du cube, qui est affichée et ajustable au début de l'exécution de l'outil Python, est ainsi donné par $\delta_{max}/n_{\delta}$ et $\alpha_{max}/n_{\alpha}$, où $n_{\delta}$ est ajustable par l'utilisateur ($n_{\alpha}$ est ensuite automatiquement ajusté par le code pour avoir la même résolution en $\alpha$),

- Le nom des fichiers de sortie (diagramme PV et distribution MV sont séparés) ainsi que le choix de sortir ou non la distribution MV,

J'ai moi-même codé l'outil Python en m'inspirant largement du schéma algorithmique de MPI-AMRVAC tel que je l'ai appris, et en définissant moi-même des classes Python pour les paramètres,
<figure><img src="image_67.png" /></figure>

FIGURE 47 – Système de coordonnées employé dans l’outil Python *Aniriniil* pour l’établissement des carte synthétiques. Le plan du ciel correspond à $(\delta, \alpha)$; $z$ est confondu avec l’axe du jet reconstruit en 3 dimensions et incliné d’un angle $i$ par rapport à $\Gamma$ qui est ici compté positivement en pointant vers nous.

pour la grille 3D reconstruite, ainsi que pour le cube de données et les cartes synthétiques en
sorties du code. Ainsi, j’ai pu exploiter les techniques de programmation orientées objet faisant
<figure><img src="image_68.png" /><figcaption>FIGURE 48 – Représentation du cube de données synthétiques déduit par l’outil Python Aniriniil.</figcaption></figure>

la force de ce langage, tout en utilisant la puissance et les performances du module NumPy pour l'échantillonnage des cellules projetées sur le plan du ciel.

1.1. Cube de données

Pour l'heure, seule la masse dM des cellules 3D projetées est utilisée pour la construction du cube de données et des cartes synthétiques (à défaut de disposer de prédicitions sur l'émissivité ou l'intensité intégrée de raies d'une espèce chimique donnée) :

$$
dM = \rho(R, z, \varphi)Rd\varphi dRdz \qquad (VI.1)
$$

Nous pouvons passer des coordonnées cylindriques $\vec{r} = R\vec{e}_R(\varphi) + z\vec{e}_z$ d'une cellule aux coordonnées projetées sur le plan du ciel $\vec{r}_{proj} = \vec{r} - \Gamma\vec{e}_\Gamma = \delta\vec{e}_\delta + \alpha\vec{e}_\alpha$, via les relations :
$$ \delta = -R \cos(i) \cos(\phi) + z \sin(i) \qquad (VI.2) $$

$$ \alpha = -R \sin(\phi) \qquad (VI.3) $$

Dans la convention adoptée en astronomie, $v_{obj}$ est compté positivement en s'éloignant de l'observateur (inversement à ma convention pour le code). Dans cette convention, la projection $v_{obj}$ le long de la ligne de visée est alors donnée par :

$$ v_{obj} = -(v_R \sin(i) \cos(\phi) - v_{\varphi} \sin(i) \sin(\varphi) + v_z \cos(i)) \qquad (VI.4) $$

Enfin, une fois les coordonnées des cellules du système 3D projetées sur le plan du ciel, les masses projetées sont échantillonnées à la résolution du cube de données final. Pour finir, le cube peut ensuite être réduit suivant une direction donnée, aboutissant à un diagramme PV et/ou à une distribution MV. L'échantillonnage et la réduction utilisent la fonction Numpy `histogramdd`.

Pendant ma programmation, j'ai eu à implémenter deux méthodes de sélection et de projection du système reconstruit en 3D, ainsi que d'échantillonnage du cube de données :

1. Le premier (paramètres `voxelling_method`: 'bazooka' et `voxelling_schema`: 'brute') calcule d'abord l'ensemble des positions $(x, y, z)$ de toutes les cellules 3D reconstruites, de leurs masses respectives, puis de leurs projections sur le plan du ciel, en utilisant les opérations (multiplication, division, addition, puissance) de diffusion (en anglais `broadcasting`) hautement optimisés du module Numpy. Le but est de minimiser l'usage de boucles itératives natives de Python, très coûteuses en temps de calcul. De la même façon, j'échantillonne ensuite l'ensemble de ces projections pour aboutir au cube de données et aux cartes synthétiques;

2. Le deuxième (paramètres `voxelling_method`: 'railgun' et `voxelling_schema`: 'progressive') effectue les mêmes opérations, mais par demi-plan successif $\varphi = \frac{2\pi}{n_{\varphi}-1}$, donc à l'intérieur d'une boucle itérative Python, construisant progressivement (par pallier) le cube de données et les cartes synthétiques. De cette manière, les parties du cube/des cartes déjà traitées n'ont pas besoin d'être retenues dans la mémoire vive avant de passer à la partie suivante ;

En minimisant le nombre de boucles for sous Python, la première méthode effectue toutes les opérations en un bloc optimisé, et est ainsi la plus rapide : sur une simulation de référence avec $n_{\varphi} = 120$, $n_{\Delta} = 60$ et $i = 2\pi/3$, la durée d'exécution du code est de 19 minutes. Cependant, pour atteindre les résolutions des observations qui nous intéressent (voir Chapitre V), cette méthode a l'inconvénient de consommer énormément de mémoire vive pour évaluer, stocker et manipuler les variables contenues ou dérivées de chacune des $N_{cells} = N_{2D}n_{\varphi}$ cellules du flot reconstruit en 3 dimensions (masses, positions, vitesses). $N_{2D}$ est le nombre total de cellules 2D en sortie de MPI-AMRVAC, et peut facilement dépasser $\sim 10^5 - 10^6$ cellules allant jusqu'au niveau 5 d'AMR. En raison de cette contrainte, même lancée sur Tycho en parallèle sur le nombre maximum de cœurs, produire des cartes synthétiques avec la même résolution que dans les observations entraîne rapidement un dépassement de mémoire : jusqu'à > 100 Go demandés, là où Tycho ne peut en allouer plus de 64 Go par tâche.

C'est pourquoi j'ai codé la deuxième méthode, qui s'exécute sans problème de dépassement de mémoire, même avec la résolution élevée recherchée. Par contre, par l'ajout de boucles itératives natives de Python (très coûteuses en temps de calcul), le gain en économie de mémoire se fait en contrepartie d'une perte très significative en vitesse de calcul. En effet, pour la même simulation de référence, mais suivant la deuxième méthode, le temps de calcul passe à $2h36$, soit $\sim 8$ fois
plus lentement que la première méthode (justifiant l'intérêt en optimisation de la méthode 1 sous Python). À l'avenir, si je devais réutiliser cet outil, je devrais optimiser la première méthode pour réduire la consommation en mémoire et bénéficier de sa vitesse de calcul, pour établir des cartes synthétiques.

## 2. Comparaison avec le flot de HH46/47

Dans cette section, nous allons effectuer une comparaison qualitative "source à simulation" avec les observations ALMA du flot bipolaire HH46-47 situé à une distance $d = 450$ pc (Reipurth et al. 2000), afin de montrer que nos simulations peuvent donner des résultats non seulement prometteurs, similaires à ce qui est observé, mais aussi sur certains points plus ressemblants aux observations qu'avec les prédictions du modèle "classique" de vent grand angle.

Nous adoptons une période de variabilité $P = 300$ ans, cohérente avec les estimations inférées à partir des observations en CO de HH46-47 (Zhang et al. 2016, 2019). Ainsi, la simulation comparée aux observations dans cette section est exactement celle de la Sect. V.3.2. Nous imposons une inclinaison $i = 55°$ par rapport à la ligne de visée, en accord avec l'inclinaison moyenne du jet de HH46-47 déterminée par Hartigan et al. (2005), à partir du rapport entre la vitesse radiale des nœuds du jet et leur mouvement propre dans le plan du ciel, mesuré avec le HST. Enfin, puisque la source de HH46-47 est considérée comme proche de la transition Classe 0/Classe 1, nous considérons un âge avancé de $t \simeq 10^4$ ans.

En plus des diagrammes PV transverses et de la distribution MV, déjà présentés dans Rabenana-hary et al. (2022), nous présentons ici nos prédictions pour la carte d'émission intégrée et le diagramme PV longitudinal, qui n'étaient pas encore disponibles lors de la soumission de notre article, et qui confirment l'accord très prometteur de notre modèle avec les observations.

### 2.1. Morphologie

D'un point de vue morphologique, nous avons déjà pu voir dans la Sect. V.3.4.1 que la couche principale du flot que nous simulons parvient à atteindre une largeur comparable au flot de HH46/47 à $z = 20\,000$ au, après $t = 10^4$ans.

La Fig. 49 vient comparer l'émission CO du flot HH46/47 intégrée sur l'intervalle $2\,\text{km/s} \le |v_{\text{obs}}| \le 10\,\text{km/s}$ (Zhang et al. 2016) avec une carte synthétique de la masse (par pixel du ciel) intégrée exactement sur le même intervalle en vitesse dans notre simulation, en considérant un angle d'inclinaison $i = 55°$ et une distance de 450 pc (voir références dans Zhang et al. 2019).

Dans les prédictions de notre simulation, il est important de noter que seule la partie située à $\Delta\delta \lesssim 20''$ de la source centrale est fiable, dans la mesure où notre domaine de simulation rectangulaire est limité à une hauteur de 20 000 au, que l'apex du flot a largement dépassé à $t = 10\,000$ ans. La troncature du flot simulé combiné à son inclinaison par rapport au plan du ciel fait qu'au-dessus de $\Delta\delta \gtrsim 20''$, une importante partie du flot en trois dimensions n'est pas simulée et s'avère ainsi manquante pour la carte de masse intégrée.

Cette précaution étant prise, en comparant prédictions et observations dans la Fig. 49, il nous est possible d'affirmer que notre simulation est qualitativement capable de reproduire la morphologie observée des cavités CO observées pour HH46/47 : nous retrouvons une cavité effectivement plus ou moins conique ou parabolique en fonction de la distance à la source (et tendant vers la même largeur à $z = 20\,000$ au).
<figure><img src="image_69.png" /><figcaption>FIGURE 50 – Comparaison des diagrammes PV transverses simulés et observés dans HH46/47. Panneaux de droite : Émission en $^{12}$CO(2–1) de HH46/47 le long de deux coupes larges de 450 au à des distances projetées de 8.1 × 10$^{16}$ cm (b) et 1.00 × 10$^{17}$ cm (d) par rapport à la source centrale. Les ellipses rouges montrent les modèles analytiques de deux couches paraboliques produites par un vent X modifié et ayant le meilleur ajustement aux observations. Figures adaptées de Zhang et al. (2019). Panneaux de gauche : Diagrammes calculés à partir de la simulation de la Fig. 43i à $t \simeq 10\,000$ ans avec une inclinaison de $i = 55^\circ$ par rapport à la ligne de visée, à une distance de 450 pc. Seule la matière originaire du coeur est représentée. Les coupes ont une largeur de 400 au et sont effectuées aux mêmes altitudes projetées que les observations.</figcaption></figure>
### 2.2.1. Résultats

Chaque diagramme synthétique (panneaux (a) et (c)) a une structure caractéristique en forme de cloche, qui pique à hautes vitesses (le jet) et s'élargit plus ou moins graduellement jusqu'à la vitesse au repos du coeur ambiant. Cette forme est une conséquence directe de la décélération des ailes des chocs d'étrave lors de leur expansion à l'intérieur de la couche, ce qui produit une chute de vitesse $v_R$ et $v_z$ loin de l'axe du jet.

Les deux structures en forme de "cornes" à la base de la cloche, au voisinage de $v_{obs} \sim 0$, tracent la "couche de cisaillement" à faible vitesse située le long des parois de la couche principale, où la majorité de la masse s'empile (voir Fig. 43).

À l'intérieur de la "cloche", nous pouvons aussi discerner des ellipses brillantes à des vitesses intermédiaires. Ces ellipses tracent les intersections de la ligne de visée avec les chocs d'étrave individuels, en expansion à l'intérieur de la couche.

Ces caractéristiques (forme de cloche avec ellipses imbriquées) sont de manière frappante qualitativement proches des diagrammes PV transverses observés dans HH46/47 avec ALMA (panneaux (b) et (d) dans la Fig. 50) par Zhang et al. 2019.

### 2.2.2. Comparaison au modèle de vent grand angle

La correspondance entre nos prédicitions et les observations est d'autant plus prometteuse que nos prédicitions semblent qualitativement mieux les ajuster que les modèles de couches produites par un vent grand angle. Ainsi, les ellipses rouges dans la Fig. 50 ont tendance à surprédire la vitesse projetée du flot à grande distance de la source, en raison de la "loi de Hubble" $V \propto R$ et de la forme parabolique qui décrivent la cinématique de la couche dans ce modèle.

## 2.3. Diagramme PV longitudinal

La Fig. 51 compare le diagramme PV de HH46/47 observé le long de l'axe du flot par Zhang et al. (2019) avec les diagrammes PV de notre simulation de référence aux temps longs, incluant

1. l'ensemble du gaz originaire du jet et du coeur (Fig. 51(b))

2. ou uniquement le gaz originaire du coeur (Fig. 51(c)).

Le pas de temps représenté ici est $t = 9900$ ans (33 périodes). Ce choix sera justifié dans la section suivante.

### 2.3.1. Résultats

— En comparant les Figures 51(b) et 51(c), nous pouvons distinguer les contributions au diagramme PV associées au gaz issu du jet et à ses chocs d'étrave imbriqués successifs. Il s'agit de :

1. une structure en forme de "dent de scie" aux plus hautes valeurs de $v_{obs}$. Elle est associée au faisceau du jet, dont la vitesse oscille dans le temps,

2. une signature en "loi de Hubble" $\Delta\delta \propto v_{obs}$ partant des faibles vitesses ($v_{obs} \sim 10$ km/s) jusqu'à se confondre avec le jet à hautes vitesses. Cette contribution provient du gaz éjecté sur les côtés du jet, qui est décéléré le long des ailes du choc d'étrave.

— Ces deuxièmes signatures, traçant les chocs d'étraves internes, demeurent visibles dans le PV même lorsque l'on ne prend en compte que la matière issue du coeur (Fig. 51(c)), bien que plus faiblement que dans la Figure 51(b). Cela est dû à l'entrainement d'une certaine quantité de
<figure><img src="image_70.png" /><figcaption>FIGURE 51 – Comparaison des diagrammes PV longitudinaux observés et simulés pour HH46/47: (a) PV observé avec ALMA en $^{12}$CO(2 – 1) par Zhang et al. (2019), dans une coupe large de 450 au. (b) PV calculé pour la simulation de la Sect. V.3.2, le long d'une coupe large de 400 au à un âge de $t = 9900$ ans (33 périodes), et une inclinaison de $i = 55^\circ$ sur la ligne de visée. (c) Comme (b) mais incluant seulement la contribution du gaz issu du coeur ambiant.</figcaption></figure>
matériau ambiant dans la cavité (comme nous l'avons vu dans la carte de fraction de mélange
au Chapitre V).

— Nous observons une forte similarité qualitative entre ces structures prédites dans la Fig. 51(c) et les structures visibles dans le PV longitudinal de HH 46-47, telles qu'observées par Zhang et al. (2019) (voir Fig. 51(a)).

2.3.2. Effet de la phase du jet

Dans la Fig. 52, nous comparons les observations de HH46/47 avec deux diagrammes PV longi-
tudinaux synthétiques qui tous deux incluent uniquement la matière originaire du coeur, et diffèrent
seulement dans la phase exacte du cycle de variation de vitesse où se trouve le jet au moment de
l'observation.

Sur cette figure, nous voyons clairement que le PV de droite (panneau c), à un âge de 10⁴
ans (c'est à dire 33 périodes plus un tiers) coïncide qualitativement moins bien avec le diagramme
observé à la base du flot (en termes notamment de tangente). La comparaison aux observations
semble davantage aller dans le sens de la configuration (b) où un seul pic de $v_{obs}$ est observé dans
le champ de vue entre $\Delta\delta = 0$ et $\Delta\delta = 1.4 \times 10^{17}$ cm (les pics associés aux signatures des autres
arcs de chocs étant alors situés hors champ). Cette configuration correspond à un âge (9900 ans)
égal à un nombre entier de périodes, c'est-à-dire à une phase où le jet est actuellement à sa vitesse
minimale et sur le point de former un nouveau choc interne.

Bien sûr, cette comparaison reste qualitative et a seulement pour objectif d'illustrer l'effet de la
phase sur le PV longitudinal prédit. Un ajustement plus fin des paramètres de simulation (période
et profil de vitesse du jet) pourrait très bien trouver une correspondance encore plus similaire avec
les observations (en "comprimant" ou "étirant" en position et en $v_{obs}$ le diagramme PV obtenu).

2.3.3. Comparaison au modèle de vent grand angle

Zhang et al. (2019) a ajusté les structures en "loi de Hubble" visibles sur le diagramme PV observé de HH46/47 par deux couches paraboliques en expansion radiale $V \propto R$, suivant le modèle analytique de vent X modifié de Lee et al. (2001). Les deux paraboles ajustées (Sr1 et Sr2) sont indiquées par les courbes rouges sur la Fig. 51(a). Pour chaque parabole, la branche H ("High") représente la face arrière de la couche (projetée à plus grand décalage vers le rouge), la branche L ("Low") représente la face avant (projetée à plus basse vitesse).

Nous pouvons voir que l'accord qualitatif de notre simulation avec les observations est compa-
rable en qualité au meilleur ajustement possible par un modèle analytique de flot poussé par un
vent grand angle.

De plus, notre modèle a l'avantage d'être entièrement auto-cohérent d'un point de vue dy-
namique, ce qui n'est pas le cas de l'ajustement de HH46/47 par le modèle analytique de Lee et al.
(2001) : la présence de deux couches exigerait que le vent grand angle ait des épisodes éruptifs
rapprochés d'environ 300 ans; mais dans cette situation, il semble difficile que la deuxième couche
se propage elle-aussi dans une loi de densité ambiante en $\sin^2\theta/r^2$, nécessaire pour obtenir une
parabole en expansion homothétique.

2.4. Distribution masse-vitesse

Avec ALMA, Zhang et al. (2016) ont pu inférer la distribution masse-vitesse de chacun des
deux lobes de HH46/47, à partir de son émission en ¹²CO, ¹³CO et C¹⁸O. Dans la Fig. 53, nous
<figure><img src="image_71.png" /><figcaption>FIGURE 52 – Comme la Fig. 51 mais cette fois les deux diagrammes PV longitudinaux synthétiques (b) et (c) incluent uniquement la matière originaire du coeur, et diffèrent seulement dans la phase du cycle de variabilité où le jet est observé: 9900 ans (exactement 33 périodes) pour (b) et $10^4$ ans (33+1/3) pour (c).</figcaption></figure>
<figure><img src="image_72.png" /></figure>

FIGURE 53 – Comparaison de la distribution masse-vitesse (MV) observée dans HH46/47 et simulée. *Symboles noirs connectés: MV* de la matière originaire du coeur, évalué à *t* ≃ 10⁴ ans à partir de la simulation de la Fig. 51, incliné de *i* = 55° par rapport à la ligne de visée. *Symboles bleus et rouges* : distributions MV dans les lobes bleu et rouge du flot HH46/47, dérivées par Zhang et al. (2016). Les carrés vides supposent une émission en CO optiquement mince, tandis que les cercles pleins incluent une correction d'épaisseur optique dépendant de la vitesse (voir texte). Les lignes tiretées montrent les pentes des lois de puissance *γ* = -3.4 (rouge) et *γ* = -2.7 (bleu).

comparons ce diagramme MV à celui prédit dans notre simulation pour le gaz entrainé originaire du coeur, à un âge $t \simeq 10\,000$ ans et en supposant le jet incliné de $i = 55^\circ$.

Nous y représentons les deux distributions MV estimées par Zhang et al. (2016) :

1. les symboles ouverts ont été calculés à partir de l'émission en $^{12}$CO sans correction de l'épaisseur optique. La faible valeur de la pente de loi de puissance $\gamma \simeq -2$ observée est typique des valeurs d'exposants précédemment mesurées (avant correction d'épaisseur optique) dans d'autres flots CO;

2. les symboles pleins ont été calculés après application d'une correction d'épaisseur optique dépendante de la vitesse. Cette correction a été évaluée sur la base des données en $^{12}$CO, $^{13}$CO, ainsi que $^{18}$CO (qui est optiquement mince à basses vitesses). Sachant que l'épaisseur optique augmente à faibles vitesses, la distribution MV inférée a une pente plus abrute d'exposant $\gamma \simeq -3.4$ (en échelle log-log). Des valeurs de pentes aussi abruptes ont aussi été inférées dans d'autres flots, après application de la même correction d'épaisseur optique (par exemple dans B5-IRS1, par Yu et al. 1999).
La Fig. 53 nous montre que notre simulation s'accorde remarquablement bien avec les distributions MV observées dans les deux lobes après correction d'épaisseur optique (qui devraient être les plus proches de la distribution en masse réelle du flot). Cet accord est moins convaincant pour le lobe bleu que le rouge, ce qui pourrait être l'effet d'une asymétrie des conditions en densité/stratification de l'enveloppe dans les deux lobes

Cette correspondance avec nos simulations renforce davantage l'attractivité du modèle de chocs d'étrave d'un jet variable lancé jusqu'à des âges avancés à travers un environnement fortement stratifié en pression, comme possible origine des flots CO.

## 3. Comparaison avec le flot de CARMA-7

Dans cette section, nous nous tournons vers un flot considéré comme très proche d'une inclinaison "edge-on", et plus jeune (Classe 0) que HH46/47 : CARMA-7. Nous illustrerons ainsi le biais très important qu'une forte inclinaison peut introduire dans les estimations de période de variabilité du jet à partir des diagrammes PV. Nous montrerons aussi en quoi les PV observés favorisent l'idée d'une source particulièrement jeune, puis discuterons d'une explication alternative qui impliquerait une précession du jet.

### 3.1. Observations et paramètres adoptés

CARMA-7 est une protoétoile très enfouie située dans la région active de formation stellaire "Serpens South", à 415 pc de nous. Elle fait partie d'un amas dense de sources très jeunes (Classe 0), dont elle est la plus brillante.

De récentes observations ALMA du flot CARMA-7 à haute résolution et haute sensibilité, par Plunkett et al. (2015), ont permis d'obtenir un diagramme PV longitudinal de grande qualité, et allant jusqu'à une distance (projetée) de $1.5 \times 10^{17}$ cm = 10 000 au (la moitié de notre domaine de simulation).

Ce diagramme PV observé est reproduit dans la Fig. 54(a). Comme noté par Plunkett et al. (2015), il exhibe une superposition bleu-rouge très significative. De ce fait, il est très probable que le flot soit incliné tout près du plan du ciel (voir Sect. V.3.5). Après quelques tests préliminaires sur les diagrammes PV prédits à différents angles pour notre simulation de référence, la valeur d'inclinaison de $i = 85°$ par rapport à la ligne de visée estimée par Plunkett et al. (2015) est celle que nous retiendrons pour comparer nos modèles à CARMA-7.

Comme noté par Plunkett et al. (2015), le diagramme PV observé de CARMA-7 exhibe aussi, sur toute sa longueur, une série quasi-périodique de 22 pics de vitesse, qui suggère une variabilité épisodique de la vitesse d'éjection. En calculant les âges dynamiques de chaque épisode selon $\tau_{dyn} \simeq D/v_{obs}$, avec D la distance de chaque pic et $v_{obs}$ sa vitesse maximale projetée, l'écart entre les pics de vitesses suggère une période de variabilité de $\Delta\tau_{dyn} \simeq 310 \pm 150$ ans, avant correction de l'inclinaison (Plunkett et al. 2015).

Cependant, la Fig. 45 du Chapitre V montre que dans notre simulation de période $P = 300$ ans, l'écart spatial entre les positions de deux pics de vitesse successifs à $i \simeq 85°$, est beaucoup plus important que l'écart entre pics observé dans le PV de CARMA-7 (Fig. 54(a)).

Ce désaccord provient des très importants effets de projection de vitesse qui affectent les âges dynamiques dans un jet vu par la tranche. Avec $v_{obs} \simeq \cos i V_j$, et $\Delta D \simeq \sin i \lambda$ (où $\lambda$ est le véritable espacement entre les nœuds avant projection) la valeur de $\Delta\tau_{dyn}$ devrait être corrigée d'un facteur ($\cos i/\sin i) \simeq 1/10$ pour $i \simeq 85°$ (Plunkett et al. 2015). De ce fait, pour comparer qualitativement
nos simulations avec les observations de CARMA-7, nous avons lancé une autre simulation à long
terme où la période de variabilité est dix fois plus courte que dans la simulation de référence, à
savoir : *P* = 30 ans.

Les diagrammes PV longitudinaux prédits par cette nouvelle simulation, pour *i* = 85°, sont comparés au PV observé de CARMA-7 dans les sections suivantes. Ils sont tous sommés sur une largeur de 419 au autour de l'axe du jet, similaire à la résolution angulaire des observations ALMA de Plunkett et al. (2015).

## 3.2. Diagrammes position-vitesse

3.2.1. Comparaison pour $t = 500$ ans

La Fig. 54 nous montre une comparaison qualitative du diagramme PV de CARMA-7 avec les diagrammes PV prédits pour notre modèle avec *P* = 30 ans, à un âge jeune *t* = 500 ans (qui sera justifié dans la section suivante). Deux versions du diagramme sont montrées :

— la Fig. 54(b) inclut l'ensemble du gaz, originaire du jet ou du coeur,

— la Fig. 54(c) n'inclut que le gaz initialement originaire du coeur.

Les résultats principaux sont les suivants :

1. Le diagramme PV de la Fig. 54 (b) est dominé par des structures ressemblant à des dents de scie, qui tracent le jet variable dans le temps ainsi que le matériel éjecté sur le côté à partir des chocs internes au jet, produits par sa variabilité en vitesse,

2. Le diagramme PV de la Fig. 54 (c) est quant à lui dominé par des structures montrant une accélération apparente de type "loi de Hubble" (Li & Shu 1996; Lee et al. 2001). Ces structures tracent le gaz ambiant repoussé et accéléré par les ailes des chocs d'étrave successifs. Les structures en dent de scie du jet demeurent visibles, mais beaucoup plus faiblement que dans la Fig. 54(b)). Ces signatures rémanentes tracent le gaz du coeur ambiant qui a été incorporé à l'intérieur de la cavité et entraîné le long des parois du jet (voir aussi la carte de mélange dans Rabenanahary et al. 2022),

3. Une comparaison révèle que les structures prédites dans la Fig. 54(c) présentent qualitative-ment une forte similarité avec les observations en CO de CARMA-7 (Fig. 54(a)), y compris l'espacement entre les pics de vitesse projetée.

En tenant compte du fait que notre modèle effectue un certain nombre de simplifications, cet accord qualitatif avec les observations est très encourageant. Il confirme la nécessité de corriger de façon adaptée les périodes de variabilité directement estimées à partir des observations, afin de tenir compte de l'inclinaison du flot. Dans le cas des flots quasiment observés par la tranche comme CARMA-7, l'effet peut être très important, avec un facteur 10 ici sur la période.

3.2.2. Facteurs favorisant un flot jeune

Ici nous justifions notre choix de comparer les observations de CARMA-7 avec une simulation
à un âge jeune de 500 ans, plutôt que *t*~10 000 ans (comme pour HH46/47).

Tout d'abord, notons que dans son étude de CARMA-7, Plunkett et al. (2015) trouvent un âge dynamique de 5000-6000 ans pour les pics de vitesse observés à 24'' = 10,000 au de la source. Après application de la correction d'inclinaison (cosi/sini) ≃ 1/10 pour i = 85°, leur âge réel d'éjection serait donc de 500-600 ans seulement.

Nous avons donc calculé une série de diagrammes PV évalués à des âges différents de notre
simulation: 500, 1000, 3000, 5000, 8000, et 10 000 ans. Nous avons trouvé une évolution claire des
<figure><img src="image_73.png" /><figcaption>FIGURE 54 – Comparaison des diagrammes PV longitudinaux observé dans CARMA-7 et simulés: (a) PV observé avec ALMA en $^{12}$CO(2-1) par Plunkett et al. (2015) dans une coupe large de 374 au. (b) PV calculé à un âge de $t = 500$ ans pour une simulation de période $P = 30$ ans (voir Sect. VI.3.1) le long d'une coupe large de 400 au, et pour un angle d'inclinaison de $i = 85°$ par rapport à la ligne de visée. (c) Comme (b) mais incluant seulement la contribution de la masse issue du coeur.</figcaption></figure>
<figure><img src="image_74.png" /><figcaption>FIGURE 55 – Comparaison du PV longitudinal observé de CARMA-7 (a) avec deux diagrammes PV synthétiques incluant uniquement la matière originaire du coeur, et différant par l’âge de la simulation: (b) $t = 500$ ans, (c) $t = 10\,000$ ans.</figcaption></figure>
diagrammes PV avec l'âge, qui est illustrée dans la Fig. 55 où nous comparons les diagrammes PV prédits à deux âges situés aux extrémités de notre échantillon : $t = 500$ ans (Fig. 55(b)) et $t = 10\,000$ ans (Fig. 55(c)).

Les diagrammes PV de la Fig. 55 nous montrent que la simulation à l'âge plus court ressemble qualitativement davantage aux observations, comparé à l'âge le plus avancé.

1. Plus l'âge de nos simulations augmente, plus les pics de vitesses ont tendance à s'aligner uniformément le long d'une même valeur $v_{obs} \simeq 35$ km/s (voir Fig. 55(c)). Celle-ci correspond à la projection sur la ligne de visée de la vitesse axiale moyenne du jet, additionnée à la vitesse transversale à laquelle le matériau du jet comprimé dans les chocs internes est éjecté sur les côtés par les forces de pression. Aux âges courts (300-500 ans), par contre, où la cavité du flot est encore de faible largeur, le gaz ainsi ejecté est contraint rapidement à décélérer à l'approche des parois, et ne peut donc pas atteindre la vitesse maximale $v_{obs} \sim 35$ km/s. Cela est d'autant plus vrai que l'on s'approche de la tête du jet, où la cavité est plus étroite, et il en résulte donc une diminution de la vitesse maximale des pics avec la distance à la source dans les flots très jeunes (voir Fig. 55(b)). Cette décélération apparente est aussi visible dans les observations de CARMA-7, qui semblent donc en meilleur accord avec notre simulation âgée de 500 ans qu'à 10 000 ans.

2. Au cours du temps, les signatures en forme d'éperon des premiers chocs d'étrave ont tendance à décélérer progressivement vers la vitesse au repos $v_{obs} \simeq 0$ km/s. Les structures en accélération apparente de type "loi de Hubble", bien visibles à $t = 500$ ans (Fig. 55(b)) deviennent ainsi de plus en plus étroites avec l'âge (en intervalle de $v_{obs}$), et indiscernables à $t = 10\,000$ ans, pour $i = 85^\circ$ (Fig. 55(c)). De ce point de vue, aussi, les observations de CARMA-7 semblent en meilleur accord avec nos simulations à 500 ans que $10^4$ ans.

Ces comparaisons nous montrent que nos simulations d'un jet dans un coeur singulier aplati vont qualitativement dans le sens d'un flot très jeune (moins de 1000 ans), en accord avec l'âge dynamique estimé par Plunkett et al. (2015) après correction d'un facteur 1/10 pour l'effet d'inclinaison.

## 3.3. Quelques réserves et alternatives

Cependant, un âge aussi jeune que 500 ans pose plusieurs difficultés.

D'une part, il semble assez peu probable d'observer une protoétoile moins de 1000 ans après le début du lancement du jet rapide, alors que la phase de Classe 0 dure quelques $10^4$ ans. D'autre part, pour obtenir une décroissance de la vitesse des pics avec la distance à la source dans notre simulation, la tête du jet doit être très proche de l'extrémité du PV (car c'est l'étroitesse de la cavité à son apex qui ralentit l'expansion latérale du gaz). Mais alors, le diagramme PV inféré possède à chacune de ses extrémités une structure en éperon (correspondant au choc d'étrave à la tête du jet), qui n'est pas observée (comparer Fig. 55(a) et (b)). En réalité, la carte de la Figure 1 de Plunkett et al. (2015) montre que le flot de CARMA-7 s'étend bien au delà des $\pm 24''$ couverts par leur diagramme PV, au moins jusqu'à $\pm 40''$, donc le ralentissement des pics de vitesse avec la distance ne semble pas associé au rétrécissement de la cavité près de son apex (comme c'est le cas dans notre simulation à $t = 500$ ans).

Une hypothèse alternative pour reproduire la décroissance de vitesse des pics sur une échelle plus courte que l'extension totale du flot, serait la précession du jet. Si la période de précession dépasse largement les 30 ans, le changement graduel d'inclinaison du jet (plus proche du plan du ciel loin de la source) induirait une diminution graduelle de la vitesse maximum projetée avec la distance. L'inclinaison proche de "edge-on" de CARMA-7 est particulièrement favorable à l'observation de
cet effet, déjà observé dans d'autres flots edge-on comme L1157 (Kwon et al. 2015; Podio et al. 2016) ou HH212 (Lee et al. 2017). En plus de cet effet de projection de la vitesse, un autre effet de la précession du jet serait que les zones de chocs internes du jet (produisant les pics de vitesse) ne se propageraient plus à vitesse constante mais seraient freinées aérodynamiquement par le gaz du cocon, car elles ne seraient plus "poussées par derrière" par les ejecta plus récents (le jet ayant changé de direction entre temps). Cet effet a déjà été invoqué pour reproduire la décélération intrinsèque des chocs d'étrave successifs avec la distance observée dans le jet de Classe 1 de HH34 (Cabrit & Raga 2000).

Nos simulations MPI-AMRVAC étant axisymétriques, nous n'avons pas pu tester plus pré-
cisément l'effet de la précession dans ce travail, mais nous la considérons comme une alternative
attractive pour expliquer la décélération apparente des pics de vitesse dans le PV de CARMA-7,
avec des âges plus typiques des objets de Classe 0.

4. Synthèse

Bien que nous n'ayons pas cherché à opérer un ajustement fin des paramètres libres à partir
d'une grille de simulations, nous pouvons être satisfaits de voir à quel point notre simulation de
référence est parvenue à produire qualitativement des prédictions comparables à ce qui est observé
dans plusieurs flots particuliers (HH46/47 et CARMA-7). De ce fait, même si nous reconnaissons
que la coincidence entre notre simulation et les observations est loin d'être parfaite et sans défauts,
les similarités qualitativement relevées montrent à quel point notre modèle très simplifié de jet
conique dans un coeur singulier stratifié s'avère prometteur pour des études plus poussées à venir
(telles qu'un ajustement plus fin entre les sources et les simulations comparées).
<figure><img src="image_75.png" /></figure>
# CHAPITRE

## VII

# CONCLUSION GÉNÉRALE, RÉSERVES ET PERSPECTIVES

### Flots poussés par un jet dans une enveloppe stratifiée en densité

Nos simulations numériques nous ont permis de confirmer que les chocs d'étraves poussés par un jet cylindrique grossissent plus vite à travers une enveloppe environnante, avec une densité qui décroît de manière abrupte, que dans une enveloppe uniforme. Cette dernière a jusqu'à présent été la configuration communément adoptée dans les modèles analytiques et les simulations de flots poussés par un jet. La différence dans les largeurs totales de flots peuvent atteindre un facteur 8. Plus le degré de stratification de l'enveloppe augmente, plus les chocs d'étraves s'élargissent rapidement.

### Largeur comparée au modèle de vent grand angle (à temps courts)

À partir de ce résultat, nous avons pu montrer dans nos simulations que pour des âges de quelques centaines d'années, le choc d'étrave externe poussé par un jet peut atteindre la même largeur totale maximale que pour un flot poussé par un vent X grand angle, tout en adoptant : la même stratification de l'enveloppe (à savoir celle d'un cœur singulier aplati), le même flux de masse injecté, et la même variabilité en vitesse pour le jet.

### Influence des paramètres libres

Par la suite, nous avons étudié l'influence de plusieurs paramètres libres dans notre modèle, sur la morphologie générale, l'angle d'ouverture, et la cinématique des chocs d'étraves poussés. Le paramètre qui a de loin le plus d'impact est la densité de l'enveloppe environnante, suivie de la densité du jet et de l'angle d'ouverture du jet. Dans les intervalles de valeurs explorés durant cette thèse, les autres paramètres libres que sont l'amplitude en vitesse de la variabilité, le profil de cette dernière, sa période et le rayon du jet n'affectent pas significativement ces résultats de nos simulations.

### Influence de l'âge et de prédicitions à 10⁴ ans

Enfin, pour la première fois, nous avons lancé des simulations de jet conique pulsé lancé jusqu'à des âges de *t* = 10 000 ans, à travers une enveloppe stratifiée cohérente avec les observations, avec une résolution maximale de ~5 au, et sur une échelle de 20 000 au. Les résultats de ces simulations montrent des changements drastiques dans les flots simulés à long terme. En effet, au bout de plusieurs centaines d'années, la structure de choc d'étrave externe initiale se divise en deux sous-structures distinctes : un choc d'étrave externe parabolique et lent composé majoritairement de gaz originaire de l'enveloppe environnante (et traçant le choc en avant); et un choc d'étrave interne, dit principal, plus rapide, traçant (au sein du flot poussé) l'interface entre le gaz originaire du jet
et celui venant de l'enveloppe. L'expansion de ce choc interne décélère et s'arrête à la base du flot après ~8000 ans, contrairement aux flots poussés par un vent grand angle.

Le choc d'étrave principal délimite les contours d'une cavité au sein de laquelle plusieurs chocs d'étrave successifs poussés par le jet pulsé sont en expansion pour ensuite ralentir en raison de leur interaction avec les chocs d'étraves imbriqués précédents. Cette décélération produit une structure de cisaillement abrupt en vitesse, parallèle aux parois (à peu près coniques) du choc d'étrave principal. À la fois la morphologie et les champs de vitesse dans nos simulations à temps longs sont très différents des prédictions analytiques des modèles balistiques de chocs d'étrave poussés par un jet à travers un milieu uniforme Ostriker et al. (2001).

## Résolution des critiques du modèle de flots poussés par un jet

Les simulations à temps longs de notre modèle basique ne montrent aucun des défauts et limites avancées comme critiques des modèles de flots poussés par des jets stationnaires dans un milieu uniforme : rapport d'aspect et collimation excessive, superposition bleu-rouge excessive, et vitesses trop faibles pour de grandes largeurs de flot.

## Des flots réalistes qualitativement proches des observations

Au contraire, elles montrent des similarités prometteuses avec les observations les plus récentes de flots de Classe 0, observés à haute résolution par ALMA. Il s'agit de similarités en termes de largeurs des chocs d'étraves prédits, d'angle d'ouverture à la base du flot (≈ 90° à 10⁴ ans), de diagrammes position-vitesse longitudinaux et transverses, et de distribution masse-vitesse (Dutta et al. 2020; Plunkett et al. 2015; Zhang et al. 2016, 2019).

Certaines de nos comparaisons montrent même des similarités aux observations encore plus proches qu'avec les prédictions du modèle de "vent X modifié" de Lee et al. (2001) largement employé à ce jour. C'est le cas de la largeur totale (modérée) du flot à 800 au, de la forme caractéristique de "cloche" dans les diagrammes PV transverses, ainsi que la relation en loi de puissance abrupte prédite entre masse et vitesse, après correction de l'épaisseur optique en CO.

De manière plus générale, les résultats de la thèse ont montré que non seulement une modélisation réaliste de la stratification en densité de l'enveloppe entourant le jet, mais aussi des temps d'intégration d'au moins 10⁴ ans, sont essentiels pour prédire avec fiabilité les propriétés de flots poussés par un jet pulsé, ainsi que pour confronter les prédictions de ce modèle avec les observations.

## Mise en garde et perspectives

### - Refroidissement

Nous pouvons aussi discuter brièvement de la physique mise de côté dans nos simulations numériques, ainsi qu'une estimation de la manière dont la prise en compte de celle-ci pourrait influencer nos résultats. Dans un premier temps, nous avons considéré une fonction de refroidissement radiatif à l'équilibre, c'est-à-dire sans suivre l'évolution de l'ionisation et du refroidissement hors équilibres derrière les fronts de choc. Par conséquent, nous n'étions pas en mesure de calculer des cartes d'émission synthétique fiables et réalistes dans les raies atomiques et ioniques, dont le flux émis est proportionnel à la densité électronique. Cette approximation est cependant suffisante, en considérant que notre objectif est de déterminer la forme globale, la dynamique et la distribution de masse dans les enveloppes denses composées de gaz post-choc refroidi. Nous avons aussi négligé les processus chimiques, en sachant que la dissociation du H₂ n'affecte que la distribution de masse observée au-dessus de 20 km/s (Downes & Cabrit 2003; Moraghan et al. 2008), tandis que nos comparaisons aux observations se font à des vitesses plus faibles. De telles approximations
permettent de réaliser pour la première fois, à moindre coût ressources numériques, des simulations
à long terme, atteignant $10^4$ans, et pouvant être comparées aux flots observés avec ALMA.

- Effondrement de l'enveloppe

Dans l'ensemble, nous avons négligé l'effondrement de l'enveloppe au cours de son interaction avec le jet, sous l'influence du champ gravitationnel protostellaire. Cela nous a permis de comparer correctement nos simulations avec les simulations de flots poussés par vents X grand angle. Ces dernières sont lancées avec une stratification d'enveloppe sous forme de cœur singulier aplati (Lee et al. 2001), et dont l'effondrement est également négligé. Dans de telles enveloppes, la région mise en mouvement par l'effondrement gravitationnel grandit et se propage de l'intérieur vers l'extérieur de l'enveloppe à la vitesse du son a atteinte dans cette dernière. Cet effondrement a pour effet de générer une distribution en densité moins abrupte, en $1/r^{1.5}$, contenue à l'intérieur du rayon sonique, $r_{inf} = a \times t$ (Terebey et al. 1984), où t est l'âge depuis le début de l'expansion de la région d'effondrement. Contrairement aux flots poussés par un vent où cet effondrement peut conduire à une couche stationnaire (Liang et al. 2020), une pente en $1/r^{1.5}$ en densité ne réduit que faiblement la largeur totale des flots poussés par un jet, par rapport à une densité d'enveloppe en $1/r^2$ (Moraghan et al. 2008). L'effet principal attendu de l'effondrement de l'enveloppe, en contact avec le flot développé, serait d'ajouter du cisaillement, de l'entraînement et de la pression dynamique supplémentaire à la base du choc d'étrave externe (cf. Liang et al. 2020, pour le cas d'un flot poussé par un vent). Pour nos paramètres de simulation, l'effondrement se produirait à l'intérieur de $r_{inf} = 1000$ au × (a/0.5 km/s) × (t/104 ans) et, par conséquent, n'affecterait qu'une petite fraction de notre domaine de simulation s'étendant jusqu'à 20 000 au.

- Période de variabilité

Nos paramètres de simulation n'explorent pas le régime de "longue période de variabilité", où l'enveloppe ambiante aurait le temps de se reconstituer partiellement entre les épisodes de sursauts liés à la variabilité du jet. Une telle situation intéressante pourrait se produire si les sursauts majeurs produisant les chocs d'étrave imbriqués ne se produisaient que tous les quelques $10^4$ ans, par exemple en raison de l'interaction de marées au sein d'un système de protoétoiles binaires avec une grande excentricité. L'espacement observé entre les nœuds du jet dans les sources de Classe 0/1 indique cependant des échelles de temps de variabilité beaucoup plus courtes, avec de multiples modes de périodes de quelques décennies, de quelques siècles, voire de quelques milliers d'années (Raga et al. 2012; Lee 2020). Le régime de "longue période de variabilité" semble assez rare chez les Classes 0/1 enfouies les plus jeunes et pour lesquelles des flots CO sont observables.

- Champ magnétique

Enfin, nous avons négligé la contribution des champs magnétiques, tant dans le jet que dans l'enveloppe ambiante. La pression magnétique dans le jet agirait pour réduire le facteur de compression atteint dans le post-choc tout en augmentant la longueur de refroidissement du gaz, typiquement d'un ordre de grandeur (Hartigan et al. 1994). Cet effet est observé directement dans les chocs internes résolus de jets stellaires (Hartigan & Wright 2015). Des simulations numériques de jets magnétisés à une haute résolution spatiale et avec refroidissement radiatif montrent que l'affaiblissement du refroidissement radiatif tend à élargir la largeur des chocs d'étraves imbriqués. Cependant, l'effet apparaît comme plutôt modeste (de Colle & Raga 2006), justifiant son omission dans nos simulations.

À l'inverse, un champ magnétique dans l'enveloppe environnante a tendance à résister à l'expansion latérale du choc d'étrave externe. Une diminution importante de la largeur de ce choc par rap-
port au cas hydrodynamique nécessite cependant des champs magnétiques intenses, près du point d'équipartition (c'est-à-dire avec un rapport $\beta \simeq 1$ entre la pression thermique et la pression magnétique). Cela est par exemple démontré dans les simulations de Shang et al. (2020), pour des flots poussés par des vents X dans des enveloppes avec différents degrés d'aplatissement et donc de magnétisation<fnref n="1" />. Pour une valeur modérée standard $n = 2$ du degré d'aplatissement de l'enveloppe, tel qu'adopté dans le présent travail, la réduction maximale de la largeur du choc d'étrave externe due au champ magnétique ambiant ne devrait être que de 25 %. Par conséquent, la forme et la dynamique des cavités que nous avons prédites devraient rester valides. Nous notons qu'un degré de complexité supplémentaire pour une enveloppe avec un degré de magnétisation plus élevé, dans le cas général d'une rotation non-nulle, permettrait probablement le lancement probable d'un vent magnétohydrodynamique massif et lent depuis le disque képlérien formé autour de la protoétoile (voir par exemple Lesur 2021 et les références qui y figurent). L'interaction d'un jet interne pulsé (ou d'un vent X grand angle) avec un vent de disque externe pourrait considérablement affecter la formation des cavités du flot, mais l'évolution à long terme n'a jusqu'à présent été examinée que dans le cas hydrodynamique, à notre connaissance (Tabone et al. 2018).

## - Chimie du gaz

Concernant les observations, nous nous sommes ici concentrés sur les signatures des chocs d'étraves ou couches en termes de formes, de cinématique et de distribution de masse à des vitesses inférieures à 20 km/s, qui sont tracées par l'émission en CO observable avec ALMA. Nos simulations peuvent également être utilisées pour évaluer la détectabilité du matériau originaire du jet et situé à de plus hautes températures. À partir de la carte de vitesse de la Fig. 43, nous pouvons estimer que les chocs d'étrave imbriquées (entraînées par chacune des pulsations du jet) génèrent des vitesses de choc typiques de ~30 km/s dans la cavité précédemment creusée. La densité pré-choc à l'intérieur de la cavité est $\simeq 100 \text{ g cm}^{-3}$ à $z \le 3000$ au, et chute à des altitudes plus élevées. En utilisant la grille du modèle de choc atomique de Hartigan et al. (1994), la brillance de surface maximale prédite en [S II] pour des abondances solaires serait alors de $\simeq 5 \times 10^{-6} \text{ erg cm}^{-2} \text{ sr}^{-1}$. L'observation d'une telle émission étendue en [S II] à l'intérieur des cavités du flot serait pertinente à rechercher, mais elle resterait difficile à isoler par rapport au jet axial bien plus brillant en ainsi que vis-à-vis de l'émission en lumière diffuse des parois de la cavité, d'autant plus que les jeunes protoétoiles avec de puissants flots brillants en CO sont souvent situées dans des régions avec une forte épaisseur optique.

Lorsque le jet est suffisamment dense pour être en partie moléculaire, un autre traceur possible pour les couches emboîtées observées à l'intérieur de la cavité de flot principale est l'émission ro-vibrationnelle en $H_2$. L'exemple le plus spectaculaire est de loin le flot de Classe 0 HH 212. Chaque choc d'étrave large observé en $H_2$ est connectée à un choc d'étrave observé séparément en CO à la base du flot. Ce choc d'étrave est observé comme étant produit par une pulsation majeure du jet, et est emboîté à l'intérieur de la cavité principale (voir Figure 2 dans Lee et al. 2015). Cette morphologie cohérente est conforme aux prédictions de nos simulations pour un flot poussé par un jet pulsé. Modéliser l'émission de raie en $H_2$ pourrait s'avérer utile afin de discriminer entre notre

<footnote n="1">1. Le rapport entre la pression thermique et la pression magnétique dans les modèles de base de Shang et al. (2020) peut être retrouvé comme valant $\beta = \alpha_b^{-2} (v_A/a)^{-2}$, où $\alpha_b = [0, 0.1, 1]$ est leur paramètre d'échelle et $(v_A/a)$ est le rapport entre la vitesse d'Alfvén et la vitesse du son dans la solution d'équilibre magnétostatique de Li&Shu (1996) (Eq. 69), et qui augmente avec le paramètre d'aplatissement $n$ de l'enveloppe. Une importante diminution de la largeur de la couche par rapport au cas hydrodynamique est observée uniquement dans les simulations avec $n \ge 4$ et $\alpha_b = 1$, correspondant aux valeurs situées dans le plan médian de l'enveloppe avec des valeurs $\beta \le 3.5$ proches de l'équipartition</footnote>
scénario et le modèle de vent grand angle pulsé proposé par Zhang et al. (2019).

En somme, les perspectives d'études futures seraient d'améliorer nos simulations en incluant l'influence de la chimie du gaz (des atomes et des molécules comme H₂ et CO), du champ magnétique, de la rotation et de la précession du jet, tout cela afin de modéliser des observations de flots spécifiques avec davantage de détails. Par exemple, à cet effet, j'ai réalisé de mon côté un projet de recherche personnel durant la thèse, et effectué en parallèle de mon travail principal : https://github.com/mrabananahary/AMRVACH2.git et https://github.com/mrabananahary/VACGRACKLE.git. Il s'agit d'un projet Fortran/C/C++, que j'ai mené sur mon temps libre durant les trois années de thèse, afin d'intégrer avec succès au code MPI-AMRVAC un module permettant de résoudre la chimie dite primordiale (incluant H, H+, He, He+, He²+, e⁻) ainsi que de la chimie de H₂ (impliquant H⁻, H₂, H₂⁺). J'y résous aussi le refroidissement du gaz en dehors de l'équilibre d'ionisation jusqu'ici supposée dans nos simulations, ainsi qu'en dehors de l'équilibre chimique. Pour cela, j'ai utilisé la bibliothèque astrophysique de chimie et de refroidissement Grackle, développée en libre accès (https://github.com/grackle-project/grackle.git. J'ai programmé un module MPI-AMRVAC de sorte à intégrer correctement la chimie dans l'algorithme du code, et à adapter cette intégration au problème et modèle physique de cette thèse. Malheureusement, ce n'est qu'au moment de finir de rédiger cette première version du manuscrit que j'ai pu obtenir des tests positifs de vérification et de comparaison avec des simulations de jets protostellaires déjà validées, établies et publiées sur d'autres codes avec processus chimiques. Je n'ai donc pas eu l'occasion de présenter ici les résultats.

## Bilan personnel

Plus généralement, cette thèse a été pour moi l'opportunité d'enrichir et renforcer mes compétences et connaissances non seulement en Astrophysique, mais aussi et surtout dans plusieurs domaines transverses tournant autour des simulations numériques, telles que le calcul haute performance ou l'intelligence artificielle. J'ai non seulement renforcé ma maîtrise des langages informatics que j'ai eu à utiliser (notamment Fortran, Python et C), mais cette thèse m'a aussi donné l'occasion de me former et d'appliquer les connaissances que j'ai acquises dans le cadre d'un projet de recherche scientifique en équipe. Si le voyage n'était pas non plus de tout repos durant cette thèse, et que j'ai aussi fait face à toutes sortes de difficultés, imprévus et épreuves (dont sur le plan scientifique), la thèse m'a offert une expérience enrichissante et intense qui m'a permis de mûrir et d'affûter mes projets d'avenir et mes ambitions sur le plan professionnel, personnel et scientifique.
<figure><img src="image_76.png" /></figure>
BIBLIOGRAPHIE

Agra-Amboage, V., Cabrit, S., Dougados, C., et al. 2014, A&A, 564, A11

Amdahl, G. M. 1967, in Proceedings of the April 18-20, 1967, Spring Joint Computer Conference, AFIPS '67 (Spring) (New York, NY, USA: Association for Computing Machinery), 483-485

Andre, P., Martin-Pintado, J., Despois, D., & Montmerle, T. 1990, A&A, 236, 180

Andre, P., Ward-Thompson, D., & Barsony, M. 2000, in Protostars and Planets IV, ed. V. Man-
nings, A. P. Boss, & S. S. Russell, 59

Antoniucci, S., Nisini, B., Giannini, T., & Lorenzetti, D. 2008, A&A, 479, 503

Arce, H. G., Borkin, M. A., Goodman, A. A., Pineda, J. E., & Halle, M. W. 2010, ApJ, 715, 1170

Arce, H. G. & Sargent, A. I. 2004, ApJ, 612, 342

Arce, H. G., Shepherd, D., Gueth, F., et al. 2007, in Protostars and Planets V, ed. B. Reipurth, D. Jewitt, & K. Keil, 245

Bachiller, R., Cernicharo, J., Martin-Pintado, J., Tafalla, M., & Lazareff, B. 1990, A&A, 231, 174

Bachiller, R., Martin-Pintado, J., & Planesas, P. 1991, A&A, 251, 639

Bally, J., Devine, D., & Alten, V. 1996, ApJ, 473, 921

Béthune, W., Lesur, G., & Ferreira, J. 2017, A&A, 600, A75

Blandford, R. D. & Payne, D. G. 1982, MNRAS, 199, 883

Bontemps, S., Andre, P., Terebey, S., & Cabrit, S. 1996, A&A, 311, 858

Bosch-Ramon, V. 2022, A&A, 660, A5

Cabrit, S. & Raga, A. 2000, A&A, 354, 667

Cabrit, S., Raga, A., & Gueth, F. 1997, in Herbig-Haro Flows and the Birth of Stars, ed. B. Reipurth & C. Bertout, Vol. 182, 163-180

Castellanos-Ramírez, A., Raga, A. C., & Rodríguez-González, A. 2018, ApJ, 867, 29

Cerqueira, A. H., Lefloch, B., Esquivel, A., et al. 2021, Astronomy and Astrophysics - A&A, 645

Colella, P. & Woodward, P. R. 1984, Journal of Computational Physics, 54, 174

Dalgarno, A. & McCray, R. A. 1972, ARA&A, 10, 375
de Colle, F. & Raga, A. C. 2006, A&A, 449, 1061

Decampli, W. M. 1981, ApJ, 244, 124

Dougados, C., Cabrit, S., Lavalley, C., & Ménard, F. 2000, A&A, 357, L61

Downes, T. P. & Cabrit, S. 2003, A&A, 403, 135

Downes, T. P. & Ray, T. P. 1999, A&A, 345, 977

Dutta, S., Lee, C.-F., Liu, T., et al. 2020, ApJS, 251, 20

Einfeldt, B., Munz, C., Roe, P., & Sjögreen, B. 1991, Journal of Computational Physics, 92, 273

El Mellah, I., Sundqvist, J. O., & Keppens, R. 2017, Monthly Notices of the Royal Astronomical Society, 475, 3240

Erkal, J., Nisini, B., Coffey, D., et al. 2021, ApJ, 919, 23

Falle, S. A. E. G. & Raga, A. C. 1993, MNRAS, 261, 573

Federrath, C., Schrön, M., Banerjee, R., & Klessen, R. S. 2014, ApJ, 790, 128

Ferreira, J. & Deguiran, R. 2013, High Energy Density Physics, 9, 67

Frank, A., Ray, T. P., Cabrit, S., et al. 2014, in Protostars and Planets VI, ed. H. Beuther, R. S. Klessen, C. P. Dullemond, & T. Henning, 451

Gojon, R. & Bogey, C. 2017, International Journal of Aeroacoustics, 16, 603

Guilloteau, S., Bachiller, R., Fuente, A., & Lucas, R. 1992, A&A, 265, L49

Guszejnov, D., Grudić, M. Y., Hopkins, P. F., Offner, S. S. R., & Faucher-Giguère, C.-A. 2021, MNRAS, 502, 3646

Habel, N. M., Megeath, S. T., Booker, J. J., et al. 2021, ApJ, 911, 153

Harten, A., Lax, P. D., & Leer, B. v. 1983, SIAM Review, 25, 35

Hartigan, P., Heathcote, S., Morse, J. A., Reipurth, B., & Bally, J. 2005, AJ, 130, 2197

Hartigan, P., Morse, J. A., & Raymond, J. 1994, ApJ, 436, 125

Hartigan, P. & Wright, A. 2015, ApJ, 811, 12

Hervet, O., Meliani, Z., Zech, A., et al. 2017, A&A, 606, A103

Ireland, L. G., Zanni, C., Matt, S. P., & Pantolmos, G. 2021, ApJ, 906, 4

Keppens, R., Meliani, Z., van der Holst, B., & Casse, F. 2008, A&A, 486, 663

Keppens, R., Meliani, Z., van Marle, A. J., et al. 2012, Journal of Computational Physics, 231, 718

Keppens, R., Teunissen, J., Xia, C., & Porth, O. 2021, Computers & Mathematics with Applications, 81, 316, development and Application of Open-source Software for Problems with Numerical PDEs
Könyves, V., André, P., Men'schikov, A., et al. 2015, A&A, 584, A91

Koren, B. 1993, A robust upwind discretization method for advection, diffusion and source terms, ed. C. Vreugdenhil & B. Koren, Notes on Numerical Fluid Mechanics (Germany: Vieweg), 117-138

Krumholz, M. R. & Federrath, C. 2019, Frontiers in Astronomy and Space Sciences, 6, 7

Kwon, W., Fernández-López, M., Stephens, I. W., & Looney, L. W. 2015, ApJ, 814, 43

Lada, C. J. 1985, ARA&A, 23, 267

Lada, C. J. & Fich, M. 1996, ApJ, 459, 638

Larson, R. B. 1969, MNRAS, 145, 271

Lee, C.-F. 2020, A&A Rev., 28, 1

Lee, C.-F., Hirano, N., Zhang, Q., et al. 2015, ApJ, 805, 186

Lee, C.-F., Ho, P. T. P., Li, Z.-Y., et al. 2017, Nature Astronomy, 1, 0152

Lee, C.-F., Mundy, L. G., Reipurth, B., Ostriker, E. C., & Stone, J. M. 2000, ApJ, 542, 925

Lee, C.-F., Mundy, L. G., Stone, J. M., & Ostriker, E. C. 2002, ApJ, 576, 294

Lee, C.-F., Stone, J. M., Ostriker, E. C., & Mundy, L. G. 2001, ApJ, 557, 429

Lesur, G. R. J. 2021, A&A, 650, A35

LeVeque, R., Steiner, O., Gautschy, A., et al. 1998, Computational Methods for Astrophysical Fluid Flow: Saas-Fee Advanced Course 27. Lecture Notes 1997 Swiss Society for Astrophysics and Astronomy, Saas-Fee Advanced Course (Springer Berlin Heidelberg)

Li, Z.-Y. & Shu, F. H. 1996, ApJ, 472, 211

Liang, L., Johnstone, D., Cabrit, S., & Kristensen, L. E. 2020, ApJ, 900, 15

Liu, T., Lacy, J., Li, P. S., et al. 2017, ApJ, 849, 25

Louvet, F., Dougados, C., Cabrit, S., et al. 2018, A&A, 618, A120

Löhner, R. 1987, Computer Methods in Applied Mechanics and Engineering, 61, 323

MacDonald, J. & Bailey, M. E. 1981, MNRAS, 197, 995

Marcowith, A., Ferrand, G., Grech, M., et al. 2020, Living Reviews in Computational Astrophysics, 6, 1

Margulis, M., Lada, C. J., & Snell, R. L. 1988, ApJ, 333, 316

Masson, C. R. & Chernin, L. M. 1993, ApJ, 414, 230

Matzner, C. D. & McKee, C. F. 1999, ApJ, 526, L109
McCaughrean, M., Zinnecker, H., Andersen, M., Meeus, G., & Lodieu, N. 2002, The Messenger, 109, 28

Meliani, Z. & Hervet, O. 2017, Galaxies, 5

Meliani, Z. & Keppens, R. 2009, ApJ, 705, 1594

Meliani, Z., Keppens, R., Casse, F., & Giannios, D. 2007, MNRAS, 376, 1189

Meliani, Z., Keppens, R., & Giacomazzo, B. 2008, A&A, 491, 321

Millas, D., Keppens, R., & Meliani, Z. 2017, MNRAS, 470, 592

Moens, N., Sundqvist, J. O., El Mellah, I., et al. 2022, A&A, 657, A81

Monceau-Baroux, R., Porth, O., Meliani, Z., & Keppens, R. 2015, A&A, 574, A143

Moraghan, A., Smith, M. D., & Rosen, A. 2006, MNRAS, 371, 1448

Moraghan, A., Smith, M. D., & Rosen, A. 2008, MNRAS, 386, 2091

Motte, F. & André, P. 2001, A&A, 365, 440

Myers, P. C. & Ladd, E. F. 1993, ApJ, 413, L47

Narayanan, G., Snell, R., & Bemis, A. 2012, MNRAS, 425, 2641

Offner, S. S. R., Lee, E. J., Goodman, A. A., & Arce, H. 2011, ApJ, 743, 91

Ogihara, M., Kokubo, E., Suzuki, T. K., & Morbidelli, A. 2018, A&A, 615, A63

Ostriker, E. C., Lee, C.-F., Stone, J. M., & Mundy, L. G. 2001, ApJ, 557, 443

Plunkett, A. L., Arce, H. G., Mardones, D., et al. 2015, Nature, 527, 70

Podio, L., Codella, C., Gueth, F., et al. 2016, A&A, 593, L4

Qiu, K., Zhang, Q., Beuther, H., & Yang, J. 2007, ApJ, 654, 361

Qiu, K., Zhang, Q., Wu, J., & Chen, H.-R. 2009, ApJ, 696, 66

Rabenanahary, M., Cabrit, S., Meliani, Z., & Pineau des Forêts, G. 2022, A&A, 664, A118

Raga, A. & Cabrit, S. 1993, A&A, 278, 267

Raga, A. C., Canto, J., Binette, L., & Calvet, N. 1990, ApJ, 364, 601

Raga, A. C., Rodríguez-González, A., Noriega-Crespo, A., & Esquivel, A. 2012, ApJ, 744, L12

Reipurth, B. & Bally, J. 2001, ARA&A, 39, 403

Reipurth, B., Heathcote, S., Yu, K. C., Bally, J., & Rodríguez, L. F. 2000, ApJ, 534, 317

Richer, J. S., Hills, R. E., & Padman, R. 1992, MNRAS, 254, 525
Richer, J. S., Shepherd, D. S., Cabrit, S., Bachiller, R., & Churchwell, E. 2000, in Protostars and Planets IV, ed. V. Mannings, A. P. Boss, & S. S. Russell, 867

Ridge, N. A. & Moore, T. J. T. 2001, A&A, 378, 495

Robert, C. M. T., Méheut, H., & Ménard, F. 2020, A&A, 641, A128

Rodriguez, L. F., Carral, P., Ho, P. T. P., & Moran, J. M. 1982, ApJ, 260, 635

Roe, P. 1981, Journal of Computational Physics, 43, 357

Roe, P. L. 1986, Annual Review of Fluid Mechanics, 18, 337

Rohde, P. F., Walch, S., Seifried, D., et al. 2019, MNRAS, 483, 2563

Santiago-García, J., Tafalla, M., Johnstone, D., & Bachiller, R. 2009, A&A, 495, 169

Schure, K. M., Kosenko, D., Kaastra, J. S., Keppens, R., & Vink, J. 2009, A&A, 508, 751

Shang, H., Allen, A., Li, Z.-Y., et al. 2006, ApJ, 649, 845

Shang, H., Krasnopolsky, R., Liu, C.-F., & Wang, L.-Y. 2020, ApJ, 905, 116

Shu, F. H. 1977, ApJ, 214, 488

Shu, F. H., Najita, J., Ostriker, E. C., & Shang, H. 1995, ApJ, 455, L155

Shu, F. H., Ruden, S. P., Lada, C. J., & Lizano, S. 1991, ApJ, 370, L31

Smith, M. D. 1986, MNRAS, 223, 57

Su, Y.-N., Zhang, Q., & Lim, J. 2004, ApJ, 604, 258

Suttner, G., Smith, M. D., Yorke, H. W., & Zinnecker, H. 1997, A&A, 318, 595

Suzuki, T. K., Ogihara, M., Morbidelli, A., Crida, A., & Guillot, T. 2016, A&A, 596, A74

Tabone, B., Raga, A., Cabrit, S., & Pineau des Forêts, G. 2018, A&A, 614, A119

Terebey, S., Shu, F. H., & Cassen, P. 1984, ApJ, 286, 529

Toro, E. F. 1997, Riemann Solvers and Numerical Methods for Fluid Dynamics A Practical Introduction

Townsend, R. H. D. 2009, ApJS, 181, 391

Tóth, G. & Odstrčil, D. 1996, Journal of Computational Physics, 128, 82

Ulrich, R. K. 1976, ApJ, 210, 377

van der Holst, B. & Keppens, R. 2007, Journal of Computational Physics, 226, 925

van Kempen, T. A., van Dishoeck, E. F., Güsten, R., et al. 2009, A&A, 501, 633

Velusamy, T. & Langer, W. D. 1998, Nature, 392, 685
Velusamy, T., Langer, W. D., & Thompson, T. 2014, ApJ, 783, 6

Völker, R., Smith, M. D., Suttner, G., & Yorke, H. W. 1999, A&A, 343, 953

Vorobyov, E. I., Elbakyan, V. G., Plunkett, A. L., et al. 2018, A&A, 613, A18

Wagner, A. Y., Bicknell, G. V., & Umemura, M. 2012, ApJ, 757, 136

Wang, L.-Y., Shang, H., Krasnopolsky, R., & Chiang, T.-Y. 2015a, ApJ, 815, 39

Wang, L.-Y., Shang, H., Krasnopolsky, R., & Chiang, T.-Y. 2015b, ApJ, 815, 39

Yee, H. 1994

Yu, K. C., Billawala, Y., & Bally, J. 1999, AJ, 118, 2940

Zapata, L. A., Arce, H. G., Brassfield, E., et al. 2014, MNRAS, 441, 3696

Zhang, Q. & Zheng, X. 1997, ApJ, 474, 719

Zhang, Y., Arce, H. G., Mardones, D., et al. 2016, ApJ, 832, 158

Zhang, Y., Arce, H. G., Mardones, D., et al. 2019, ApJ, 883, 1

Zinnecker, H., McCaughrean, M. J., & Rayner, J. T. 1998, Nature, 394, 862
# APPENDICES
<figure><img src="image_77.png" /></figure>
A. Force de contrebalancement d'un cœur singulier aplati

En coordonnées cylindriques de repère (O, $\vec{e}_R$, $\vec{e}_z$), la force de contrebalancement $\vec{F}_p$ des gradients de pression initiaux s'écrit :

$$ \vec{F}_p = \nabla [p(t=0)] = \vec{e}_R \frac{\partial p}{\partial R} + \vec{e}_z \frac{\partial p}{\partial z} \qquad (VII.1) $$

(les gradients de pression et la force de contrebalancement sont contraints dans un plan à angle azimutal constant : $\frac{\partial p}{\partial \varphi} = 0$)

L'Eq. (V.6) nous permet d'obtenir la pression thermique d'un cœur singulier, dans l'hypothèse d'un gaz parfait :

$$ \begin{aligned} p &= \frac{\rho k_B T}{\mu m_H}, \\ &= \frac{\rho a_0 k_B T}{\mu m_H} \frac{r_w^2 R^2}{(R^2 + z^2)^2} \equiv \alpha \frac{R^2}{(R^2 + z^2)^2}, \end{aligned} \qquad (VII.2) $$

En dérivant cette pression, nous obtenons l'expression analytique à employer pour la composante radiale cylindrique et axiale de $\vec{F}_p$ :

$$ (\vec{F}_p)_R = \frac{\partial p}{\partial R} = - \frac{2\alpha R (R^2 - z^2)}{(R^2 + z^2)^3}, \qquad (VII.3) $$

$$ (\vec{F}_p)_z = \frac{\partial p}{\partial z} = - \frac{4\alpha R^2 z}{(R^2 + z^2)^3}, \qquad (VII.4) $$

$$ \text{où } \alpha = \frac{\rho a_0 k_B T r_w^2}{\mu m_H}, \qquad (VII.5) $$

B. Article publié en premier auteur
# Wide-angle protostellar outflows driven by narrow jets in stratified cores

M. Rabenanahary<fnref n="1" />, S. Cabrit<fnref n="1" />,<fnref n="2" />, Z. Meliani<fnref n="3" />, and G. Pineau des Forêts<fnref n="1" />,<fnref n="4" />

Received 17 January 2022 / Accepted 31 March 2022

## ABSTRACT

Most simulations of outflow feedback on star formation are based on the assumption that outflows are driven by a wide angle “X-wind,” rather than a narrow jet. However, the arguments initially raised against pure jet-driven flows were based on steady ejection in a uniform medium, a notion that is no longer supported based on recent observations. We aim to determine whether a pulsed narrow jet launched in a density-stratified, self-gravitating core could reproduce typical molecular outflow properties, without the help of a wide-angle wind component. We performed axisymmetric hydrodynamic simulations using the MPI-AMRVAC code with optically thin radiative cooling and grid refinement down to 5 au, on timescales up to 10 000 yr. Then we computed the predicted properties for the purposes of a comparison with observational data. First, the jet-driven shell expands much faster and wider through a core with steeply decreasing density than through an uniform core. Second, when blown into the same singular flattened core, a jet-driven shell shows a similar width as a wide-angle wind-driven shell in the first few hundred years, but a decelerating expansion on long timescales. The flow adopts a conical shape, with a sheared velocity field along the shell walls and a base opening angle reaching up to $\alpha \simeq 90^{\circ}$. Third, at realistic ages of ~10 000 yr, a pulsed jet-driven shell shows fitting features along with a qualitative resemblance with recent observations of protostellar outflows with the Atacama Large Millimeter Array, such as HH46-47 and CARMA-7. In particular, similarities can be seen in the shell widths, opening angles, position-velocity diagrams, and mass-velocity distribution, with some showing a closer resemblance than in simulations based on a wide-angle “X-wind” model. Therefore, taking into account a realistic ambient density stratification in addition to millenia-long integration times is equally essential to reliably predict the properties of outflows driven by a pulsed jet and to confront them with the observations.

**Key words.** stars: formation – stars: pre-main sequence – methods: numerical – ISM: jets and outflows – shock waves – hydrodynamics

## 1. Introduction

The most spectacular, and often the first observed, signature of the birth of a new star is the formation of a slow bipolar outflow of molecular gas. This phenomenon starts in the early protostellar phase of stellar mass assembly (Class 0), persists during the envelope dispersion phase (Class 1), and is ubiquitous across all masses (for a review, see e.g., Frank et al. 2014). Given their large sizes and high mass and momentum fluxes, ubiquity, and duration, molecular outflows are believed to play a key role in star formation on both small and large scales: recent numerical simulations (see e.g., Krumholz & Federrath 2019, for a review) suggest that they could be the main feedback agent setting the final stellar mass and core-to-star efficiency (via the removal of parent core material), and regulating the IMF peak, multiplicity fraction, and star formation efficiency at cluster scales (via the disruption of infall streams and replenishment of turbulence).

The exact effect of outflow feedback, however, depends on the assumed structure for the “primary” protostellar wind sweeping up the slow outflow. Two wind models are currently in use, both involving a fast and dense jet along the flow axis (as commonly observed in atomic or molecular tracers, see Frank et al. 2014) but strongly differing in the momentum injected at wider angles.

The first and most frequently used wind prescription in feedback simulations is that of Matzner & McKee (1999, hereafter MM99). It assumes a wide-angle wind radially expanding at a constant speed ($\simeq 100 \text{ km s}^{-1}$) over all angles, with a steep density increase towards the axis responsible for the appearance of an axial “jet.” This asymptotic structure was first derived for an “X-wind” magnetically launched from the disk inner edge (Shu et al. 1995) and only applies to hydromagnetic winds launched radially from a narrow region (MM99). Towards the equator, the wind momentum flux is still a sizeable fraction of the isotropic wind case<fnref n="1" /> and can directly impact equatorial infall.

The second wind model, motivated by more recent MHD simulations and observations, assumes that the fast axial jet is surrounded by a slower disk wind, ejected within a limited solid angle (Federrath et al. 2014; Rohde et al. 2019). Feedback is then dominated by the jet, with much lesser impact on equatorial regions than in the MM99 prescription. Time variability in the form of episodic outbursts was also shown to affect outflow feedback (Rohde et al. 2019).

In principle, realistic MHD simulations of protostellar wind launching should provide the best wind prescription to adopt. However, the simulated wind structure depends on complex

<footnote n="1">¹ Observatoire de Paris, PSL University, Sorbonne Université, CNRS, LERMA, 75014 Paris, France
e-mail: mialy.rabananahary@obspm.fr</footnote>

<footnote n="2">² Université de Grenoble Alpes, CNRS, IPAG, 38000 Grenoble, France</footnote>

<footnote n="3">³ Observatoire de Paris, PSL University, Université de Paris, CNRS, LUTH 5 Place Jules Janssen, 92190 Meudon, France</footnote>

<footnote n="4">⁴ Université Paris-Saclay, CNRS, Institut d'Astrophysique Spatiale, 91405 Orsay, France</footnote>

<footnote n="1">¹ The fraction is $1/\ln(2/\theta_0) \simeq 1/5$ for a “jet” collimation angle $\theta_0 \simeq 0.01 \text{ rad, cf. Eq. (2) in MM99.}$</footnote>
effects that are still the subject of intense research and debate, such as the magnetosphere-disk interaction, the distribution of magnetic flux retained in the disk long after its formation, the turbulent viscosity and resistivity, and non-ideal effects (see e.g., Ireland et al. 2021; Ferreira & Deguiran 2013; Béthune et al. 2017).

An independent approach to determining the most realistic wind model for outflow feedback studies is to simulate the swept-up outflow properties on protostellar core scales $\le 0.1$ pc (where the ambient density structure is dominated by self-gravity and not yet perturbed by cloud inhomogeneities) and see which wind model best reproduces the observed outflow shapes and kinematics. Such a comparison was performed early on for two extremes in wind collimation: the wide-angle X-wind model (Shu et al. 1995; Matzner & McKee 1999) and a pure jet driving the outflow through large bowshocks (Masson & Chernin 1993; Raga & Cabrit 1993). Successes and caveats were identified in each case, based on early observations (see Cabrit et al. 1997; Lee et al. 2001; Arce et al. 2007) and these are briefly summarized and updated below.

First, in models of outflows driven by a wide-angle X-wind, the ambient medium is assumed to have a steep $1/r^2$ density decrease and a moderate degree of magnetic flattening. It is, in addition, assumed to mix instantaneously with shocked wind material. The swept-up shell then expands radially in a self-similar fashion that can reproduce several observed features of molecular outflows: “Hubble-law” kinematics $V \propto z$, mass-velocity distribution with a power-law slope $\gamma \approx -2$ (before opacity correction), and parabolic shapes with a wide base opening angle (Shu et al. 1991; Li & Shu 1996; Lee et al. 2001; Shang et al. 2006, 2020). An intrinsic caveat of this model, however, is that the uniform wind speed over all angles predicts much flatter internal shocks than observed in shock-excited $H_2$ along outflow axes. The observed curved $H_2$ bows requires a sharp drop of wind ram pressure away from the axis (see discussions in Lee et al. 2001; Arce et al. 2007).

Second, simulations of jet-driven bowshocks, in contrast, reproduce the curved morphology of internal shock fronts seen in $H_2$ very well (e.g., Suttner et al. 1997; Völker et al. 1999) as well as the associated “spur-like” features in CO (Lee et al. 2000). They can also reproduce the observed mass-velocity relations (Downes & Cabrit 2003; Moraghan et al. 2008) and apparent “Hubble-laws” when the jet is variable and precessing (Völker et al. 1999; Rohde et al. 2019). The jet-driven model for outflows was strongly criticized, however, for predicting too highly elongated cavities on long timescales (Ostriker et al. 2001), too much overlapping blueshifted and redshifted emission over a wide range of inclinations (Lee et al. 2001) as well as overly low velocities $\approx 0.03 \text{ km s}^{-1}$ when the bowshock had expanded to typical outflow widths of 10 000 au (Arce et al. 2007). As a result, it is commonly believed that jets alone cannot explain outflows with a wide opening angle, as reported for several evolved Class 1 protostars.

In order to combine the strengths of each model, a “dual wind” structure has been invoked with both a fast jet and a slower wide angle wind, where the latter would increasingly dominate at later times (Yu et al. 1999; Velusamy & Langer 1998; Arce et al. 2007; Zapata et al. 2014). There are several good reasons, however, to reconsider pure jet-driven shells as the potential origin for molecular outflows.

First, the critiques of Ostriker et al. (2001) and Lee et al. (2001) were based on models of jet-bowshocks in a uniform, or quasi-uniform, ambient medium (with at most a factor of 2 in density variation over the computational domain). In contrast,

a steep radial $1/r^2$ density decrease is assumed in wide-angle wind-driven models to yield the apparent “Hubble-law” acceleration (Shu et al. 1991; Lee et al. 2001). Such a decrease is expected on protostellar core scales $\le 0.1$ pc as a result of self-gravity. Raga & Cabrit (1993) and Cabrit et al. (1997) showed that it could produce a wider opening angle for jet-driven shells, more similar to observed flows, a result confirmed for steady jets on $<1000$ yr timescales by Moraghan et al. (2008). Now that a steep density stratification is widely confirmed by observations of protostellar cores (e.g., van der Marel et al. 2009, for HH46-47) and that computational capabilities have greatly improved, it is important to explore the predicted effect on jet-driven shells over a broader parameter space and longer timescales than was feasible in the early study of Moraghan et al. (2008).

Second, the issue raised by Arce et al. (2007) related to insufficient bowshock speed at large widths no longer applies with a pulsed jet. New internal bowshocks generated by the jet variability will replenish a slow jet-driven shell with faster material at observable speeds (Raga & Cabrit 1993; Völker et al. 1999). In addition, the interaction between successive bowshocks will decrease their transverse speed, possibly alleviating the excessive blue and red overlap predicted at early times by Lee et al. (2001). Recent Atacama Large Millimeter Array (ALMA) observations revealed multiple $H_2$ bowshocks along outflow axes, connected to nested CO cavities along their flanks (e.g., HH212 and HH46-47 in Lee et al. 2015; Zhang et al. 2016). While some have been modeled with wind-driven shells (Zhang et al. 2019), it is important to have similar predictions for nested jet bowshocks at realistic ages $\approx 10^4$ yr.

Third, recent observational studies show that CO outflows are more elongated and collimated than initially thought: cloud-wide CO maps, as well as optical and infrared imaging surveys, show that at least 40% of outflows are more than a parsec long (Frank et al. 2014, and references therein). A striking example is the B5-IRS1 outflow, driven by a Class 1 protostar. While it exhibits a (projected) full opening angle $\alpha \approx 100^\circ$ at its base, argued as evidence for a wide-angle wind, unbiased CO maps of its parent cloud reveal that each lobe extends (at least) up to 2.2 pc from the source (Arce et al. 2010). The length-to-width ratio is then $q > 11$ (see maps in Frank et al. 2014). This is inconsistent with current models of shells driven by a wide angle X-wind (which predict an aspect ratio $q < 3.5$ for a base opening angle $\alpha \ge 70^\circ$, see Table 2 in Shang et al. 2020). The parsec size of many outflows appears more suggestive of jet-driven flows.

Estimated outflow opening angles are also affected by several biases: angular resolution, the height at which they are measured (Velusamy & Langer 1998; Velusamy et al. 2014), and inclination (angles appear wider in flows seen closer to pole-on). Two recent studies in Orion minimize these biases by providing uniform measurements at the same (high) linear resolution and projected height over randomly selected samples. In an ALMA survey of the widths of 22 (mostly Class 0) CO outflows (Dutta et al. 2020), 50% subtend projected full-opening angles in the range $\alpha = 25^\circ-65^\circ$ at a projected altitude $z_{\text{proj}} = 800$ au. In a sample of 29 older outflow cavities (mostly Class 1) imaged in scattered light with HST (Habel et al. 2021), we see that 50% are in the range $\alpha = 8^\circ-46^\circ$ at $z_{\text{proj}} = 8000$ au, and the fraction of point sources (viewed down the cavity interior) suggests a maximum deprojected opening angle $\alpha_{\text{deproj}} \le 70^\circ$. Therefore, CO outflows seem to be more collimated on average than previously believed and it is necessary to investigate whether a pure jet (in a stratified core) could reproduce typical observed widths, before drawing any conclusions on a dominant contribution from a wider angle wind.
Here, we examine this issue by presenting the first high-resolution simulations of pure jet-driven shells in strongly stratified cores, up to ages of 10 000 yr and physical scales of 0.1 pc. For the first time in jet simulations, we consider the same flattened core structure as in the competing model of outflow driven by a wide-angle X-wind (Li & Shu 1996; Lee et al. 2001). We show that pulsed conical jets propagating in this environment sweep up a wider outflow cavity than in a uniform medium, with a width and opening angle that are compatible with recent outflow surveys. The predicted position-velocity diagrams and mass-velocity relation also show a promising qualitative agreement with recent ALMA outflow observations at high resolution and sensitivity, without any of the caveats noted previously for jet bowshocks in uniform media.

In Sect. 2, we present our numerical method and generic set up. In Sect. 3, we present the effect of a 1D and 2D density stratification on a jet-driven shell. In Sect. 4, we introduce a small jet opening angle and explore the effect of various free parameters on the cavity shape and kinematics. In Sect. 5, we present a simulation up to 10 000 yr and compute predicted flow widths, position-velocity diagrams, and mass-velocity relation, finding excellent qualitative agreement with recent ALMA observations. Section 6 summarizes our main results and conclusions.

## 2. General numerical setup

### 2.1. Governing equations, code, and numerical method

We performed axisymmetric 2D hydrodynamic simulations in cylindrical coordinates (R, z), using the Message Passing Interface-Adaptive Mesh Refinement Versatile Advection Code (MPI-AMRVAC; Keppens et al. 2021). The hydrodynamics module of this finite volume, cell-centered code solves the hydrodynamic equations of mass, momentum, and energy conservation described, respectively, by:

$$ \frac{\partial \rho}{\partial t} + \nabla \cdot (\rho \mathbf{v}) = 0, \qquad (1) $$

$$ \frac{\partial(\rho\mathbf{v})}{\partial t} + \nabla \cdot (\rho\mathbf{v}\mathbf{v}) + \nabla p = \mathbf{F}_p, \qquad (2) $$

$$ \frac{\partial e}{\partial t} + \nabla \cdot (e\mathbf{v} + p\mathbf{v}) = -n_H^2 \Lambda(T), \qquad (3) $$

where $\rho$ is the mass density, $\mathbf{v}$ is the velocity vector, $p$ is the thermal pressure, and $e = p/(\gamma - 1) + \rho \mathbf{v}^2/2$ is the total (thermal and kinetic) energy density, with $\gamma$ the adiabatic index (taken here as 5/3). Two source terms are introduced on the right-hand side: following Lee et al. (2001), an inward-directed force field, $\mathbf{F}_p = \nabla p(t = 0)$, is imposed to maintain the unperturbed stratified ambient core in hydrostatic equilibrium at any time. In addition, optically thin equilibrium radiative cooling is included as a source term $-n_H^2 \Lambda(T)$ in the energy equation, with $n_H$ the number density of hydrogen nuclei and $T$ the gas kinetic temperature (van Marle & Keppens 2011). We consider an atomic gas with a standard helium fraction $x(He) = n(He)/n_H = 0.1$, such that $n_H = \rho/(1.4m_H)$ with $m_H$ the mass of a proton. The temperature is inferred from gas pressure using the perfect gas law, $T = p/(n_{tot}k_B)$, with $n_{tot}$ as the total number of particles per unit volume. Given the moderate shock speeds encountered in our simulation, we assume that hydrogen and helium remain mostly neutral, so that $n_{tot} \approx n(H) + n(He) = 1.1n_H$. The cooling curve $\Lambda(T)$ depends on the local temperature, $T$, and the gas metallicity. Here, we use an atomic cooling curve with solar metallicity. In Sect. 3, we use for comparison purposes the

same cooling function as Lee et al. (2001), namely the curve of MacDonald & Bailey (1981) at temperatures above $10^4$ K and Dalgarno & McCray (1972) below $10^4$ K (with an ionization fraction $x_{ion} = 10^{-3}$). In the rest of the paper (Sects. 4 and 5) we keep the same low-temperature cooling curve but we adopt above $10^4$ K the updated cooling function $\Lambda_{HD}(T)$ from Schure et al. (2009) generated with the SPEX code (Kaastra & Mewe 2000), which takes into account more up-to-date atomic parameters than MacDonald & Bailey (1981). In each simulation, the minimum temperature for radiative cooling is set at the initial ambient core temperature. This prevents the non-realistic radiative cooling of unshocked material inside the jet beam and the ambient core.

To solve Eqs. (1)–(3), we chose a Harten–Lax–van Leer contact (HLLC) scheme (Li 2005) with a minmod limiter. This combination is extremely robust in handling the shocks encountered in our problem.

As for the boundary conditions, the jet axis ($R = 0$) is treated as an axisymmetric boundary; the boundary conditions at $z = 0$ are fixed inside the jet inlet ($R < R_j$) and open in the surrounding ambient core. The outer limits of the computational domain are treated as open boundaries. The jet inlet (inside which the density and velocity vectors are prescribed at each time step) is a conical domain of radius $R = R_j$ at $z = 0$, height $z_i$, and semi-open angle $\theta_j$.

Due to the high dynamic ratio between inner jet shock scale and propagation scales studied in this paper, it is crucial to use adaptive mesh refinement (AMR) to resolve the shocks in a cost-effective manner. Apart from Sect. 3 (where we use the same fixed grid resolution as Lee et al. 2001 for comparison purposes), all simulations performed in this paper have a base grid that is allowed to be locally refined up to four times, doubling the resolution at each new level of refinement (i.e., with a highest grid resolution at level 5 that is $2^4$ finer than at level 1). The refinement-derefinement criterion is based on a Lohner error estimator on the quantity $n_H^2 \Lambda(T)$. We further ensure that the highest level of refinement is always present inside the jet beam.

### 2.2. Velocity variability of the jet

Pulsed jets are simulated by time-varying the jet velocity around its initial value $v_0$, as

$$ v_j(t) = \|\mathbf{v}_j(t)\| = v_0 h(t), \qquad (4) $$

where $h(t)$ is the dimensionless variability profile defining how the jet will pulse.

Following Lee et al. (2001), we keep the jet mass-flux constant over time in the present paper, by imposing an inverse variation of the jet density:

$$ \rho_{j0}(t) = \frac{\rho_{j0}}{h(t)} = \frac{\rho_{j0}v_0}{v_j(t)}, \qquad (5) $$

where $\rho_{j0}$ is the initial jet density at $z = 0$, $R = R_j$, and $t = 0$.

### 2.3. Tracers

Following Porth et al. (2014), the code also solves the advection equation of two passive fluid tracers, tracer_a and tracer_j, which trace fluid parcels originating from the ambient core and from the jet, respectively. They evolve according to:

$$ \frac{\partial \text{tracer}_k}{\partial t} + \mathbf{v} \cdot \nabla \text{tracer}_k = 0, \text{ for } k = a, j. \qquad (6) $$
<table><caption>**Table 1.** Parameters of the pulsed cylindrical jet simulations (Sect. 3).</caption><thead><tr><th colspan="5">Fixed parameters (from Lee et al. 2001)</th></tr></thead><tbody><tr><td>One-sided jet mass-loss rate</td><td colspan="4">M&#x0307; = 6.0 &times; 10<sup>-8</sup> M<sub>&odot;</sub> yr<sup>-1</sup></td></tr><tr><td>Jet semi-opening angle</td><td colspan="4">&theta;<sub>j</sub> = 0&deg;</td></tr><tr><td>Jet radius</td><td colspan="4">R<sub>j</sub> = 2.5 &times; 10<sup>15</sup> cm</td></tr><tr><td>Mean jet velocity</td><td colspan="4">v<sub>0</sub> = 120 km s<sup>-1</sup></td></tr><tr><td>Jet velocity variation</td><td colspan="4">v<sub>j</sub>(t) = v<sub>0</sub> + &Delta;V sin (2&pi;t/P)</td></tr><tr><td>Semi-amplitude</td><td colspan="4">&Delta;V = 60 km s<sup>-1</sup></td></tr><tr><td>Jet density variation</td><td colspan="4">&rho;<sub>j</sub>(t) = &rho;<sub>j0</sub> [v<sub>0</sub>/v<sub>j</sub>(t)] (constant mass-flux)</td></tr><tr><td>Jet density at t = 0</td><td colspan="4">&rho;<sub>j0</sub> = 1.6 &times; 10<sup>-20</sup> g cm<sup>-3</sup></td></tr><tr><td>Jet temperature</td><td colspan="4">T<sub>j</sub> = 270 K</td></tr><tr><td>Ambient core temperature</td><td colspan="4">T<sub>a</sub> = 30 K</td></tr><tr><td>Radiative cooling function</td><td colspan="4">&Lambda;(270 K &le; T &lt; 10<sup>4</sup> K) from Dalgarno &amp; McCray (1972)<br/>&Lambda;(T &ge; 10<sup>4</sup> K) from MacDonald &amp; Bailey (1981)</td></tr><tr><td>Simulation domain</td><td colspan="4">(R, z) = (7.0 &times; 10<sup>16</sup> cm, 1.4 &times; 10<sup>17</sup> cm) = (4679 au, 9358 au)</td></tr><tr><td>Number of cells</td><td colspan="4">n<sub>R</sub> &times; n<sub>z</sub> = 336 &times; 672</td></tr><tr><td>Resolution</td><td colspan="4">&Delta;R = &Delta;z = 2.08 &times; 10<sup>14</sup> cm = 14 au</td></tr><tr><th colspan="5">Variable parameters</th></tr><tr><td>Model name</td><td>PCJ-U <sup>(a)</sup><br/>(Fig. 1a)</td><td>PCJ-Z<br/>(Fig. 1b)</td><td>PCJ-R<br/>(Fig. 1c)</td><td>PCJ-RW <sup>(b)</sup><br/>(Fig. 2)</td></tr><tr><td>Core density profile</td><td>Uniform</td><td>z-stratified</td><td colspan="2">Flattened singular core</td></tr><tr><td>&rho;<sub>a</sub>(r)</td><td>&rho;<sub>a0</sub></td><td>&rho;<sub>a0</sub>(1 + z/z<sub>c</sub>)<sup>-2</sup><br/>z<sub>c</sub> = 10<sup>16</sup> cm</td><td colspan="2">&rho;<sub>a0</sub> sin<sup>2</sup> &theta; r<sub>0</sub><sup>2</sup>(R<sup>2</sup> + z<sup>2</sup>)<sup>-1</sup><br/>r<sub>0</sub> = 2.5 &times; 10<sup>15</sup> cm</td></tr><tr><td>Core density at z = 0, R = 2.5 &times; 10<sup>15</sup> cm</td><td colspan="2">&rho;<sub>a0</sub> = 1.6 &times; 10<sup>-20</sup> g cm<sup>-3</sup></td><td colspan="2">1.6 &times; 10<sup>-18</sup> g cm<sup>-3</sup></td></tr><tr><td>Core/jet density contrast at z = 0</td><td colspan="2">&eta; = &rho;<sub>j0</sub>/&rho;<sub>a0</sub> = 1</td><td colspan="2">0.01</td></tr><tr><td>Jet variability period</td><td colspan="2">P = 310 yr</td><td colspan="2">115 yr</td></tr><tr><td>Simulation age</td><td colspan="2">t = 610 yr</td><td colspan="2">296 yr</td></tr></tbody></table>

Notes. (a)Identical parameters to pulsed jet simulation of Lee et al. (2001) (see their Fig. 7). (b)Identical parameters to pulsed wind simulation of Lee et al. (2001) (cf. their Fig. 12), except that the wide-angle wind is replaced here by a cylindrical jet of same mass-flux, injection radius, and velocity variability.

These tracers are dimensionless numerical quantities. To improve numerical precision and accurately track mixing between the jet and ambient core material, they are taken to cover a wide range [0, 10<sup>7</sup>]. Thus, at t = 0, tracer<sub>j</sub> (resp. tracer<sub>a</sub>) is initialized to 10<sup>7</sup> inside the jet beam (resp. the surrounding core) and is set to zero elsewhere.

From those tracers, we may compute the local fraction of material originating from the ambient core and from the jet, f<sub>a</sub> and f<sub>j</sub>, inside each mesh cell and at any time step, as

$$f_k = \frac{\text{tracer}_k}{\text{tracer}_a + \text{tracer}_j}, \quad \text{for } k = a, j. \qquad (7)$$

### 3. Pulsed cylindrical jet in a non-uniform medium

In this section, we investigate how the shape of the shell driven by the same cylindrical jet as in Lee et al. (2001) is influenced by a more realistic, steep density decrease in the ambient core (Sect. 3.2). We then compare with the wide-angle wind-driven shell modeled in Lee et al. (2001) for the same age and core stratification (Sect. 3.3).

#### 3.1. Setup and choice of density stratifications

To make a comparison with Lee et al. (2001), we adopted their pulsed cylindrical jet model (thereafter PCJ) with the same jet

radius, R<sub>j</sub>, temperature, T<sub>j</sub>, sinusoidal variation of jet velocity, v<sub>j</sub>(t), inversely varying jet density, ρ<sub>j</sub>(t), and constant one-sided jet mass-loss rate, M&#x0307;. We also adopted the same uniform spatial resolution of 14 au, radiative cooling law (using MacDonald & Bailey 1981 above 10<sup>4</sup> K), domain size, ambient temperature, and simulation ages. The particular setup presented in this section is summarized in Table 1.

We then investigated three different ambient core stratifications: first, just as in the jet-driven simulations of Lee et al. (2001), we considered a uniform core, with

$$\rho_a(R, z) = \rho_{a0}. \qquad (8)$$

Second, we investigated a z-stratified core following the prescription in Cabrit et al. (1997):

$$\rho_a(R, z) = \rho_{a0} \left(1 + \frac{z}{z_c}\right)^{-2} \quad \text{with} \quad z_c = 10^{16} \text{ cm}, \qquad (9)$$

where the value of z<sub>c</sub> is justified by observations of protostellar cores showing a steep stratification from 20 000 au down to (at least) 1000 au (see e.g., Motte & André 2001; van der Marel et al. 2009). Lee et al. (2001) investigated (see their Sect. 3.6) a somewhat similar z-stratification ρ<sub>a</sub>(z) ∝ 1/[1 + (z/z<sub>c</sub>)<sup>2</sup>] but with a flattening scale z<sub>c</sub> = 1.25 × 10<sup>17</sup> cm that is ten times greater than in Eq. (9) and equal to the size of their computational box. The ambient density thus dropped only by a factor of 2 across
<figure><img src="image_78.png" /><figcaption>**Fig. 1.** Maps of hydrogen nuclei density $n_H$ from simulations of the same cylindrical jet as in Lee et al. (2001) (see Table 1). The three columns confront three different ambient core stratifications, and the two rows show the map at two ages $t = 296$ yr (top) and 610 yr (bottom). The jet is launched in Figs. 1a and d through a uniform ambient core with profile $\rho(R, z) = \rho_{a0}$; in Figs. 1b and e through an $z$-stratified ambient core $\rho(R, z) = \rho_{a0}/(1+z/z_c)^2$ with $z_c = 10^{16}$cm; and in Figs. 1c and f through a flattened singular core $\rho(r, \theta) = \rho_{a0} \sin^2 \theta(r_0/r)^2$, where $r$ is the spherical radius and $r_0 = 2.5 \times 10^{15}$ cm. All core density profiles have the same base density at $R = R_j$, $\rho_{a0} = 1.6 \times 10^{-20}$ g cm$^{-3}$, density-matched with the jet at $t = 0$. W note how the jet-driven shell expands faster and wider through an increasingly stratified core, whereas the nested shells also grow wider.</figcaption></figure>

the whole domain and no strong effect on shell morphology was reported.

Third, we considered the same flattened singular core profile as in the wide angle wind-driven simulations of Lee et al. (2001), namely:

$$ \rho_a(r, \theta) = \rho_{a0} \sin^2 \theta \left(\frac{r_0^2}{r^2}\right) \quad \text{with} \quad r_0 = 2.5 \times 10^{15} \text{cm}, \quad (10) $$

where $r = (R^2 + z^2)^{1/2}$ denotes the spherical radius and $\theta$ is the polar angle from the jet axis. This $1/r^2$ decrease is appropriate for a self-similar singular isothermal core. The $\sin^2 \theta$ dependence is an approximation for the flattened magnetostatic equilibrium solution of Li & Shu (1996) with $n = 2$ (where this parameter measures the degree of magnetic support and equatorial flattening). When a core is stratified according to this particular solution and is then swept-up by an X-wind, it can reproduce the typical collimation, parabolic shape, and mass-velocity distribution in CO outflows, while also being consistent with the mean observed flattening of prestellar cores (Li & Shu 1996). Hence, analytical models and numerical simulations of X-wind driven shells (Lee et al. 2000, 2001; Zhang et al. 2019) widely adopt the expression in Eq. (10) as a “standard” ambient medium.

However, to our knowledge, it was never used in jet-driven outflow simulations until now.

## 3.2. Effect of ambient core stratification on jet-driven shells

Figure 1 shows the results of our simulations at the same ages and for the same cylindrical jet propagating into three different density distributions: uniform (model PCJ-U, Figs. 1a and d), $z$-stratified (model PCJ-Z, Figs. 1b and e), and flattened singular core (model PCJ-R, Figs. 1c and f). All three have the same value of $\rho_{a0} = 1.6 \times 10^{-20}$ g cm$^{-3}$ in Eqs. (8), (9), and (10), respectively. Since $R_j = 2.5 \times 10^{15}$ cm = $r_0$, the ambient density at the jet base ($z = 0, R = R_j$) is also the same in all three cases (equal to $\rho_{a0}$), as well as the initial ($t = 0$) jet-ambient density contrast $\eta = 1$ at this point. All model parameters are summarized in Table 1.

Going from left to right in Fig. 1, the ram pressure constraints exerted on both the main and nested shells are relaxed along the $z$-axis, and then along both the radial and polar directions. As we may see, this change leads to a main shell expanding faster along $z$ and wider over time.

As already noted by Lee et al. (2001), we confirm that a uniform medium produces a narrow and roughly cylindrical jet-driven shell, unlike observed outflows. However, introducing a
steep $1/z^2$ stratification here leads to a more conical leading shell as time proceeds, confirming the analytical and numerical predictions of Raga & Cabrit (1993) and Cabrit et al. (1997). For example, at 610 yr, the shell full width at the top of the domain is increased by a factor 4 between PCJ-U and PCJ-Z (from $1.3 \times 10^{16}$ cm to $5.6 \times 10^{16}$ cm). In the flattened singular core PCJ-R, the rarefied polar holes create an even wider shell with a parabolic shape and a twice larger full width than in PCJ-Z ($1.0 \times 10^{17}$ cm).

Even though the ambient density at the jet inlet ($z = 0, R = R_j$) was kept the same, a steeper density stratification also increases the shell full width near its base. From Fig. 1, we can measure full widths at $z = 800$ au ($= 1.2 \times 10^{16}$ cm) of $W_{800} = 0.95 \times 10^{16}$ cm, $1.16 \times 10^{16}$ cm and $1.75 \times 10^{16}$ cm for the uniform medium, $1/z^2$ decrease and flattened singular core, respectively.

In addition to the leading shell carved by the supersonic jet head, jet variability produces successive and periodic internal working surfaces (thereafter IWS) where high-pressure shocked material is ejected sideways, forming bowshocks expanding inside the leading shell and producing nested “secondary” shells, visible in the $t = 610$ yr snapshot in Fig. 1. As the leading shell expands faster and wider in a stratified core, its inner density distribution (and, hence, the pressure locally exerted on the nested shells) drops more steeply than through an uniform core. This allows the nested shells to expand more widely as well. In parallel, each IWS still moves along the jet beam at the same velocity, independently of the core stratification. This is because the IWS propagation speed only depends on the jet velocity and density conditions upstream and downstream of the forming working surface (Raga et al. 1990). These conditions are entirely determined by the jet variability properties, which remain unmodified for each of our core density profiles.

### 3.3. Comparison with a wide-angle wind-driven shell

Here, we adopt the same setup as in PCJ-R (Figs. 1c, f), but with a flattened singular core that is 100 times denser and identical to that considered in the wide-angle wind models of Lee et al. (2001), with $\rho_{a,0} = 1.6 \times 10^{-18}$ g cm$^{-3}$. We also adopt a shorter variability period $P = 115$ yr than in Sect. 3.2, so that our pulsed cylindrical jet has the exact same mass-loss rate, injection radius, and velocity variability as their pulsed wide-angle wind. The parameters of this new model (PCJ-RW) are summarized in the last column of Table 1.

Figure 2 shows the resulting density map of the PCJ-RW simulation at an age $t = 296$ yr. This map can directly be compared with the wide-angle wind simulation at the same age in Fig. 12 of Lee et al. (2001). At this early age, the jet-driven shell opens as wide as for the pulsed wide-angle wind in the same ambient stratification; namely, we measure similar maximum shell widths of $2.6 \times 10^{16}$ cm (~$10 R_j$) for both models, and a similar full opening angle<fnref n="2" /> at $z = 800$ au of $\alpha_{800} = 70^\circ$ for PCJ-RW and $\alpha_{800} = 80^\circ$ for the wide-angle wind model.

### 3.4. Summary

The usual criticism of jet-driven shells producing overly narrow opening angles appears no longer valid when a realistic stratified ambient medium is considered. The morphology of the shell

<figure><img src="image_79.png" /><figcaption>Fig. 2. Density snapshot at $t = 296$ yr of pulsed cylindrical jet model PCJ-RW, with same injected mass-loss rate, velocity variability, and ambient density distribution as the pulsed wide-angle wind model in Fig. 12 of Lee et al. (2001). The resulting shell size and opening angle are the same as for the wind-driven shell, at this young age.</figcaption></figure>

driven by a pulsed jet is strongly affected by a steep stratification in density of the ambient core surrounding the jet. Wider shells are formed with increased opening angle near the base of the outflow and conical or parabolic shapes on large scales are highly reminiscent of observed CO outflows. Furthermore, in a standard flattened singular core, the jet-driven shell initially expands as wide as with a wide-angle wind. In the rest of this paper, we examine which factors affect the jet-driven shell shape (see Sect. 4) and how it evolves on longer timescales ~10 000 yr, which more closely resemble the actual outflow ages (see Sect. 5).

## 4. Shells driven by a conical jet in a flattened singular core

In this section, we investigate the influence that several physical parameters (expected to vary among observed sources) have on the morphology and kinematics of jet-driven shells. For consistency, we keep the singular flattened core stratification in Eq. (10) (as in wide-angle wind-driven models), but we adopt a more realistic conical jet geometry. High-resolution jet observations suggest jet half-opening angles $\theta_j$ of a few degrees on the scales of our simulations, for example $\theta_j = 2^\circ$ out to 800 au in the atomic jet of RW Aur (Dougados et al. 2000) and $\theta_j \approx 5^\circ$ from 1000–12 000 au in the CO jet of IRAS04166+2706 (Santiago-García et al. 2009). We explore a similar range in our simulations.

### 4.1. Setup: Jet spray nozzle, radiative cooling, and AMR

Similar to Völker et al. (1999), we introduce a spray angle in the jet inlet by taking a velocity vector that is radially diverging from a virtual point:

$$ (R, z) = (0, -z_0), \quad \text{with } z_0 = R_j / \tan \theta_j, \qquad (11) $$

<footnote n="2">² Defined following Dutta et al. (2020) as $\alpha_{800} = 2 \arctan(W_{800}/[2 \times 800 \text{ au}])$.</footnote>
<table><caption>**Table 2.** Parameters of pulsed conical high-density jet simulations with the resulting opening angles and full widths (Fig. 3).</caption><thead><tr><th colspan="7">Fixed parameters</th></tr></thead><tbody><tr><td>Mean jet velocity</td><td colspan="6">v<sub>0</sub> = 120 km s<sup>-1</sup></td></tr><tr><td>Jet density variation</td><td colspan="6">ρ<sub>j</sub>(t) = ρ<sub>j0</sub> [v<sub>0</sub>/v<sub>j</sub>(t)] × (R<sub>j</sub><sup>2</sup> + z<sub>0</sub><sup>2</sup>)(R<sup>2</sup> + [z+z<sub>0</sub>]<sup>2</sup>)<sup>-1</sup>, with z<sub>0</sub> = R<sub>j</sub>/ tan θ<sub>j</sub> (constant mass-flux)</td></tr><tr><td>Core density profile</td><td colspan="6">Flattened singular core ρ<sub>a</sub>(r) = ρ<sub>a0</sub> sin<sup>2</sup> θ r<sub>0</sub><sup>2</sup> (r<sup>2</sup>)<sup>-1</sup>, with r<sub>0</sub> = 2.5 × 10<sup>15</sup> cm</td></tr><tr><td>Jet temperature</td><td colspan="6">T<sub>j</sub> = 100 K</td></tr><tr><td>Ambient core temperature</td><td colspan="6">T<sub>a</sub> = 100 K</td></tr><tr><td>Radiative cooling function</td><td colspan="6">Λ(100 K ≤ T < 10<sup>4</sup> K) from Dalgarno & McCray (1972)<br>Λ(T ≥ 10<sup>4</sup> K) from Schure et al. (2009)</td></tr><tr><td>Simulation domain</td><td colspan="6">(R, z) = (7.0 × 10<sup>16</sup> cm, 1.9 × 10<sup>17</sup> cm) = (4679 au, 12 700 au)</td></tr><tr><td>Number of cells</td><td colspan="6">n<sub>R</sub> × n<sub>z</sub> = 56 × 152 for the full grid at AMR level 1</td></tr><tr><td>Maximum resolution</td><td colspan="6">7.8 × 10<sup>13</sup> cm = 5.2 au at AMR level 5</td></tr><tr><td>Snapshot age</td><td colspan="6">700 yr</td></tr><tr><th>Parameter</th><th>Reference model<br>H_REF</th><th>Modified<br>parameter<sup>(a)</sup></th><th>Modified<br>Model name</th><th>α<sub>800</sub><sup>(b)</sup><br>(°)</th><th>W<sub>800</sub><sup>(c)</sup><br>(10<sup>16</sup> cm)</th><th>W<sub>12 700</sub><sup>(d)</sup><br>(10<sup>16</sup> cm)</th></tr><tr><td>Core base density ρ<sub>a0</sub></td><td>1.6 × 10<sup>-18</sup> g cm<sup>-3</sup></td><td>1.6 × 10<sup>-20</sup> g cm<sup>-3</sup></td><td>H DenSA</td><td>112</td><td>3.5</td><td>13.8</td></tr><tr><td>Jet semi opening angle θ<sub>j</sub></td><td>3°</td><td>7°</td><td>H THETA</td><td>95</td><td>2.6</td><td>10.5</td></tr><tr><td>Jet base initial density<sup>(e)</sup> ρ<sub>j0</sub></td><td>1.8 × 10<sup>-17</sup> g cm<sup>-3</sup></td><td>1.8 × 10<sup>-15</sup> g cm<sup>-3</sup></td><td>H DEnSJ</td><td>86</td><td>2.2</td><td>8.4</td></tr><tr><td>Semi-amplitude ΔV</td><td>60 km s<sup>-1</sup></td><td>90 km s<sup>-1</sup></td><td>H VARAMP</td><td>91</td><td>2.4</td><td>8.0</td></tr><tr><td>Jet variability period P</td><td>115 yr</td><td>300 yr</td><td>H PER</td><td>90</td><td>2.4</td><td>7.4</td></tr><tr><td>Jet variability profile h(t)</td><td>1 + ΔV/v<sub>0</sub> sin 2πt/P</td><td>1 + ΔV/v<sub>0</sub> {1 - 2 · mod(t/P, 1)}</td><td>H SAWT</td><td>86</td><td>2.2</td><td>7.2</td></tr><tr><td>Jet radius R<sub>j</sub></td><td>7.5 × 10<sup>14</sup> cm</td><td>3.0 × 10<sup>14</sup> cm</td><td>H RAD</td><td>84</td><td>2.2</td><td>6.9</td></tr></tbody></table>

**Notes.** <sup>(a)</sup>In each modified model, only one parameter at a time is changed with respect to the reference model. <sup>(b)</sup>Shell full opening angle at z = 800 au. We obtain α<sub>800</sub> = 88° for the reference model. <sup>(c)</sup>Full shell width at z = 800 au. We obtain W<sub>800</sub> = 2.3 × 10<sup>16</sup> cm for the reference model. <sup>(d)</sup>Full shell width at z = 12 700 au (top of the grid). We obtain W<sub>12 700</sub> = 7.6 × 10<sup>16</sup> cm for the reference model. <sup>(e)</sup>High jet density case, yielding a one-sided mass-flux M̄ = 6 × 10<sup>-6</sup> M<sub>☉</sub>/yr for the reference values of θ<sub>j</sub> and R<sub>j</sub>.

where R<sub>j</sub> is the jet radius at z = 0 and θ<sub>j</sub> the jet semi-opening angle.

In order to conserve mass-flux throughout the jet inlet, we set a jet density profile ρ<sub>j</sub>(t) that decreases as the inverse square distance to the above-mentioned virtual point and varies inversely in time with the velocity modulus, as

$$ \rho_j(t) = \rho_{j0} \left( \frac{R_j^2 + z_0^2}{R^2 + (z+z_0)^2} \right) \frac{v_0}{v_j(t)}, \quad (12) $$

where ρ<sub>j0</sub> is the jet “base density” at z = 0, R = R<sub>j</sub>, and t = 0. The (constant) jet mass-flux is then given by

$$ \dot{M} = 2\pi (1 - \cos \theta_j) (R_j^2 + z_0^2) \rho_{j0} v_0, \quad (13) $$

$$ = 2\pi (1 - \cos \theta_j) (R_j / \sin \theta_j)^2 \rho_{j0} v_0. \quad (14) $$

Hereafter, and in Sect. 5, we use the updated radiative cooling function from Schure et al. (2009) (see Sect. 2.1). Furthermore, the resolution and accuracy are improved by enabling AMR up to level 5, with a 84 au minimal and 5.2 au maximal spatial resolution. The computational grid is also expanded to 1.9 × 10<sup>17</sup> cm in z.

## 4.2. Model parameters

First, we defined a set of typical parameters for the reference models. For consistency, we kept the same ejection velocity law and ambient density distribution (singular flattened core) as in our model PCJ-RW in Table 1. We adopted a reference value

θ<sub>j</sub> = 3° for the jet semi-opening angle, and a reference jet radius R<sub>j</sub> = 50 au at z = 0. As for the jet base density ρ<sub>j0</sub>, we considered a “high-density” value, giving a one-sided jet mass-flux M̄ = 6 × 10<sup>-6</sup> M<sub>☉</sub> yr<sup>-1</sup> typical of very active Class 0 protostars, for our reference values of θ<sub>j</sub> and R<sub>j</sub>.

We then launched seven modified models, each having only one parameter changed with respect to the reference model (H_REF):

(i) H DenSA: with an ambient density scaling ρ<sub>a,0</sub> = 1.6 × 10<sup>-20</sup> g cm<sup>-3</sup> instead of 1.6 × 10<sup>-18</sup> g cm<sup>-3</sup>;

(ii) H THETA: with a jet semi-opening angle θ<sub>j</sub> = 7° instead of 3°;

(iii) H DEnSJ: with a jet base density ρ<sub>j0</sub> 100 times larger than the reference model;

(iv) H VARAMP: with a variability semi-amplitude ΔV = 90 km s<sup>-1</sup> instead of 60 km s<sup>-1</sup>;

(v) H PER: with a variability period P = 300 yr instead of 115 yr;

(vi) H SAWT: with a sawtooth velocity variability profile instead of a sinusoidal one;

(vii) H RAD: with an initial jet radius R<sub>j</sub> = 20 au instead of 50 au.

We note that all the models have the same value of jet mass-flux as the reference model H_REF, except for H_THETA, H_DEnSJ, and H RAD, since M̄ independently varies with θ<sub>j</sub>, ρ<sub>j0</sub>, and R<sub>j</sub> according to Eq. (13). All input parameters for the reference and modified models in this “high-density” case are summarized in Table 2.

To check the robustness of our conclusions, we also computed a second sequence of models with the same parameter