© The Authors 2024

1. Introduction

Driven by the unprecedented observations of the Parker Solar Probe (PSP; Fox et al. 2016) and Solar Orbiter (Solar Orbiter; Müller et al. 2020; García Marirrodriga et al. 2021) space missions, scientific questions linking solar coronal phenomena and structures with in situ measurements in the inner heliosphere are receiving a tremendous amount of attention (Raouafi et al. 2023). Among those outstanding questions is the possible link between coronal jet-like events with solar wind (SW) switchbacks (SBs).

The SB phenomena have attracted significant attention from the community since the launch of PSP due to their ubiquitousness in the PSP measurements at a few tens of solar radii (Raouafi et al. 2023), while they were only seldomly observed at the Earth orbit and beyond by previous missions (e.g. Yamauchi et al. 2004; Landi et al. 2006; Gosling et al. 2009; Matteini et al. 2014). SBs are intermittent events characterised by a sharp deflection of the magnetic field vector away from the ambient direction and back. While most have deflections < 90° (Dudok de Wit et al. 2020; Bandyopadhyay et al. 2022), some do lead to full reversals, and are called full-reversal SB hereafter. We diverge from definitions given in other studies that restrict SBs to instances with Alfvénic inversion (i.e. magnetic deflections exceeding 90°). This is supported by studies such as Dudok de Wit et al. (2020) and Bizien et al. (2023), which provide evidence that SBs are self-similar structures. This characteristic implies that the properties of SBs–such as the type of boundary or their thickness–are not significantly affected by the magnitude of the magnetic deflection. The SBs are Alfvénic features in which the deflection is most noticeable in the radial magnetic field component (but keeping a roughly constant normal magnetic field intensity |B|∼ constant). They are associated with an increase in the SW radial velocity, v_r. The SBs have thus also been associated with plasma jets (Bale et al. 2019; Kasper et al. 2019). A central open question about SBs is their formation mechanism. Two main families of scenarios have emerged: whether SBs are formed in situ in the SW (Ruffolo et al. 2020; Squire et al. 2020, 2022; Toth et al. 2023), or whether they find their source low in the solar atmosphere (Drake et al. 2021; Wyper et al. 2022; Bale et al. 2023). Building on the latter scenario, the question is then which solar features and events would be the SBs’ progenitors. Solar jet-like features appear as a putative SB source. Similarly to SB, solar jets are both intermittent and ubiquitous phenomena, prominent in the low solar atmosphere. Jets are a middle range aspect of solar activity, in terms of temporal, spatial, and energetic scales. They can be generally defined as impulsive collimated outwardly propagating transient features, observed in the solar atmosphere in emission or absorption in a wide range of spectral lines (Hα through soft X-rays), temperatures (10⁴ − 10⁷ K), and spatial scales (10³ − 10⁸ m). Depending on their size, locations, and waveband of observation, solar jet-like events have historically received a variety of names, such as spicules, photospheric jets, chromospheric jets, surges, jetlets, macrospicules, coronal jets, extreme-ultraviolet (EUV) jets, X-ray jets, among others (cf. reviews of Raouafi et al. 2016; Shen 2021). In terms of Sun-heliosphere relationships, key questions regarding solar jets are whether some by-products of some of these jet-like events are able to reach the inner heliosphere, what would then be the in situ signature of these events, and would those signatures correspond to the SBs’ properties. It is evident that jets observed in the closed corona, such as those present in the vicinity of active regions (Nishizuka et al. 2011), have little chance to reach the SW. There is, however, a tremendous amount of spicules, jetlets, and jets present in coronal holes (e.g. Savcheva et al. 2007; Raouafi et al. 2008; Raouafi & Stenborg 2014; Shimojo & Tsuneta 2009; Kumar et al. 2022; Skirvin et al. 2023), which are thus within the open magnetic field regions magnetically connected to the (fast) SW. Some large-scale jets have been observed to indeed propagate over a few solar radii away from the solar surface (Wang et al. 1998; Patsourakos et al. 2008; Nisticò et al. 2009; Pucci et al. 2013; Moore et al. 2015). The S-shaped structure in coronagraphic data observed by Telloni et al. (2022) may also be the signature of a particularly large jet event propagating upwards. However, such a coronagraphic detection has been relatively infrequent, regarding the jet frequency. Using 70 minute multi-site observations of solar eclipse white light images, hence with an unprecedented image quality and cadence compared to standard coronagraphic images, Hanaoka et al. (2018) show that all EUV jets whose brightnesses are comparable to ordinary soft X-ray jets and that occurred in the polar regions near the eclipse period, could be seen in eclipse images propagating further than 2 R_⊙ (in heliocentric distance). They conclude that ordinary polar jets shall generally reach high altitudes and escape from the Sun as part of the SW. Their results tend to indicate that our capacity to detect jet propagation in coronal images (either direct or by coronagraph) is strongly instrument limited. The absence of smaller jets (e.g. jetlets) being detected by imaging instruments likely does not mean that such events do not propagate far away towards the SW. Based solely on solar observations, some studies speculate that SBs and microstreams can be induced by jet-like events (Neugebauer et al. 1995; Neugebauer 2012; Sterling & Moore 2020; Huang et al. 2023). Eruptions and coronal mass ejections (CME), which, being major perturbations, can be easily tracked from their triggering site in the low solar corona by EUV imagers, through the middle corona and the inner heliosphere thanks to (heliospheric) coronagraphs, and detected by in situ instruments. In contrast, jets, even the largest of them, are too small to be unambiguously associated with a specific in situ feature. To date, the study by Parenti et al. (2021) is the sole study to have achieved a deterministic one-to-one association between specific coronal jet observations and in situ measurements. This underscores the inherent challenges of performing such a direct comparison, even for large jets, particularly in terms of connectivity.

Due to gaps in observations and the community-limited capacity to model, sufficiently precisely, the magnetic connection between a particular region of the heliosphere with its counterpart on the Sun, the link between jet-like events and SBs can, at best, be done statistically, and even this approach is not straightforward. Several authors have compared SBs’ statistical properties during PSP perihelion, trying to link them with the statistical distribution of structures and events in the low solar atmosphere (Fargette et al. 2021; de Pablos et al. 2022; Raouafi et al. 2023; Kumar et al. 2023; Hou et al. 2024; Bizien et al. 2024). Using wavelet transform of magnetic field measurements, Fargette et al. (2021) note that SB occurrence and spectral properties appear to depend on the source region of the SW rather than on the radial distance of PSP. They also found periodic spatial modulations consistent with solar granulation and super granulation, suggesting an influence of the low solar corona on the SB formation mechanism. Bale et al. (2023) show that the mixed-polarity photospheric radial magnetic field distribution is on the same spatial scale as microstreams’ events observed by PSP, suggesting the importance of interchange reconnection in the formation process. Focusing on a corotation interval, Kumar et al. (2023) identified recurring jets, arising from interchange or breakout reconnection at coronal hole bright points and plume bases, as the probable origins of microstreams and SBs observed in the SW. Up to now, the state of the art of numerical modelling is similarly limited in its capacity to deterministically link low corona dynamics with SBs. ’All-inclusive’ simulations of the SW formation are only nascent (Iijima et al. 2023) and are not yet able to resolve all the necessary structures. With two-dimensional (2D) particle-in-cell (PIC) simulations, Bale et al. (2023) point out the potentially important role of interchange reconnection dynamics in explaining SBs. Using three-dimensional (3D) magnetohydrodynamic (MHD) simulations, Wyper et al. (2022) emulated flythroughs within an interchange reconnection-generated curtain of propagating and interacting torsional Alfvénic waves and found Alfvénic patches that closely resemble observations of SBs, although with relatively small magnetic field deflections. In these models, it is directly the properties of interchange reconnection that dictate the properties of the SBs. However, the modelled current sheets are unrealistically large, and how the properties of the Alfvénic patch vary with more realistic coronal current sheets remains to be explored.

Focusing on jet-like events as precursors for SBs, only a few simulations have been produced and analysed (Lionello et al. 2016; Karpen et al. 2017; Szente et al. 2017). This is surprising given the fact that jet-like event modelling is a mature field that has benefited from numerous studies (cf. reviews Raouafi et al. 2016; Shen 2021; Skirvin et al. 2023; Tziotziou et al. 2023). Jet modelling has mostly been focused on understanding the trigger process of jets on the complex dynamics of the multi-thermal and multi-flows’ dynamics, in comparison to coronal observational features, and particle accelerations (see e.g. Pariat et al. 2016; Wyper et al. 2019; González-Avilés et al. 2021; Pallister et al. 2021; Chen et al. 2022; Zhu et al. 2023, for recent models). Jet modelling has demonstrated that the observed features are induced by several simultaneously acting physical processes such as reconnection-driven outflows (exhaust from the reconnection site; e.g. Yokoyama & Shibata 1995; Moreno-Insertis et al. 2008; Archontis et al. 2010), evaporation flows induced by heating (Yokoyama & Shibata 1996; Shimojo et al. 2001; Miyagoshi & Yokoyama 2003, 2004), and the non-linear Alfvénic magnetic untwisting waves induced by reconnection between twisted field lines with untwisted ones (Pariat et al. 2009; Archontis & Hood 2013; Lee et al. 2015; Wyper et al. 2018; Doyle et al. 2019). Three-dimensional MHD models tend to indicate that the magnetic untwisting is the energetically dominant process (Pariat et al. 2009, 2016; Fang et al. 2014). To a large extent, jets can be seen as eruptive processes induced by reconnection-driven full destruction of a twisted structure or flux rope (Pariat et al. 2015; Wyper & DeVore 2016; Wyper et al. 2017).

The importance of the magnetic untwisting mechanism is further enhanced in the jet propagation simulations of Lionello et al. (2016), Karpen et al. (2017, hereafter K17), and Szente et al. (2017), despite their marked differences. The jet propagation simulations of Lionello et al. (2016) modelled the corona using a spherically symmetric, steady-state solar wind resulting from the relaxation of an initial thermodynamic solution (based on Lionello et al. 2009). A series of jets is generated thanks to interchange between an emerging flux twisted flux rope and the ambient open coronal field. Both in Szente et al. (2017) and K17, the jet was impulsively and self-consistently generated following photospheric-boundary shearing motions (based on the model of Pariat et al. 2009). While the simulation of Szente et al. (2017) relies on the two-temperature (protons and electrons) Alfvén Wave Solar atmosphere Model (AWSoM) framework (Holst et al. 2014), in which the atmosphere and SW is self-consistent induced by low-Alfvén wave turbulence, in K17 the atmosphere is more directly initiated with a uni-thermal Parker solar wind solution. Local coronal simulations by Lionello et al. (2016), Pariat et al. (2016), Uritsky et al. (2017), K17, and Szente et al. (2017) have identified that the dominant outcome of jet initiation and the fastest propagating structure is the Alfvénic untwisting magnetic wave. Moreover, Pariat et al. (2016) and Uritsky et al. (2017) found that the propagation of these waves leads to significant complexity, characterised by various structures such as the wave front, a shear Alfvén turbulence domain, shear and compressible turbulence, and ultimately a dense jet. Roberts et al. (2018) carried out an analysis of the possible signatures of the K17’s jet as it was upwardly propagating, and produced synthetic signatures that would be observed by in situ instruments measuring the magnetic and velocity fields. The result being that the untwisting wave front does indeed correspond to a local magnetic deflection, with a decrease in the radial magnetic field intensity, |B_r|, while the norm of the magnetic field |B| remained roughly constant, and that the untwisting wave was also associated with enhanced radial velocity (cf. Fig. 9 of Roberts et al. 2018).

Although only a single jet was analysed, the study by Roberts et al. (2018) suggests that some SBs and microstreams can indeed result from magnetic untwisting waves generated along solar jet-like events. However, the simulation of K17 analysed by Roberts et al. (2018) corresponds to a very specific atmospheric condition in which the atmosphere is magnetically dominated (β ≪ 1) throughout the domain. At the top of the domain, at 9 R_⊙, the value of β in the SW is still lower than 10⁻². The wind is sub-Alfvénic everywhere in the domain. The Alfvén surface, defined as the radially increasing SW speed, equals the local Alfvén speed (Cranmer et al. 2023), and through which the solar wind transition from sub- to super-Alfvénic is thus absent from the setup of K17. Furthermore, the simulation case of K17 corresponds to a rather extreme case of an open corona, in which the Alfvén surface would be tens of solar radii away (at 45 − 50 R_⊙), which is larger than the most recent models and observations of the range of location of the Alfvén surface (cf. Sect. 3 of the review Cranmer et al. 2023, and reference therein).

The objective of the present study is thus to build upon the works of Pariat et al. (2016), K17, and Roberts et al. (2018) to analyse the propagation of coronal jets into the solar wind, by performing parametric simulations using coronal profile that includes an Alfvénic surface inside the simulation domain. Our objective is to determine whether and how the outputs of the self-consistent generation of a coronal jet are influenced by the plasma β properties of the background atmosphere. Such an analysis was partly carried by Pariat et al. (2016), but restricted to a uniform atmosphere (without solar wind and gravity) and in a domain limited to a few fractions of the solar radius above the surface. Pariat et al. (2016) find that the background plasma β significantly influences the outward propagation of the wave, with a clear separation between the Alfvénic magnetic untwisting wave and the dense jet at low β, whereas the two were embedded within each other when the background was uniform and with β ∼ 1. When studying the propagation over several solar radii, the question now is how the propagation of the jet evolves with a radially changing background plasma β.

In the present paper, we want to address the following questions:

• The generation and propagation of self-consistent jets towards the super-Alfvénic wind,

• The influence of different Parker solar wind profiles on jet propagation,

• The investigation into whether the magnetic untwisting wave, triggered by a jet-like event, results in signatures typically associated with SBs or full reversal SBs.

Our paper is organised as follows. In Sect. 2 we first summarise the concept and properties of the numerical experiments. In Sect. 3 we study the dynamics of different propagating structures revealed through parametric simulations. In Sect. 4 we focus on one specific structure of the jet, a leading Alfvénic wave, and examine whether its characteristics align with those of a SB. Finally, in Sect. 5 we summarise our results and discuss them in the broader context of the relation between jets and a SB.

2. Simulation description

This section begins with an introduction to the numerical domain (Sect. 2.1), laying the groundwork for our consecutive analyses. Subsequently, we will outline the initial conditions (Sect. 2.2), setting the stage for three distinct parametric simulations (Sect. 2.2.3). Finally, we will present the first main phases pivotal to the generation of the jet, offering a comprehensive overview of our scientific study (Sect. 2.3).

2.1. Numerical model

The three-dimensional numerical simulations are performed using the Adaptively Refined Magnetohydrodynamic Solver (ARMS, DeVore 1991; DeVore & Antiochos 2008) to solve time-dependent, isothermal, compressible, ideal MHD equations within a spherical coordinate framework, akin to the approach presented in K17 (cf. Sect. 2 and the model therein).

Moreover, the gas follows a ideal gas law for fully ionised hydrogen:

$$\begin{array}{c}\hfill P=2\frac{k_B}{m_{\mathrm{p}}}\rho T\end{array}$$(1)

with thermal pressure P, mass density ρ, proton mass m_p, and the Boltzmann constant k_B. ARMS also evolves the magnetic field, B, and the velocity field, v. We impose a constant and uniform temperature of the plasma, T, and do not solve any temperature equation.

The ARMS code employs the flux-corrected transport (FCT) algorithm derived from DeVore (1991). The FCT algorithm ensures that the magnetic field remains divergence-free within the precision limits of the machine. This approach also prevents the emergence of non-physical outcomes, such as negative mass densities. Additionally, it minimises numerical oscillations associated with intense gradients that occur at the grid scale. This code uses the PARallel Adaptive MESH refinement (PARAMESH) toolkit (MacNeice et al. 2000) that offers parallel adaptive mesh refinement, dynamically adjusting the grid throughout the computation to accommodate the evolving solution.

The domain of the simulation is a spherical wedge that spans radially from 1 R_⊙ (the solar surface) to 15 R_⊙, and ±9° in both co-latitude, θ, and longitude, ϕ, from a central spine axis of the domain of coordinates (θ = 90°, ϕ = 0°; see Fig. 1). On this domain, there is an adaptive mesh with a block-structured and dynamic grid. The computational domain exhibits a radial-to-perpendicular ratio of 45 : 1, with the outer boundary extending significantly beyond the sonic point denoted as r_s, defined in Sect. 2.2.1. The grid is uniformly spaced in angles and exponentially in radius. This arrangement ensures a consistent ∂r/r ratio throughout various levels of refinement. The domain is built upon a framework of 27 × 3 × 3 grid root blocks, each containing 8 × 8 × 8 cells. Up to four extra levels of grid refinement were allowed during the simulation, doubling the spacial resolution at each level (see Fig. 1).

Fig. 1. Simulation domain and magnetic topology. Top left: 3D volume of the domain of the simulation with the root blocks of the grid. Bottom left: A 2D cross section at constant angle of ϕ = 0° of the velocity colour at t = 5500 s. This snapshot highlights the spatial distribution of the velocity. The grid block boundaries are delineated by grey lines. Each block contains 8 × 8 × 8 cells. Right: Initial magnetic topological structure. The isosurface of plasma β = 20 (red spheroid) indicates the location of the 3D magnetic null point. Magnetic field lines of distinct connectivity bounding the separatrix surface are represented, either closed (white field lines) or open to the heliosphere (blue field lines).

A volume, extending from the solar surface (entire coronal base) to a height of 1.77 × 10⁹ cm (2.5 × 10⁻² R_⊙ ), is imposed to refine to the level 4, having a grid resolution of 192 × 192 pixels, corresponding to a spatial resolution about 1.1 × 10⁸ cm (1.6 × 10⁻³ R_⊙) on the solar surface. The same maximum refinement level was imposed on a central region in the volume where the jet forms and propagates, a domain which extends up to 9.8 × 10¹¹ cm (14 R_⊙) in height and over [−3°, +3°] in angle with a spatial resolution going from 1.1 × 10⁸ cm (1.6 × 10⁻³ R_⊙) on the solar surface to 1.6 × 10⁸ cm (2.3 × 10⁻³ R_⊙) at 14 R_⊙. The outermost perimeter of root blocks (within 3° of the side boundaries) was restricted to only one additional level of refinement (level 2) beyond a height of 3.61 × 10⁹ cm (5.2 × 10⁻² R_⊙) above the surface, to avoid overly resolving regions that are near the side boundaries and far from the central, jet-containing portion of the domain. In the transitional zone between these two regions, the blocks undergo refinement up to level 3, ensuring a more gradual transition in refinement levels across the distinct areas (see Fig. 1).

At the radial inner boundary, a line-tying condition is imposed. The tangential velocity, v_b is specified throughout the coronal base and the radial magnetic field, B_r, is specified throughout the domain. In the guard cells beneath the open radial inner boundary, both mass density and pressure have fixed values. At the radial outer boundary, we extrapolate the mass density using the ratio of outer to inner values in the Parker steady wind. All three velocity components are extrapolated, assuming zero-gradient conditions (free flow-through and free slip). The side boundaries (in θ and in ϕ) are closed, reflecting with respect to the normal flow velocities (v_θ and v_ϕ, respectively) and free slip with respect to the tangential flow velocities (v_ϕ, v_θ, and in both cases, v_r). Beyond all six boundaries, the mass density and the magnetic field vector are extrapolated using zero-gradient conditions.

2.2. Initial conditions

This section provides an in-depth analysis of the initial conditions that lay the groundwork for our study. We begin by investigating the hydrodynamic conditions, utilising a solar wind model tailored for a stratified atmosphere (Sect. 2.2.1). Subsequently, we turn our attention to the initial magnetic configuration (Sect. 2.2.2). The intricate interplay between these initial conditions yields three distinct parametric simulations, contributing to a comprehensive understanding of the complex dynamics at play (Sect. 2.2.3).

2.2.1. Solar wind model for a stratified atmosphere

Building on Masson et al. (2013) and K17, we initialised the atmosphere using Parker’s steady, spherically symmetric, supersonic wind model from Parker (1958):

$$\begin{array}{c}\hfill \frac{v_r^2}{c_{\mathrm{s}}^2}exp(1-\frac{v_r^2}{c_{\mathrm{s}}^2})=\frac{r_{\mathrm{s}}^4}{r^4}exp(4-4\frac{r_{\mathrm{s}}}{r})\end{array}$$(2)

where v_r is the radial velocity, $c_{\mathrm{s}}=\sqrt{2k_B{T}_{\mathrm{b}}/m_{\mathrm{p}}}$ is the isothermal sound speed, and r_s = G M_⊙m_p/4k_BT_b is the radius of the sonic point. We assumed a constant temperature for the three parametric simulations : T_b = 10⁶ K (respectively 2 × 10⁶ K for medium and high β simulations) giving a sound speed $c_{\mathrm{s}}=128\mathrm{km}{\mathrm{s}}^{-1}$ (respectively 182 km s⁻¹ for medium and high β simulations) and r_s = 5.7 R_⊙ (respectively 2.9 R_⊙ for medium and high β simulations). This isothermal approximation simplifies the representation of the solar wind but disregards the heating effects caused by reconnection-driven outflows of the jet. As a result, we cannot make precise predictions about the observable phenomena that depend on the detailed thermodynamic behaviour of the coronal plasma. To fully determine the atmosphere, we set the pressure at the bottom boundary, P_b (cf. values in Table 1), which fixes the photospheric mass density thanks to Eq. 1. The density in the rest of the domain is determined by maintaining a constant mass flux, ρ_⊙v_r, ⊙ r_⊙² = cst. This provides an isothermal atmosphere stratified in density with a radial plasma flow.

2.2.2. Initial magnetic configuration

Following Pariat et al. (2009) and K173, to self-consistently model the onset and generation of a jet occurring in small opposite sign magnetic polarity within a predominantly uniform coronal hole, we initially combine two potential magnetic fields. To model the uniform background magnetic field of a coronal hole, we adopted a Sun-centred monopole configuration given by the equation

$$\begin{array}{c}\hfill {\mathit{B}}_{\mathit{m}}=-|B_{\mathrm{m}}|\frac{R_{\odot }^2}{r^2}\widehat{\mathit{r}}\end{array}$$(3)

with $\widehat{\mathit{r}}$, the radial unit vector. The value of the background magnetic field for each simulation can be found in Table 1. While magnetic monopoles are theoretical entities that do not exist in reality, a practical workaround is employed to simulate the open field characteristics of a coronal hole region within a confined computational domain. This involves utilising a monopole to generate a unipolar radially decreasing magnetic field within the simulation volume.

To fit a magnetic concentration within the coronal hole like a parasitic polarity that typically serve as the origin of jets, we used a point dipole located at a depth d = 1 × 10⁷ cm beneath the surface, outside the computational domain, and aligned in the radial direction. These elements result in an anemone-like structure (see Fig. 1, right panel) characterised by a dome-shaped separatrix surface, hosting a 3D null point accompanied by its corresponding fan surface and two spine lines (Longcope 2005; Pariat et al. 2009). This magnetic configuration, featuring an embedded-dipole arrangement, aligns closely with observations of jet-like events (see Nisticò et al. 2009; Schmieder et al. 2013; Liu et al. 2016; Li & Yang 2019; Joshi et al. 2020). The maximum radial field on the surface is B_d = +35 G. This dipole represents a compact embedded polarity placed at the equator (colatitude, θ = 90°) and central meridian (longitude, ϕ = 0°) in our spherical coordinate system to ensure optimal grid resolution.

2.2.3. Parametric simulations

To explore the multifaceted outcomes of this interplay, we systematically varied specific parameters, leading to the derivation of three distinct parametric simulations. The varying parameters and their corresponding values are summarised in the table below:

The varying choices of these parameters result in three distinctly different atmospheric models. One of these simulations, using parameters from the K17 paper, was designated as the reference model, that is, the low β simulation. In order to characterise these new atmospheres, we looked at the spatial evolution of different quantities as a function of the solar radius. This visual analysis is integral to our study as it enables a direct comparison of the simulated atmospheres, highlighting similarities and biases that may offer insights into common jet behaviours. To achieve optimal visualisation of the atmospheric evolution, we extract a specific line corresponding to θ = 97° and ϕ = 7°. Selecting the central part would introduce disturbances caused by the central polarity, which could hinder the clarity of the atmospheric profile. After the relaxation was achieved, we selected key atmospheric parameters, including mass density, magnetic field, plasma β, solar wind velocity and Alfvén velocity. The subsequent paragraphs will delve into a detailed interpretation of Fig. 2, elucidating the observed trends and their implications for our understanding of atmospheric phenomena.

Fig. 2. Comparison of the background radial profiles of diverse quantities in the three parametric simulations: low β (blue curves), medium β (green curves), and high β (orange curves). The figure panels correspond to: mass density, ρ (top left), magnetic field intensity, |B| (top right), solar wind velocity intensity, |vSW| (bottom left), and Alfvén speed, vA (bottom right).

Simulations reveal a consistent mass density pattern, displayed in the top left corner of Fig. 2. Initially, there is a rapid decline within the first solar radii, followed by a slower decrease, indicating an exponential decay trend on the logarithmic scale, thus it is a double exponential curve. Medium β and high β curves exhibit a slower rate of decrease compared to the low β curve, with a consistent gap between the medium and high β curves. Despite different initial mass densities, all simulations converge to the same profile, consistent with theoretical expectations. The observed behaviours align with established models, such as ρv_rr² = cst and v_r ∝ ln(r)^1/2, hence giving the dependence $\rho \propto 1/r^2\sqrt{ln(r)}$. The similarities between the medium and high β curves are attributed to a common temperature, while the differences between the low β curve and the other two may result from temperature variations.

Similarly, the simulations depict a consistent magnetic field pattern, illustrated in the top right corner of Fig. 2. At the solar surface, each curve aligns with the initial value of the background environment. Despite the logarithmic scale, the 1/r² decay trend is evident across all curves, consistent with Eq. (3).

Across all three simulations, the solar wind velocity curves exhibit a similar behaviour characterised by a logarithmic dependence, approximately resembling ln(r)^1/2. This trend aligns with the approximation of the radial velocity v_r using Eq. (2. However, near the solar surface, the behaviour of the low β curve differs. The velocity seems to follow an exponential law: initially, the curve rises slowly, but then its slope increases.

Regarding the Alfvén velocity, the simulations demonstrate a consistent pattern: within the first few solar radii, there is an increase, followed by a subsequent decrease. These variations can be elucidated by the formula for the Alfvén velocity: $v_{\mathrm{A}}=B/\sqrt{{\mu }_0\rho }$, where B ∝ 1/r² and $\rho \propto 1/r^2\sqrt{ln(r)}$, leading to $v_{\mathrm{A}}\propto {ln(r)}^{1/4}/r$. Near r = 1 R_⊙, v_A adheres to the ln(r)^1/4 curve, swiftly transitioning to a 1/r dependence. The disparity between the curves primarily stems from variations in mass density: notably, in the low β simulation, the values at higher solar radii are significantly lower than those of the medium and high β simulations.

Figure 3 displays the radial evolution of the plasma β and explains the rationale behind the names given to each simulation. Across all three simulations, the plasma β = 2ρk_BTμ₀/B² parameter exhibits similar variations: an initial decrease followed by an increase. These variations primarily stem from disparities in mass density values. Notably, the medium and high β simulations display parallel trends in both mass density and plasma β curves. The radius at which the plasma β parameter exceeds unity is particularly significant. It marks the boundary between the region where magnetic pressure dominates (plasma β ≪ 1) and the region where gas pressure dominates over magnetic pressure (plasma β > 1). In the so-called low β simulation (hereafter referred to as Lβ), this curve does not intersect within the domain, indicating that the dynamics are dominated by magnetic pressure and tension. Conversely, in the medium and high β simulations (hereafter referred to as Mβ and Hβ), the curves intersect the plasma β = 1 threshold at radii of 9.6  R_⊙ and 4.8  R_⊙ respectively. Beyond these radii, the dynamics are primarily governed by environmental plasma dynamics rather than magnetic field effects.

Fig. 3. Radial distribution of plasma β for the low β (blue), medium β (green), and high β (orange) simulations.

Figure 4 displays, for each simulation, the radial distribution of the velocity relative to the background Alfvén and sound speeds. Notably, in Lβ, there is no intersection between the local Alfvén velocity and the local solar wind speed, indicating a consistent sub-Alfvénic wind regime. Conversely, in Mβ and Hβ, there exists a radius (effectively a surface) where these velocities converge. Below this threshold, the simulation remains sub-Alfvénic, while above it, the solar wind transitions to a super-Alfvénic regime. These Alfvén surfaces occur at 6.76 and 4.59  R_⊙ respectively for Mβ and Hβ, values consistent with the plasma β thresholds, albeit in proportion.

Fig. 4. Radial distributions of the atmospheric velocity, Alfvén speed, and sonic speed for Lβ (top panel), Mβ (middle panel), and Hβ (bottom panel) simulations. The radius of the Alfvén and sonic surfaces are respectively indicated by dashed and dotted vertical lines.

The study explores the interplay between hydrodynamic and magnetic initial conditions, resulting in three distinct atmospheric simulations. Visual analysis reveals consistent patterns in mass density, magnetic field, and solar wind velocity across simulations, while variations in Alfvén velocity and plasma β highlight the influence of different parameters. These findings underscore the diversity of atmospheric behaviours and offer insights into the dominant factors shaping atmospheric dynamics.

2.3. Simulation steps

In this section we briefly describe the three main phases of the simulation: the initial relaxation phase (Sect. 2.3.1), the bottom boundary driving phase (Sect. 2.3.2) and the jet propagation phase (Sect. 2.3.3).

Before describing the different steps of the simulations, we aim to examine the evolution of energy, focusing particularly on the assessment of each energy component within each simulation and its alignment with prior research findings.

Figure 5 illustrates the distribution of energy components (gravitational, kinetic, internal, and magnetic) for each simulation, along with the combined energy budget for each scenario. The different energy components are defined as:

$$\begin{array}{cc}& E_{\mathrm{mag}}={\displaystyle \int \int {\int }_V}\frac{1}{2{\mu }_0}{\mathit{B}}^2\mathrm{d}V,\hfill \end{array}$$(4)

$$\begin{array}{cc}& E_{\mathrm{kin}}={\displaystyle \int \int {\int }_V}\frac{1}{2}\rho {\mathit{v}}^2\mathrm{d}V,\hfill \end{array}$$(5)

$$\begin{array}{cc}& E_{\mathrm{int}}={\displaystyle \int \int {\int }_V}\frac{P}{\gamma -1}\mathrm{d}V,\hfill \end{array}$$(6)

$$\begin{array}{cc}& E_{\mathrm{grav}}={\displaystyle \int \int {\int }_V}{g}_{\odot }{R}_{\odot }^2\rho r\mathrm{d}r\mathrm{d}\theta \mathrm{d}\varphi \text{.}\hfill \end{array}$$(7)

Fig. 5. Stacked bar chart illustrating the distribution of energy components defined by Eqs. (4)–(7): gravitational (indigo), kinetic (blue), internal (yellow), and magnetic (orange) for the three simulations. Each bar represents the aggregated energy budget for a specific simulation, where different colours denote the contributions of the respective energy components.

The analysis of the figure reveals several noteworthy observations. Firstly, there is a discernible trend wherein the total energy sum appears to increase with higher plasma β values. Particularly, the combined energies exhibit closer proximity in simulations with Mβ to Hβ values. Additionally, distinct trends emerge indicating a decrease in magnetic energy and concurrent increase in internal, gravitational, and kinetic energies as plasma β rises. These trends align with previous findings: Lβ, indicative of a plasma β parameter below 1, exhibits dynamics predominantly governed by magnetic field lines, thereby resulting in magnetic energy dominance over other energy components. In contrast, Mβ and Hβ demonstrate a transition towards dynamics dominated by thermodynamic effects beyond a certain solar radius. Notably, closer proximity of the plasma β = 1 surface (and the Alfvén surface) to the sun, in Hβ, accentuates the prominence of internal, gravitational, and kinetic energy proportions relative to Mβ.

Focusing on Hβ, Fig. 6 shows the temporal evolution of gravitational, kinetic, internal, and magnetic energies. These fluctuations are computed as the energies at time t subtracted by the energies at the initial time (t = 0). The reference time (t = 0) is defined as 100 seconds preceding the commencement of magnetic injection. This plot reveals four distinct phases: (1) a period of relaxation, (2) a gradual accumulation of energy induced by photospheric forcing, (3) a brief, explosive onset of reconnection driven by instability, resulting in energy release and the generation of a jet; and (4) an extensive propagation phase that begins with the jet generation, each of which will be delineated subsequently.

Fig. 6. Temporal evolution of gravitational (indigo), kinetic (blue), internal (yellow), and magnetic (orange) energy variations for the Hβ simulation. Energies are computed relative to their values at time t = 0: (E(t)−E(t = 0)). The reference time, t = 0, is defined as 100 seconds before the onset of magnetic injection.

2.3.1. Relaxation phase

Due to the influence of the solar wind flow, the system initially experiences an imbalance in forces. To attain a quasi-steady state where the magnetic forces and kinetic pressure of the wind counteract this initial imbalance, the system requires a relaxation phase. In our approach, we initially executed the simulation without dynamic grid refinement nor photospheric forcing, allowing the system to relax until both magnetic and kinetic energy reached a nearly constant state. The duration of the relaxation phase varies among simulations. For Lβ, it is only 100 s, while for Mβ it is 71 100 s and for Hβ it is 65 100 s. Figure 6 shows a part of this phase with a period spanning from t = −200 s to t = 100 s, with the reference time, t = 0, corresponding to 100 s before the magnetic injection. We can see in Fig. 6, the energy variations are negligible compared to the variations in the later phases.

2.3.2. Photospheric forcing and twisting motions

In order to energise the system and provide the magnetic energy to the close-field system that will be fuelling the jet, a quasi-steady rotational photospheric flow was applied to the plasma on the bottom surface over the entire parasitic polarity, that is, the one of the embedded dipole. This mimics the slow shearing and rotational flows observed at the magnetic polarities on the Sun. The spatial profile of the photospheric flow, tangential velocity, v_b, depends on the gradient of the radial magnetic field, B_r, and is given by the following equation (K17):

$$\begin{array}{c}\hfill {\mathit{v}}_{\mathit{b}}=\mathit{v}{|}_{r=R_{\odot }}=v_{0}f(t){\lambda }_{\mathrm{b}}\frac{B_2-B_1}{B_r}tanh\left({\lambda }_{\mathrm{b}}\frac{B_r-B_1}{B_2-B_1}\right)\widehat{\mathit{r}}\times {\mathbf{\nabla }}_t{B}_r{|}_{r=R_{\odot }}\end{array}$$(8)

with v₀ = 20 (25 for Hβ) ×10¹² cm²s⁻¹G⁻¹, B₂ and B₁, user defined upper and lower bounds on where the flow is applied. This yields a peak flow of v_b = 88 km s⁻¹ (respectively 59 km.s⁻¹ and 114 km s⁻¹), about 69 % (respectively 32% and 62 %) of the coronal sound speed and only 1.3 % (respectively 0.7% and 1.5%) of the peak Alfvén speed at the surface for the Lβ (respectively Mβ and Hβ).

The velocity flow is applied by adding a temporal ramp profile using the cosine function below:

$$\begin{array}{c}\hfill f(t)=\frac{1}{2}[1-cos\left(\pi \frac{t-t_{\max }}{t_{\max }-t_{\min }}\right)]\end{array}$$(9)

with t ∈ [t_min = 100 s, t_max = 1100 s]. t_min, t_max represent the initial and final time of application of the flow after the end of the relaxation time. The imposed flow v_b is in the clockwise direction and imparts a right-handed twist to the field. Contrary to K17, the flow is smoothly decelerated to rest at t = t_max, similarly to Pariat et al. (2009).

The fan is not directly stressed. The flow applied adheres to the contours of B_r, thus maintaining its surface distributions unchanged over time. Consequently, the potential field remains temporally roughly constant (Linan et al. 2020), facilitating the tracking of the system’s temporal evolution of free magnetic energy, by calculating the difference between the total magnetic energy and its initial value.

During this phase, magnetic energy increases while internal energy decreases (cf. Fig. 6). The imposed clockwise flow induces a counter-clockwise twist to the closed magnetic field lines, boosting the total field strength and reducing the pressure within the expanding separatrix dome. This behaviour explains the magnetic energy’s upward trend and the internal energy’s decline seen in Fig. 6. However, variations in kinetic energy due to tangential flow and separatrix dome expansion, as well as gravitational energy, remain negligible (see Fig. 6).

In Fig. 7, the temporal evolution of the magnetic (top panel) and kinetic energy (bottom) variations are displayed for the three simulations. In the three simulations, during the relaxation, the magnetic energies stay roughly constant. Then, during the photospheric forcing, the magnetic energy increase with the same trend for the three simulations: even though the simulations have three different atmospheres, they are consistent in terms of photospheric forcing, twist injection. Whereas Lβ and Mβ have the same peak value, at 7 × 10²⁸ erg, Hβ exhibits a peak at higher value: 8.4 × 10²⁸ erg. This difference can be explained by the difference of speed of the photospheric flow.

Fig. 7. Temporal evolution of magnetic (top) and kinetic (bottom) energy variations for the three simulations: Lβ (blue), Mβ (green), and Hβ (orange). Energies are computed relative to their values at time t = 0: (E(t)−E(t = 0)). The reference time, t = 0, is defined as 100 seconds before the onset of magnetic injection.

2.3.3. Magnetic reconnection and propagation of the jet

As the energy-build-up phase concludes, the evolution suddenly transitions to an explosively dynamic phase. This phase begins with the onset of instability and the initiation of fast reconnection. Rapid conversion of magnetic energy into kinetic energy follows. The process culminates in the generation of the coronal jet. It is important to note that the jet generation process is entirely self-consistent; the amount of twist is not a parameter we control or choose. The initially symmetrical separatrix dome undergoes asymmetry at t = 900 s, causing the buckling of one side of the separatrix tower. This results in the rapid flattening of the null-point, dispersion of magnetic field into a current patch, and subsequent fragmentation. This process results in the unwinding of the magnetic field lines, effectively mitigating the accumulated twist during the energy build-up, thus elucidating the abrupt decline in magnetic energy variations seen in Fig. 7. Notably, this unwinding entails the reconnection of closed, twisted magnetic field lines with open magnetic field lines, facilitating the propagation of the twist towards both the outer corona and the inner heliosphere.

Throughout this phase, the newly opened magnetic field lines, exemplified by the one depicted in pink in Fig. 8, exhibit substantial bending and entanglement, forming a ‘U-loop’ configuration (defined in Sect. 3.1). The imparted twist to the reconnected open field lines initiates the generation of nonlinear Alfvén waves, which in turn drives plasma outflow through wave-induced pressure gradients. Consequently, at the onset of the instability, kinetic energy escalates (see Fig. 7), marked by the immediate initiation of a wave propagating at a velocity of 100 km s⁻¹ for Hβ. Subsequently, after 1300 s, heightened pressure at the apex of the former dome separatrix triggers the release of plasma outflow, reaching velocities close to 300 km s⁻¹, close to the local Alfvén velocity in Hβ (see Fig. 2). This event initiates the upward propagation of plasma outflows, resulting in a significant increase in pressure. As a consequence, the internal energy, previously in decline, undergoes a sudden and substantial increase (see Fig. 6).

Fig. 8. Snapshots of the evolution of representative magnetic field lines and the high velocity region in Hβ for four different times: 1500 s, 2500 s, 4000 s, and 5500 s. The colour-coded group of magnetic field lines are plotted from fixed points along the bottom boundary from areas where no flow was applied. The bottom boundary corresponding to the solar surface has a size of 18° in both θ and ϕ. In the right column, the isosurface of constant velocity (v = 300 km s−1) is colour-coded according to the ϕ component of the velocity vϕ: blue indicates out-of-plane velocity and red into-plane velocity. Zoomed images are provided for the first two time intervals (1500 s and 2500 s). In the zoomed images, the magnetic fieldlines are visible up to 1.8 R⊙ above the solar surface whereas in the bottom pictures, the magnetic fieldlines are represented up to 4.6 R⊙ above the solar surface. An animation of this figure is available online.

In Fig. 7, we can see that this phase begins with the concurrent reaching of a maximum in the magnetic energy variations and an onset of a steep increase in the kinetic energy variations for the three simulations at t = 1200 s (respectively t = 1900 s and t = 1100 s) for Lβ (respectively Mβ and Hβ). During this phase, the decrease in magnetic energy and steep increase in kinetic energy suggest a rapid conversion of magnetic to kinetic energy and can be seen for the three simulations.

After the magnetic reconnection, the jet is generated from the emergence of nonlinear Alfvén waves and plasma flows triggered by intense kink-driven reconnection events. This jet extends its propagation into both the outer corona and the inner heliosphere. As the wave pressure exerts its force, it compresses the plasma, thereby augmenting the density of the jet material trailing behind the Alfvénic wave front, surpassing the local ambient value of the solar wind. Progressing at the coronal Alfvén speed, which can reach up to 420 km s⁻¹ (Fig. 4), the front can cross a solar radius in about 1300 s, as observed in Hβ. The propagation of the wave pressure results in increased compression of the plasma, consequently causing a sharp rise in internal energy initially. Subsequently, the internal energy becomes stable (see Fig. 6). Furthermore, the compression of the plasma by the wave pressure enhances the density of the jet material behind the Alfvénic wave front, leading to an augmentation of gravitational energy as the denser portions of the jet propagate away from the solar surface. Finally, as the jet progresses, its various segments accelerate, resulting in a increase in kinetic energy variations (see Fig. 7).

This section demonstrates the consistent presence of distinct steps in all three simulations, indicating an auto-consistent generation of the jet. Additionally, it underscores the similarity in the injection of twisting motions across the three simulations. This injection initiates magnetic reconnection, facilitating a partial conversion of magnetic energy into kinetic energy and the jet generation and propagation. Some of the limitations of the model will be discussed in Sect. 5.2. In the next section we delve deeper into the specifics of jet propagation in the three parametric simulations, each characterised by distinct Alfvén surfaces.

3. Comparative analysis of the dynamics of the propagating jet structures in the parametric simulations

In this section we analyse the dynamics of the different propagating structures identified in the parametric simulations. We begin by conducting a detailed analysis and initial identification of the substructures in the Hβ case (Sect. 3.1). Subsequently, we compare the dynamics of these structures in the three different simulations, varying according to the radial atmospheric profiles (Sect. 3.2).

3.1. Identification of structures of the jet

In this section we first focus on the specific case of Hβ, to analyse the jet propagation and dynamics. Figure 8 presents snapshots of the 3D jet propagation at different times with the representation of magnetic field lines and a high velocity region.

We first examine the magnetic field lines in the left columns of Fig. 8. At t = 1500 s, which is 400 s after the onset of kinetic energy increase (see Fig. 6), both the cyan and green field lines are open. The formation and propagation of the untwisting magnetic wave has just started. The green lines appear relatively straight, indicating they have not yet reconnected. They represent the background magnetic environment. In contrast, the cyan lines show some bending, suggesting they reconnected slightly earlier, resulting in their opening. The bending is the signature of a propagating untwisting magnetic wave. The orange lines and most of the pink ones are closed and exhibit significant twisting. The closed pink lines have a section in closer proximity, and orthogonal with a segment of the green lines, in the vicinity of the parasitic polarity. This indicates that they are about to interchange reconnect with each other. One pink line is open, having recently reconnected, and shows pronounced twisting at its lower section. This open pink line also features a ’U-loop’ shape segment. A ’U-loop’ is a segment of a magnetic field line that forms a U-shape, where the bottom of the U is parallel to the surface. It can be defined as a magnetic deflection greater than 90°. This configuration can be clearly observed in the zoomed view on the top left corner of Fig. 8, with the open pink magnetic field line.

At t = 2500 s, the shape of the field lines and the connectivity of some of them have changed, the latter due to ongoing magnetic reconnection. The previously open green lines are now closed, displaying moderate twist. They have undergone interchange reconnection with the orange and pink lines, which are now open. All groups of open lines exhibit a bent section. The cyan and pink lines show moderate twisting, with the bent section of the cyan lines now present further up (1.6  R_⊙). This shift is the main evidence of the propagation of an untwisting wave, as discussed in Pariat et al. (2009, 2016). The newly open orange lines also display a twisted section. Although there appears to be a U-shaped segment in the orange lines, this is merely a projection effect. Other viewpoints (not shown here) confirm that while the orange field lines are twisted, they do not present any magnetic field inversion this time.

At t = 4000 s, further magnetic reconnection is observed: the previously open cyan lines are now closed, while the previously closed green lines are now open. The magnetic deflections along the orange, green, and pink field lines continue to move upwards, signalling the ongoing propagation of an untwisting wave.

At t = 5500 s, such propagation is still ongoing. The bent segment of the orange, green, and pink lines have moved further up, as the result of the ongoing radial propagation of an untwisting wave along their lengths. The orange and green field lines appear as twisted but still do not exhibit U-shaped segments. The bend of the pink lines seems increasingly marked, suggesting that the magnetic deflection is steepening as the magnetic wave front propagates upwards.

The velocity isosurface displayed in the right column of Fig. 8 reveals that the jet exhibits a horizontally rotating structure propagating through the domain over time. A careful examination reveals that the regions with twisted and bent magnetic field lines do correspond to areas of strong transverse velocities. This indicates that the rotating motion of the high-velocity region is associated with the propagation of twisted and bent structures of the magnetic field lines. Therefore, the magnetic deflection of the field lines is closely associated with a velocity enhancement, highlighting the Alfvénic nature of the wave.

Regarding the presence of magnetic inversion, i.e. the presence of U-shaped segment along the field lines, as previously noted, it is observed that immediately after the interchange reconnection of the strongly twisted field lines, the newly formed open field lines display a U-shaped segment in the vicinity of the magnetic polarity (as with the open pink line at t = 1500 s). However, this magnetic inversion disappears rapidly as the untwisting wave propagates upwards. For example, by t = 2500 s, the pink lines no longer present any ‘U-loop’. This observation holds true for all groups of field lines, with none exhibiting U-shaped segments higher than 1.11 R_⊙. This enables us to state that a ’U-loop’ does not survive in the corona.

As previously discussed in Pariat et al. (2009, 2015), Fig. 8 illustrates that the jet generation is a continuous process characterised by the sequential reconnection of magnetic field lines. The newly formed magnetic field lines propagate untwisting torsional waves upwards. While U-shaped loops, indicative of magnetic field inversion, are commonly noted at the coronal base, at the jet source, they quickly disappear as the torsional wave ascends.

To deepen our analysis, we investigate the temporal evolution of various variables as the jet propagates. Given the radial propagation of the jet, we perform a cut along a specific angle within the domain: ϕ = 0°. The 2D cuts at this angle allow us to examine the temporal evolution of the variables, as displayed in Fig. 9. To clearly distinguish the variation of the physical quantities associated with the jet propagation, the background was removed by subtracting the initial values of the quantities (at t = 100 s, i.e. base difference).

Fig. 9. 2D cuts at ϕ = 0° of the evolution of selected quantities in Hβ at t = 2900 s, t = 5400 s, and t = 12 500 s: mass density variations, ρ − ρrel (top row panels), velocity variations v − vrel (middle row panels), and ϕ velocity component, vϕ (bottom row panels), perpendicular to the plane of the cut with blue (respectively red) colour scales indicating observer-directed (respectively oppositely directed) velocities. An animation of this figure is available online.

The base difference is defined as the difference between a quantity at a given time and the same quantity at a previous reference time. In this subsection, the reference time is the end of the relaxation period, t = 100 s. Mathematically, the base difference for a quantity q at time t and position (r, θ, ϕ) is given by:

$$\begin{array}{c}\hfill \delta q(r,\theta ,\varphi ,t)=q(r,\theta ,\varphi ,t)-q_{\mathrm{rel}},\forall \theta \in [81°,99°]\end{array}$$(10)

where q_rel ≡ q(r, θ, ϕ, t = 100 s). As explained in Sect. 2.2.3, the atmosphere is stratified, resulting in increasing background velocities and decreasing mass density with radial distance (cf. Fig. 2). During jet propagation, the apparent velocity of a structure may become comparable to the surrounding environment, preventing the structure from being seen. Hence, removing the background is essential for visualising the different structures. This method is analogous to base-difference imaging in observational data, as used in Kumar et al. (2019) article. We choose a cut at constant ϕ = 0° because the centre of the jet is approximately at this location.

Figure 9 presents the spatial distribution of the mass density variations, δρ = ρ(t)−ρ_rel, the velocity variations, δv = v(t)−v_rel, and the ϕ component of the velocity, v_ϕ. The mass density variations, represented in the top panel of Fig. 9, show that there is a dense structure that extends through the domain. Indeed, at t = 2900 s, a dense structure (plotted in darker red) is present from the solar surface to about 1.97  R_⊙, whereas at t = 5400 s, the structure is extending further to nearly 3 R_⊙. At t = 12 500 s, the dense structure seems to extend only to 2.2 R_⊙ with the colour scale that was selected. However, at t = 18 000 s, the dense structure extends to 11.5 R_⊙ with a density of δρ = 7 × 10⁻¹⁷ kg m⁻³. Moreover, the mass density seems to decrease temporally close to the solar surface: ρ(r = 1.27 R_⊙,  θ = 89.97° , ϕ = 0° , t = 5400 s)−ρ(r = 1.27 R_⊙,  θ = 89.97° , ϕ = 0° , t = 2900 s) = 8.5 × 10⁻¹⁵ kg m⁻³. Thus, the flow of mass density seems to scatter through the simulation domain. In addition to that, it seems that there is a launch of denser plasma at t = 12 500 s: there is a local temporal increase close to the solar surface, almost 4 × 10⁻¹⁵kg m⁻³ between t = 5400 s and t = 12 500 s.

In the middle panels in Fig. 9, the variations in velocity magnitude are depicted. We can observe two different structures. At t = 2900 s, a region with a high increment of velocity (with a maximum value of δv = 357 km s⁻¹) is observed at 1.97 R_⊙, preceded, at a smaller radius, by areas with lower velocity increment (δv ∼ 120 km s⁻¹). This pattern persists for the two other time steps. Indeed, at t = 5400 s, there is a region of higher velocity with a maximum δv = 301 km s⁻¹ near 4 R_⊙ followed by regions of lower δv, around 80 km s⁻¹. At t = 12 500 s, the region of higher velocity with a maximum δv = 108 km s⁻¹ is located at 9.42 R_⊙ followed by regions of lower δv, around 50 km s⁻¹. In addition to that, we can see that the prominent structure with the higher δv has a significant temporal decrease: at t = 5 400 s, the maximum δv is 301 km s⁻¹, whereas at t = 12 500 s, the maximum δv is 108 km s⁻¹. The other areas with lower velocity increment have decreased almost as markedly: a temporal decrease of around 50% between t = 5400 s and t = 12 500 s. By comparing the top and middle panels, we can notice that the areas with a small increase in velocity correspond to areas with an increase in mass density but the high velocity region is on the other hand not correlated with a region of sensitive mass density increment. This will be investigated more in depth in the following paragraphs (see Figs. 10, 11).

Fig. 10. 2D distribution (at ϕ = 0°) of the velocity variations, δv = v − vrel, in green with isocontours for mass density variations, δρ = ρ − ρrel, in red. The isocontours for mass density variations range from 2 × 10−15 to 10−14 kg m−3 (respectively 10−16 to 10−15 kg m−3 and 5 × 10−17 to 8 × 10−16 kg m−3) in Lβ (respectively Mβ and Hβ). On the right part, radial probes (at ϕ = 0° and θ = 90°) are displayed for both quantities. An animation of this figure is available online.

Fig. 11. Radius-time diagrams, at θ = 90° and ϕ = 0°, of the velocity variations v − vrel, the logarithm of the mass density, and the variations of the transverse component of the magnetic field, Bt − Bt, rel, with $B_t=\sqrt{B_{\theta }^2+B_{\varphi }^2}$ for the three different simulations.

The ϕ component of the velocity highlights the rotational motions of the jet. The blue component is predominantly located below θ = 90°, with a small patch above this angle. Conversely, the red component is primarily situated above θ = 90°. Despite their propagation, these patches appear to maintain roughly the same angle locations: the higher positive ϕ velocity component is located above the probe θ = 90° while the lower negative ϕ velocity component is located below the probe θ = 90°. At the latest time step (t = 12 400 s), a new patch of positive velocity appears close to θ = 90°.

When comparing the ϕ component of the velocity and velocity structures, we observe that the leading edge of the high-velocity structure aligns with regions exhibiting the highest rotational velocities. These rotational structures do not appear to correspond with variations in mass density. However, let us note that regions with slower rotational velocities (less than 25 km s⁻¹) are associated with changes in mass density (not discernible in Fig. 9 because of the choice of the range of values).

Figure 9 suggests the presence of multiple distinct structures within the jet, which can be categorised into two types:

• a high velocity, very low mass density structure: this region exhibits high velocity increase without a corresponding significant increase in mass density, and associated with the jet’s rotational component. This pattern may indicate the presence of (near-) incompressible waves. It appears to correspond to areas observed in Fig. 8, where a portion of the isosurface is linked with a high-velocity rotation and a deflection in the magnetic field lines.

• high-mass density, low-velocity structures: these regions have higher mass density but lower velocity, with little or no rotating motions, suggesting a dense bulk flow of plasma.

These results are consistent with findings in Pariat et al. (2016), Uritsky et al. (2017), Roberts et al. (2018). To assess the consistency of these structures across different simulations, we will now investigate if similar substructures are present in the other two simulations.

3.2. Comparative analysis of jet dynamics in three different parametric simulations

In this section we examine the three parametric simulations to determine if the previously observed structures are present in the three simulations. We focus on two quantities: mass density, δρ, and velocity, δv, variations, by removing the background radial profile (see Eq. (10)). In Fig. 10, we compare the spatial distribution of mass density and velocity in order to determine the relative location of the high velocity region and the dense flow region in each simulation. We analyse each simulation when the front of the high velocity structure has reached a given radius of 6.3 R_⊙. In each simulation, this radius is specifically reached at different instants: t = 2800 s (respectively t = 6500 s and t = 7700 s) for Lβ (respectively Mβ and Hβ). This variation in timing is due to the differing speeds at which the front travels, influenced by the properties of the atmospheric stratification. This indicates that the high velocity structure propagates significantly faster in Lβ and slightly faster in Mβ compared to Hβ. To the right of the spatial distribution of mass density and velocity variations, we can see radial probes of these two quantities. In this figure, a radial probe corresponds to the radial distribution of the given quantities along the central spine (θ = 90° and ϕ = 0°). The radial probe effectively reveals the variations in each quantity and their fronts.

We can note in Lβ simulation of Fig. 10, that the structure of increased mass density is primarily located between 1 and 1.8 R_⊙, whereas the high velocity structure extends from 4 R_⊙ to 6.3 R_⊙. Below 4 R_⊙, the increase in plasma velocity is lower than 100 km s⁻¹, while it is above 800 km s⁻¹ in the high velocity front. The radial probe shows that the mass density is higher below 2 R_⊙, location where the velocity enhancement is only moderate. Moreover, between 1.8 R_⊙ and 4 R_⊙, there is a region with minimal mass density and a nearly constant velocity around 100 km s⁻¹. The previously identified substructures are present: a high velocity structure with no increase in mass density and a bulk flow of plasma between 1 and 1.8 R_⊙.

In the Mβ simulation (middle panel of Fig. 10), mass density variations are an order of magnitude lower than in Lβ, spanning from 10⁻¹⁶ to 10⁻¹⁵ kg m⁻³. The high density region reaches 4.3 R_⊙ but thanks to the isocontours, we can observe that the relative mass overdensity decreases radially. The leading edge of the high velocity region has a velocity around 300 km s⁻¹ and is not associated with a significant relative increase in mass density. We can observe that in Mβ, unlike in Lβ, the high velocity region is not located on the central spine (i.e. at θ = 90° and ϕ = 0°), but surrounds it, having a 3D helical shape. In the right part of the middle panel, the radial probe shows, close to the Sun, a high mass density structure with a negative velocity variation relative to the background. Then, between 2.3 and 4.3 R_⊙, there are three mass density variations peaks that coincide with velocity variations peaks. Above 4.3 R_⊙, there is the high velocity structure that coincide with a really low mass density variations region.

For the Hβ case (lower panel of Fig. 10), the velocity increase relative to the background atmosphere can reach up to 214 km s⁻¹ which is close to the values of Mβ. However, the relative mass density increase is about one magnitude lower. The high density region is close to the solar surface and this sensitive increase in mass density can be noted up to 5.5 R_⊙ (faintest isosurface). Some mass density variations isocontours match patches of velocity between 3 and up to 5.5 R_⊙. The leading edge of the high velocity enhancement structure at 6.3 R_⊙ has lower velocities than the previous simulations. Similarly to the two previous β cases, the front of this high velocity region is not associated with a noticeable increase in mass density. Moreover, the high velocity region is again not centred on the spine, but rather adopt a helical shape around the central axis.

Thus, the two substructures identified in Hβ, that are the high velocity front with no mass density enhancement, and the high density region associated with moderate velocity, are also present in the other two simulations. With increasing plasma β between the simulations, the spatial separation between structures decreases.

In order to analyse more precisely the propagation speed of the different structures, Fig. 11 presents radius-time (r–t) diagrams of the base difference (see Eq. (10) of the velocity, the logarithm of mass density and the base difference of the transverse magnetic field (being the square root of the sum of squared θ and ϕ components), for each parametric simulation. These diagrams were built by taking the radial profiles of the quantities along the central spine of the domain (at θ = 90° and ϕ = 0°) for each output of the simulations. This approach means that one radial probe at a given time corresponds to a vertical line in one plot in Fig. 11. These kinds of plots enable the propagation of structures to be seen and for their apparent propagation speed to be estimated directly thanks to the slope of the patterns present in the r–t diagram. This method is commonly used in observations (e.g. Sheeley et al. 1999; Morton et al. 2012; Schmieder et al. 2013).

The top panel of Fig. 11 corresponds to Lβ. The first diagram displays velocity variations and reveals a first structure with a high velocity increase. The mean value of the velocity increase (above the solar wind) is δv ∼ 1600 km s⁻¹. The mean linear speed estimated from the slope indicates an apparent radial propagation speed about 2300 km s⁻¹. In contrast, the mean value of the radial velocity yields v_r ∼ 1 200 km  s⁻¹. The apparent propagation speed is thus significantly larger than the actual radial velocity of the plasma, revealing the wave nature of this leading propagation pattern. The slope in the r-t diagram thus corresponds to the phase speed of the wave and the actual plasma velocity the group speed. This pattern does not appear in the plot of the logarithm of the mass density, suggesting it is thus not associated with a significant density increase. The wave is thus incompressible or weakly compressible. Conversely, this pattern is clearly discernible in the transverse magnetic field increase, highlighting its magnetic nature. All these elements thus suggest that this leading high-propagation-velocity structure is an Alfvénic wave. This structure is the torsional wave induced by the magnetic untwisting mechanism of the jet, as previously discussed by Pariat et al. (2016), Uritsky et al. (2017), Roberts et al. (2018).

In the plot of the velocity variations, we can see a second propagation pattern, having a less inclined slope (about 400 km s⁻¹). This pattern is associated with a weaker plasma velocity increase (above the solar wind), with a mean value of δv ∼ 270 km s⁻¹. The plasma radial velocity, v_r, of this pattern ranges between 190 and 550 km s⁻¹, and with a mean value of 370 km s⁻¹. The real radial plasma velocity and the apparent radial propagation speed are thus corresponding. The velocity is mainly radial with no significant transverse component of the velocity field (not shown here). This propagation pattern appears also in the logarithm of mass density plot, with the same slope (hence same speed), but not in the transverse magnetic field. Thus, we can deduce that this structure is a bulk flow of plasma. This structure was previously considered in Pariat et al. (2016), Uritsky et al. (2017), Roberts et al. (2018).

Moreover, in the plot of the velocity variations, we can see a third structure with a low velocity variation: a mean value around 120 km s⁻¹ above the solar wind and a slope of 80 km s⁻¹. This structure appears also in the logarithm of mass density plot, but not in the transverse magnetic field.

For the other two simulations displayed in the lower panels of Fig. 11, the same structures are present:

• a propagation pattern corresponding to a high velocity structure, having a mean velocity increase, δv, higher than 130 km s⁻¹ above the solar wind (respectively 80 km s⁻¹) for Mβ (respectively Hβ), with radial values of the plasma velocity that range between 140 (respectively 130) and 480 (respectively 530) km s⁻¹. On the r-t diagram, the apparent propagation speed is given by the slope, which value is about 780 km s⁻¹ (respectively 580 km s⁻¹). Moreover, the slope gives an apparent propagation speed close to the Alfvén velocity. This structure is not discernible in the logarithm of the mass density but this structure is also present in the transverse magnetic field increase plot. Thus, this propagation pattern is the signature of a torsional Alfvénic wave induced by the magnetic untwisting mechanism of the jet.

• multiple structures characterised by a velocity variation, δv, of : 80 km s⁻¹ (respectively 70 km s⁻¹) for Mβ (respectively Hβ) above the solar wind. The plasma radial velocity, v_r, of this pattern ranges between 140 (respectively 130) and 520 km s⁻¹ for both simulations, and with a mean value of 480 (respectively 475) km s⁻¹. These structures have a slope giving an apparent propagation velocity of 510 km s⁻¹ (respectively 515 km s⁻¹). These propagation patterns appear also in the logarithm of mass density plot with the same slope, but not in the transverse magnetic field. Hence, we can deduce that these patterns are induced by bulk flows of propagating plasma.

Therefore, in the three simulations there is a leading wavefront, propagating at Alfvénic velocity, near incompressible and associated with intense transverse velocity and magnetic field. This part is the torsional wave induced by the magnetic untwisting of the jet, as previously discussed by Pariat et al. (2016), Uritsky et al. (2017), Roberts et al. (2018). The other structures are bulk flows of plasma, associated with an increased mass density and mainly radially directed with almost no transverse velocity. These findings are consistent with the structures discussed in Pariat et al. (2016), Uritsky et al. (2017), Roberts et al. (2018).

However, among the three simulations, there are some differences: the phase speeds of the wave front decreases for increasing plasma β simulations while the velocity variations of the dense bulk flow of plasma increase for increasing plasma β simulations. The speeds of these structures are influenced by the atmospheric characteristics. Thus, it leads to a propagation ratio of two structures mentioned above decreasing closer to one for increasing plasma β simulations.

4. Switchback in situ signature of a magnetic untwisting jet

Although the U-shaped loops present at the onset of the jet quickly straighten in the lower solar corona (cf. Sect. 3.1), the Alfvénic magnetic deflection continues to propagate radially. Notably, the front part of the jet consists of a torsional Alfvénic wave packet induced by the jet’s magnetic untwisting mechanism, characterised by high rotational speed and transverse magnetic field. This wave front is a potential candidate for inducing typical switchback magnetic deflection signatures when measured in situ. This prompts several key questions: Can the leading part of the jet provide a signature consistent with a switchback ? Furthermore, how does the jet front’s torsional magnetic wave evolves over time?

Figure 12 illustrates the transverse distribution (i.e. at a given radius) of the radial component of the magnetic field, B_r, in the vicinity of the leading edge of the torsional Alfvén wave packet at three different instants during its propagation. At t = 5400 s (respectively at t = 10 100 s and 13 900 s), the wave front is studied at a radius of r = 4.02 R_⊙ (respectively r = 7.73 R_⊙ and r = 10.71 R_⊙), where the background plasma β is equal to 0.76 (respectively 1.94 and 3.45) and the background Alfvénic Mach number, v_SW/v_A, is 0.84 (respectively 1.95 and 2.88). Thus, we study the wave front’s evolution, below, near and above the transition from low-to-high β and the Alfvén surface.

Fig. 12. Distribution at constant radius of the radial magnetic field, Br, in the vicinity of the front of the jet torsional magnetic wave, at three different times: t = 5400 s (top left panel), t = 10 100 s (top right panel), and t = 13 900 s (bottom left panel) for the Hβ simulation. For the one-dimensional probe delineated in blue in the 2D radial magnetic field cut at t = 10 100 s, synthetic data for various quantities B, Br, Bt, and V, Vr, Vt are presented in two separate graphs (bottom right panels).

The colourmap values differed for each radial cut, showing a consistent decrease in the overall intensity of the radial magnetic field, aligned with its definition (cf. Sect. 2.2.3 and Fig. 2). At the three time steps, a similar distribution pattern of the radial magnetic field, B_r, is present. A dark-purple ring like-shape structure of weaker |B_r| is present inside θ ∈ [87° ,93° ], ϕ ∈ [ − 2° ,3° ]. There, the magnetic field has in average 49% (respectively 48% and 42%) of the background radial magnetic field intensity at 5400 s (respectively 10 100 s and 13 900 s). In particular, the weakest radial magnetic field intensity can be found at around θ = 92° and ϕ = 0° for the three instants. However, there are some differences for the distribution of the radial magnetic field for the three instants: the structure narrows after the plasma β equals one, with visible turbulent cells forming in the distribution : darker cells that appear in the background radial magnetic field meaning a weaker radial magnetic field intensity, still higher than the one of the darkest purple structure.

At t = 10 100 s, we produce synthetic in situ like measurements akin to those of the PSP instrument of the FIELDS (Bale et al. 2016) and Solar Wind Electrons Alphas and Protons (SWEAP) (Kasper et al. 2016) suites. The lower right panels of Fig. 12 present the co-latitudinal profiles (at r = 7.73 R_⊙ and θ = 91.35°) of the intensity, radial and transverse components of both the magnetic and velocity fields, respectively normalised by the ambient solar wind values of the wind speed, v_SW, and the magnetic field intensity, B_SW. The later quantities are obtained by taking the mean of the solar wind speed and the magnetic field intensity in the environment, here for ϕ in [−9°, −3°]∪[3°, 9°]:

$$\begin{array}{cc}\hfill q_{\mathrm{SW}}(r,\theta ,t)& =\overline{q}(r,\theta ,\varphi \in [-9°,-3°]\cup [3°,9°],t)\hfill \end{array}$$(11)

$$\begin{array}{cc}\hfill q_{\mathrm{SW}}(r,\varphi ,t)& =\overline{q}(r,\theta \in [81°,87°]\cup [93°,99°],\varphi ,t).\hfill \end{array}$$(12)

The spacecraft trajectories considered in our study assume an idealised probe with unrealistically rapid motion and simplified path characteristics. Specifically, we consider that the spacecraft has infinite speed, allowing it to sample a region spanning 2.43 R_⊙ = 1.7 × 10⁹ m in just one second. Additionally, the trajectory is constrained to vary along only one angle while the other remains constant, simplifying the initially three-dimensional motion to a one dimensional motion along either θ or ϕ angles. This abstraction, although not reflective of actual spacecraft trajectories, enables us to isolate and analyse specific dynamical phenomena without the complexities introduced by realistic motion and temporal limitations.

The synthetic profiles of Fig. 12 show no inversion (change of sign) of the magnetic field, as already found earlier. Despite no full reversal, the absolute value of the normalised radial magnetic field drops from 1 to about 0.5 within 1° in ϕ, then increases back up. The normalised transverse magnetic field concurrently rises from 0 to nearly 1, indicating that the magnetic vector evolves from purely radial to predominantly horizontal. During this magnetic deflection, the value of the norm of the magnetic field undergoes small variations, with the maximum value of |B|/|B_SW| is: 1.12. This is less than 15% of the magnetic field norm of the background solar wind between ϕ = −0.6° and 1.5°. This is consistent with the observations of SB (e.g. Larosa et al. (2011)). Velocity measurements indicate a 40% increase in the radial velocity from the ambient solar wind radial velocity, v_SW, simultaneous with the magnetic field deflection. There is also a simultaneous increase in the transverse component (being the square root of the sum of squared θ and ϕ components) leading to a higher increase in the velocity than the radial component during the magnetic deflection. The velocity enhancement is also consistent with the SB statistics (see Fig. 4 of Larosa et al. 2011). The magnetic deflection thus clearly displays an Alfvénic behaviour, with a link between the magnetic field and velocity field variations. These characteristics are consistent with a switchback signature.

Thus, the leading part of the jet, associated with an Alfvénic wave, exhibits signatures consistent with switchbacks. We estimate the maximum angle of magnetic deflection by plotting the magnetic field vector on a hodogram in the Bt − Br plane (see Fig. 13). The magnetic field vector B, evolves on a sphere of nearly constant field intensity of the background solar wind, |B_SW|. We define the deflection angle as the maximum deviation of the magnetic field vector, B from the radial direction, which represents the background magnetic field direction in our parametric simulations. The deflection angle for the synthetic in situ data of the lower panel of Fig. 9, is approximately 62°. This magnetic deflection, while not representing a full reversal, is thus substantial. Indeed, this deflection stands out from the background. Furthermore, it suggests the possibility of observing deflections without reversal – smaller deflections that exhibit a signature consistent with that of a switchback. This finding, along with other results presented in Fig. 14, indicates that smaller deflections may actually be predominant (see Dudok de Wit et al. (2020)).

Fig. 13. Magnetic field evolution in the Bt − Br plane of the synthetic measurements of Fig. 12. This normalised hodogram displays the change of the magnetic field vector orientation, colour coded according to the varying values of the ϕ coordinate. The quarter of constant ambient solar wind magnetic field intensity, |BSW|, is plotted with an orange dashed curve. The vector of maximum magnetic deflection is represented in grey.

Next, we briefly examine the magnetic deflections across different simulations and their evolution at different radial distances. Figure 14 displays representative hodograms of magnetic deflections as a function of θ (i.e. at given r and ϕ), near the front of the torsional magnetic wave. The hodograms are obtained from both Mβ and Hβ simulations, where the torsional wave propagates through domains with specified plasma β. In the Mβ case, the leading part of the jet, the untwisting magnetic wave crosses layers where the ambient plasma β is ∼0.5 at t = 6500 s. The hodogram taken at r = 6.08 R_⊙ and ϕ = 0.72° and displayed in the upper left panel of Fig. 14, shows a deflection angle of about 42°. Later, the wavefront reaches the height where the background plasma β ∼ 1.5 at around t = 12 800 s. The hodogram at 12.35 R_⊙ and ϕ = 0.57° (cf. Fig. 14 lower left panel) gives a deflection angle of ∼36°. This suggests that in Mβ, the magnetic deflection remains relatively constant, if not diminishing slightly. For the specific condition of this simulation there seems to be now growth of the deflection. In the Hβ case, the leading part of the jet, the untwisting magnetic wave crosses layers where the ambient plasma β is ∼0.5 at t = 3400 s. The hodogram recorded at r = 2.66 R_⊙ and ϕ = 0.63° (upper right panel of Fig. 14) indicates a deflection angle of about 31°. The wavefront subsequently reaches the plasma β ∼ 1.5 layer, around t = 6500 s. The hodogram obtained at 6.52 R_⊙ and ϕ = 0.9° (cf. Fig. 14 lower right panel) shows a deflection angle of ∼38°. Figure 13 shows that when the wavefront is at r = 7.73 R_⊙, corresponding to an ambient plasma β of ∼1.7, the deflection angle is even higher (∼62°). Thus, in Hβ, it seems that the deflection increases as the untwisting wave propagates upwards and encounters layers with higher plasma β. This behaviour is not observed in Mβ.

Fig. 14. Characteristic hodograms of the magnetic field vector in the Bt − Br plane of Mβ (left panels) and Hβ (right panels) simulations at radial layers of two given values of the ambient plasma β: β = 0.5 (top panels) and β = 1.5 (bottom panels). The hodograms are built similarly to the one in Fig. 13, with the difference that the dependence is on θ instead of ϕ.

To conclude, the analysis of magnetic deflection during the jet propagation in three parametric simulations reveals significant insights into its evolution over time and under varying plasma β conditions. The leading part of the jet, associated with an Alfvénic torsional wave packet, demonstrates properties consistent with switchbacks, including a notable deflection in the radial magnetic field and a velocity enhancement. The observed deflection angles vary depending on the plasma β value. For Mβ, the observed deflection angle seems to decrease with increasing ambient plasma β value, whereas for Hβ, the observed deflection angle increases with increasing ambient plasma β value. These findings suggest that the dynamics of magnetic deflection are intricately linked to the plasma environment. A more comprehensive study of the deflection angle evolution in relation to atmospheric profiles, extending beyond the scope of this current study, shall be conducted in a future work. Future studies should also account for the role of the Alfvén speed and Mach number.

5. Conclusions

5.1. Summary

The present study investigates the dynamics of coronal jet propagation into the solar wind using adaptive 3D MHD simulations. Previous research by Pariat et al. (2016) examined the influence of the plasma β on jet propagation, focusing on a domain constrained to a small fraction of the solar radius above the surface and assuming a uniform atmosphere, without gravity and solar wind effects. In contrast, K17 and Roberts et al. (2018) investigated the propagation of a single jet within a stratified atmosphere over an extended domain, although in these studies, the jet remained within the sub-Alfvénic region, within a plasma β ≪ 1 region. Building on these simulations, we conducted a parametric analysis to explore how different coronal radial plasma β profiles can affect the jet propagation in a stratified atmosphere. The jet was generated impulsively and self-consistently through shearing motions at the photospheric boundary, in reference to the model of Pariat et al. (2009). Additionally, the atmosphere was initialised using a uni-thermal Parker solar wind solution (Sect. 2.2.1).

We confirm the results of Pariat et al. (2016) that a self-consistent jet can indeed be generated and propagate in different low corona plasma β conditions (Sect. 2.2.3). After imposing a relaxation period and some photospheric shearing flows, we observed distinct stages such as the magnetic helicity and energy accumulation in the close parasitic polarity, the impulsive release of energy permitted by magnetic reconnection, and the sequential launch of untwisting magnetic torsional waves along newly formed field lines - across all simulations (Sect. 2.3). This suggests a robust and self-consistent generation model for coronal jets. Unlike prior works (Pariat et al. 2016; Karpen et al. 2017; Roberts et al. 2018), which do not consider an Alfvénic surface, we explored the propagation of coronal jets from the sub-Alfvénic wind to the super-Alfvénic wind regimes in two simulations (Lβ simulation serving as a baseline).

As demonstrated by various numerical models (Lionello et al. 2016; Szente et al. 2017, K17), our investigation confirms that the output of coronal jets can propagate easily over tens of solar radii away from the solar surface. Our parametric study shows that jet propagation can occur over a wide range of radial profiles of plasma β, with the Alfvén surface located at different distances, hence encompassing a wide range of conditions existing in the solar atmosphere (West et al. 2023). Notable structures, such as the leading Alfvénic torsional wave and the trailing dense jet, were consistently observed, aligning with previous studies (Pariat et al. 2016; Uritsky et al. 2017; Roberts et al. 2018). The leading structure corresponds to an Alfvénic torsional wave, characterised by a sudden increase in the velocity magnitude and a phase speed close to the local Alfvén speed, with a slight variation in mass density, indicative of its Alfvénic nature and strong rotational speeds. Conversely, patches of lower velocity, but higher mass density variations denote a dense bulk flow of plasma in the trailing dense-jet region (Sect. 3).

While all structures are present in the parametric simulations, we observed differences in their propagation properties depending on the stratification profiles. Specifically, the phase velocity associated with the leading structure decreases from the low to the high β simulations, reflecting differences in the radial ambient Alfvén velocity profiles. Additionally, the radial propagation speed ratio between the torsional magnetic wave and the dense bulk flow diminishes from the low to the high β simulations. As a result, the two structures follow each other more closely with a higher β profile (Sect. 3).

To address whether the magnetic untwisting wave induced by a jet-like event can produce signatures akin to SBs or full reversal SBs, our three parametric simulations consistently show that the magnetic untwisting wave persists along magnetic field lines in the leading part of the jet (Sect. 3.2). Moreover, this jet-induced untwisting torsional magnetic wave exhibits characteristics typical of switchbacks:

• minor variations in the magnetic field norm (less than 15%);

• a deflection in the radial magnetic field, a decrease in the radial component, and an increase in the transverse component; and

• a simultaneous increase in radial velocity.

The angle of magnetic deflection ranges from 17° to 67° away from the ambient magnetic field, which is radial (Sect. 4).

However, the three simulations show that the ‘U-loops’ observed close to the solar surface do not persist higher up in the corona (Sect. 4). Moreover, in none of our parametric simulations were full reversal SBs’ signatures observed, despite the presence of U-shaped loops at the jet onset (Sect. 3).

Furthermore, the evolution of the magnetic deflection varies between simulations: in the Mβ simulation, the deflection angle decreases as plasma β radially increases from 0.5 to 1.5. Conversely, in the Hβ simulation, the deflection angle increases over time, with propagation distance (radius) and ambient solar wind plasma β (Sect. 4). Our parametric study thus deterministically demonstrates the ability of the jet-associated untwisting torsional magnetic wave to exhibit in situ characteristics measurements akin to switchbacks under diverse solar wind conditions.

5.2. Discussion

While our simulations provide valuable insights into the dynamics of the jets and the jet-associated magnetic untwisting waves that induce typical switchbacks signatures, several limitations should be acknowledged.

The first point to consider pertains to the simulation’s setup. In our study, the magnetic sources generating the jets have a spatial width of about 50 Mm, corresponding to large-scale coronal jets. However, preliminary observational statistical studies indicate that such large coronal jets are likely too infrequent to fully explain the occurrence of SBs (Huang et al. 2023). Smaller-scale jet-like events, such as jetlets or spicules (Kumar et al. 2023; Raouafi et al. 2023; Lee et al. 2024), are more likely candidates. To investigate this further, we would need to reduce the size of the jets to determine if the spatial scaling significantly influences the propagation properties and the ability of the jet-like outputs to reach the super-Alfvénic wind. This would require higher spatial resolution simulations, which is challenging due to the already substantial data volume of our runs. However, it is important to note the robust ability of the Pariat et al. (2009) model to generate jets under various conditions, as demonstrated by numerous studies (see Sect. 1). This model’s capability in a broad range of plasma β environments suggests that smaller-scale jets are likely to behave similarly to larger ones, indicating that changing the size may not affect the underlying mechanisms of jet formation.

Another limitation of the present simulations is the reliance on the uniform Parker model. The approach does not include an explicit energy equation, inhibiting modelling of certain physical processes. Our model does not account for other jet-driving mechanisms, such as evaporation flows (Shimojo et al. 2001). Unlike the simulations by Lionello et al. (2016) and Szente et al. (2017), which include more self-consistent solar wind generation mechanisms, our model employs only a simplified background solar wind. The advantage of this approach is the high level of tunability, which allows multiple parametric simulations to be carried out. Future development should incorporate more comprehensive solar wind model, such as polytropic ones (Keppens & Goedbloed 1999).

Additionally, it is important to consider the limitations of the ideal MHD paradigm used in our simulations. The fluid approximations employed are not ideally suited for the plasma conditions in the upper part of the domain. In the upper corona and inner heliosphere, the solar wind becomes collisionless, with particle mean free paths of the order of 1 AU (Verscharen et al. 2019). As a result, kinetic effects, which are known to be significant in the early solar wind, are completely neglected in our simulations.

The MHD treatment is considered both consistent and reliable for modelling the dynamics of large-scale eruptive events in the low solar corona. Our simulations highlight critical aspects of the solar origin scenario for switchbacks. Given that switchbacks are characterised as Alfvénic magnetic deflections with no characteristic deflection angle (Dudok de Wit et al. 2020), a primary outcome of our study is that the untwisting magnetic torsional wave packet induced through the generation of jets exhibits an in situ signature consistent with switchbacks (see Sect. 4). Our simulations deterministically demonstrate a possible link between jet-like events and switchbacks, with magnetic deflections exceeding 62° observed in the super-Alfvénic wind.

However, our analysis remains preliminary. In Sect. 4 we observed distinct radial variations in the properties of the magnetic deflection across two simulations with different radial profiles for β. Further investigation is needed to determine the driving forces behind the propagating magnetic deflection and whether the steepening of the torsional Alfvénic wave is a systematic phenomenon, as suggested by Alfvén wave growth through expansion mechanisms (Squire et al. 2020, 2022; Shoda et al. 2021; Johnston et al. 2022), or if it varies, as indicating by Mβ simulation. Additionally, a more detailed examination of switchback characteristics, such as the nature of their boundaries as rotational or tangential, is warranted (Larosa et al. 2011; Akhavan-Tafti et al. 2021, 2022; Bizien et al. 2023).

Focusing on the formation of full-reversal switchbacks, and in agreement with and complementing the simulations by Wyper et al. (2022), our analysis rules out the direct-injection scenario, where magnetic inversion or ‘U-shaped’ loops are transported directly from the low atmosphere into the super-Alfvénic solar wind. In none of our simulations, which cover a relatively broad range of stratified atmospheric profiles, did we observe the advection of ‘U-shaped’ loops – defined as magnetic deflections greater than 90° – as depicted in several conceptual models of switchback formation (e.g. Fisk & Kasper 2020; Akhavan-Tafti et al. 2022). Although, ‘U-loops’ are present at the onset of jets and during the very early stages of untwisting magnetic waves, the strong Lorentz forces prevalent in the low-β corona immediately straighten these ’U-loops’, preventing their propagation into the upper corona. Thus, our model suggests that full-reversal switchbacks are unlikely to be observed, hence infrequent, in the low-β/sub-Alfvénic solar wind, which is consistent with initial observational analyses of the sub-Alfvénic region reported by PSP (Bandyopadhyay et al. 2022; Jagarlamudi et al. 2023; Cheng et al. 2024).

Our simulations suggest a two-step scenario for the formation of full-reversal SBs. Initially, jet-like events generate propagating torsional magnetic waves within the low-beta corona. These waves induce Alfvénic magnetic deflections, which can directly account for the majority of SBs observed. Some of these deflections may then serve as perturbation seeds for secondary in situ processes, leading to the steepening and eventual reversal of the magnetic field, thereby forming full-reversal SBs. Potential secondary mechanisms include Alfvénic expansion growth (Squire et al. 2020, 2022; Shoda et al. 2021; Johnston et al. 2022) or velocity shears (Toth et al. 2023). This scenario integrates a solar origin for the initial seeds, suggesting that the statistical properties of SBs are influenced by the distribution and frequency of solar-like jets, while in situ mechanisms are responsible for generating full-inversion SBs and determining some of their specific properties.

To better understand the properties, dynamics, and radial evolution of these seed Alfvénic waves, we need to enhance the production of synthetic in situ data from our simulations. In this paper, we assumed that the torsional wave is intersected by an infinitely fast spacecraft. To improve the realism and accuracy of comparisons between synthetic and observed data, we must implement more realistic spacecraft trajectories, similar to those described by Lynch et al. (2021). This likely necessitates rerunning the simulations with a higher graphics output cadence. Such improvements will enable a more accurate validation and refinement of our simulation models against observations from PSP and Solar Orbiter.

Overall, a significant gap in understanding the propagating dynamics of jets arises from observational limitations, including the sampling of the middle corona, spatial resolution, cadence, and access to lower regions (West et al. 2023). New coronagraphic observations, such as those enabled by METIS/Solar Orbiter (Antonucci et al. 2020), can provide valuable new insights. The S-shaped structures observed in METIS data (Telloni et al. 2022) may correspond to signatures of large-scale propagating jets. The study of propagation is particularly important for statistical analyses attempting to link low coronal observations with SB observations. To further investigate these phenomena, observational campaigns could be designed to detect similar signatures near the solar surface with METIS and track their propagation using instruments on PSP and Solar Orbiter. Additionally, the upcoming PROBA-3 (Zhukov 2016) and Polarimeter to UNify the Corona and the Heliosphere (PUNCH) missions (Deforest et al. 2022; Cranmer et al. 2023), with their next-generation coronagraphic capabilities, are expected to offer groundbreaking insights into the propagation of intermediate and small-scale jet events towards the super-Alfvénic wind.

Data availability

Movies associated to Figs 8-10 are available at https://www.aanda.org

Acknowledgments

We thank the anonymous referee for her/his thorough review of the manuscript. JT acknowledges funding from the Initiative Physique des Infinis of Sorbonne Université. This work was supported by the French Programme National PNST of CNRS/INSU, co-funded by CNES and CEA, and by financial support from the French national space agency (CNES) through the APR program. Essential tests and preparatory simulation runs were made possible thanks to the granted access to the HPC resources of MesoPSL, financed by the Region Île-de-France and the project Equip@Meso (reference ANR-10-EQPX-29-01) of the programme Investissements d’Avenir, supervised by the Agence Nationale de la Recherche. This work was also granted access to the HPC resources of IDRIS and CINES under the allocations A0130406331 and A0160406331 made by GENCI, which enabled the production simulation runs. CF acknowledges funding from the CEFIPRA Research Project No. 6904-2. PW was supported by an STFC (UK) consortium grant ST/W00108X/1 and a Leverhulme Trust Research Project grant. The authors thank Rick DeVore, Benjamin Lynch and Judy Karpen for their physical insights and useful discussions.

References

All Tables

All Figures

Fig. 1. Simulation domain and magnetic topology. Top left: 3D volume of the domain of the simulation with the root blocks of the grid. Bottom left: A 2D cross section at constant angle of ϕ = 0° of the velocity colour at t = 5500 s. This snapshot highlights the spatial distribution of the velocity. The grid block boundaries are delineated by grey lines. Each block contains 8 × 8 × 8 cells. Right: Initial magnetic topological structure. The isosurface of plasma β = 20 (red spheroid) indicates the location of the 3D magnetic null point. Magnetic field lines of distinct connectivity bounding the separatrix surface are represented, either closed (white field lines) or open to the heliosphere (blue field lines).

In the text

	Fig. 2. Comparison of the background radial profiles of diverse quantities in the three parametric simulations: low β (blue curves), medium β (green curves), and high β (orange curves). The figure panels correspond to: mass density, ρ (top left), magnetic field intensity, |B| (top right), solar wind velocity intensity, |vSW| (bottom left), and Alfvén speed, vA (bottom right).
In the text

	Fig. 3. Radial distribution of plasma β for the low β (blue), medium β (green), and high β (orange) simulations.
In the text

	Fig. 4. Radial distributions of the atmospheric velocity, Alfvén speed, and sonic speed for Lβ (top panel), Mβ (middle panel), and Hβ (bottom panel) simulations. The radius of the Alfvén and sonic surfaces are respectively indicated by dashed and dotted vertical lines.
In the text

	Fig. 5. Stacked bar chart illustrating the distribution of energy components defined by Eqs. (4)–(7): gravitational (indigo), kinetic (blue), internal (yellow), and magnetic (orange) for the three simulations. Each bar represents the aggregated energy budget for a specific simulation, where different colours denote the contributions of the respective energy components.
In the text

	Fig. 6. Temporal evolution of gravitational (indigo), kinetic (blue), internal (yellow), and magnetic (orange) energy variations for the Hβ simulation. Energies are computed relative to their values at time t = 0: (E(t)−E(t = 0)). The reference time, t = 0, is defined as 100 seconds before the onset of magnetic injection.
In the text

	Fig. 7. Temporal evolution of magnetic (top) and kinetic (bottom) energy variations for the three simulations: Lβ (blue), Mβ (green), and Hβ (orange). Energies are computed relative to their values at time t = 0: (E(t)−E(t = 0)). The reference time, t = 0, is defined as 100 seconds before the onset of magnetic injection.
In the text

Fig. 8. Snapshots of the evolution of representative magnetic field lines and the high velocity region in Hβ for four different times: 1500 s, 2500 s, 4000 s, and 5500 s. The colour-coded group of magnetic field lines are plotted from fixed points along the bottom boundary from areas where no flow was applied. The bottom boundary corresponding to the solar surface has a size of 18° in both θ and ϕ. In the right column, the isosurface of constant velocity (v = 300 km s−1) is colour-coded according to the ϕ component of the velocity vϕ: blue indicates out-of-plane velocity and red into-plane velocity. Zoomed images are provided for the first two time intervals (1500 s and 2500 s). In the zoomed images, the magnetic fieldlines are visible up to 1.8 R⊙ above the solar surface whereas in the bottom pictures, the magnetic fieldlines are represented up to 4.6 R⊙ above the solar surface. An animation of this figure is available online.

In the text

Fig. 9. 2D cuts at ϕ = 0° of the evolution of selected quantities in Hβ at t = 2900 s, t = 5400 s, and t = 12 500 s: mass density variations, ρ − ρrel (top row panels), velocity variations v − vrel (middle row panels), and ϕ velocity component, vϕ (bottom row panels), perpendicular to the plane of the cut with blue (respectively red) colour scales indicating observer-directed (respectively oppositely directed) velocities. An animation of this figure is available online.

In the text

Fig. 10. 2D distribution (at ϕ = 0°) of the velocity variations, δv = v − vrel, in green with isocontours for mass density variations, δρ = ρ − ρrel, in red. The isocontours for mass density variations range from 2 × 10−15 to 10−14 kg m−3 (respectively 10−16 to 10−15 kg m−3 and 5 × 10−17 to 8 × 10−16 kg m−3) in Lβ (respectively Mβ and Hβ). On the right part, radial probes (at ϕ = 0° and θ = 90°) are displayed for both quantities. An animation of this figure is available online.

In the text

	Fig. 11. Radius-time diagrams, at θ = 90° and ϕ = 0°, of the velocity variations v − vrel, the logarithm of the mass density, and the variations of the transverse component of the magnetic field, Bt − Bt, rel, with $B_t=\sqrt{B_{\theta }^2+B_{\varphi }^2}$ for the three different simulations.
In the text

Fig. 12. Distribution at constant radius of the radial magnetic field, Br, in the vicinity of the front of the jet torsional magnetic wave, at three different times: t = 5400 s (top left panel), t = 10 100 s (top right panel), and t = 13 900 s (bottom left panel) for the Hβ simulation. For the one-dimensional probe delineated in blue in the 2D radial magnetic field cut at t = 10 100 s, synthetic data for various quantities B, Br, Bt, and V, Vr, Vt are presented in two separate graphs (bottom right panels).

In the text

Fig. 13. Magnetic field evolution in the Bt − Br plane of the synthetic measurements of Fig. 12. This normalised hodogram displays the change of the magnetic field vector orientation, colour coded according to the varying values of the ϕ coordinate. The quarter of constant ambient solar wind magnetic field intensity, |BSW|, is plotted with an orange dashed curve. The vector of maximum magnetic deflection is represented in grey.

In the text

	Fig. 14. Characteristic hodograms of the magnetic field vector in the Bt − Br plane of Mβ (left panels) and Hβ (right panels) simulations at radial layers of two given values of the ambient plasma β: β = 0.5 (top panels) and β = 1.5 (bottom panels). The hodograms are built similarly to the one in Fig. 13, with the difference that the dependence is on θ instead of ϕ.
In the text