The First Ionization Potential (FIP) effect is an enrichment of heavy elements with low-FIP such as Fe, Si, and Mg compared with photospheric abundances. It was initially measured in the solar wind and Solar Energetic Particles (SEPs) and later inferred from spectroscopic observations of the corona (see Meyer, 1985a, Meyer, 1985b; Bochsler et al., 1986; Gloeckler and Geiss, 1989; Feldman, 1992, and references therein). The FIP bias, i.e., the ratio of coronal to photospheric abundances, is moreover mass independent. This means that processes below the transition region are strongly affecting the hydrostatic balance of the partially ionized chromosphere. Early on, explanations of the FIP effect involved diffusion, flows or some sort of turbulent mixing in the chromosphere to prevent any gravitational settling (see, e.g., Marsch et al., 1995; Schwadron et al., 1999).

Schwadron et al. (1999) favored the hypothesis of turbulent wave heating in the chromosphere as a way to both prevent a mass-dependant fractionation and to obtain a low-FIP bias. Later on, Laming (2004) proposed the ponderomotive acceleration, i.e., the time averaged gradient of electromagnetic fluctuations, as the origin of the FIP effect. The ponderomotive acceleration has this advantage that it may change sign and could explain the inverse FIP effect observed in low-mass stars (Wood and Linsky, 2010; Laming, 2015). Both processes rely on Alfvén waves propagating parallel and anti-parallel to the magnetic field, to trigger a turbulent cascade through non-linear interactions and heating. These wave populations naturally arise in coronal loops, where footprints motions excite the loop on both ends, but they are also expected in the open solar wind where reflection on large scale gradients (Velli et al., 1989; Zhou and Matthaeus, 1989) or compressible instabilities (Tenerani and Velli, 2013; Shoda et al., 2018b; Réville et al., 2018) create an sunward component from a purely anti-sunward wave packet.

The Ulysses spacecraft have shown clear composition differences between the fast and the slow wind components, the slow wind showing FIP biases close to the one observed in coronal loops (Geiss et al., 1995). After all, decades of observations (Belcher and Davis, 1971; Tu and Marsch, 1995; Bruno and Carbone, 2013) and modeling of fluctuations in the solar wind have now brought compelling evidence that turbulence is likely to be a fundamental ingredient of coronal heating and solar wind acceleration (see, e.g., Verdini and Velli, 2007; Perez and Chandran, 2013; Shoda et al., 2018b; Réville et al., 2020b). Perturbations, or waves are also essential to provide the additional acceleration giving birth to the fast wind ( ≥ 600  km/s), through the ponderomotive force (Alazraki and Couturier, 1971; Belcher, 1971; Jacques, 1977; Leer et al., 1982). The low-FIP bias of the slow wind is often understood as the proof that the it originates in coronal loops, through exchange and magnetic reconfiguration (or reconnection) in the low corona (see, e.g., Antiochos et al., 2012). Yet, if turbulence is both responsible for the solar wind acceleration and the FIP effect, can we rule out a scenario where FIP fractionation occurs in purely open regions?

We investigate this question using a coupled modeling approach. First, we extract unidimensional profiles of open solar wind flux tubes and coronal loops using a multi-dimensional MHD code (Réville et al., 2020b). We then use the Shell-Atm code (Buchlin and Velli, 2007; Verdini et al., 2009) to compute the propagation, cascade and dissipation of purely transverse perturbations of velocity and magnetic field along these profiles. We perform a parameter study where we vary the geometry, the initial perturbation amplitude and the initial injection scale, and discuss their effect on the turbulence properties in the chromosphere and transition region. We estimate the resulting FIP biases with an analytical fractionation model.

In this section, we describe the global MHD simulation used to extract the different wind profiles later input in our turbulence model. The code itself is a global MHD solver based on PLUTO (Mignone et al., 2007) and our current implementation is described in details in Réville et al. (2020b) and Réville et al. (2020a). The main purpose of the global simulation is to provide a realistic structure of the transition region, which plays an essential part in the present work. We rely on a single run, of a dipolar solar field of 5 G, which is typical of solar minimum configurations (see, e.g., DeRosa et al., 2012), and a uniform coronal heating prescribed with the following function: Q h = F h / H ( R ⊙ r ) 2 exp ( − r − R ⊙ H ) (1) where H = 1 R ⊙ , and F h = 10 5   erg.cm − 2 s − 1 . The value of F h corresponds to the global energy output of the solar wind.

The result of the simulation is shown in Figure 1. The left panel is a zoom on the transition region as a function of height and latitude. Contrarily with previous studies (Réville et al., 2020a; Réville et al., 2020b), whose inner boundary was located at the base of the corona, we start the simulation at the photosphere with a temperature of 6000 K. We initialize the simulation from a one dimensional profile of the solar atmosphere obtained with the 1D version of the code described in Réville et al. (2018) and the boundary conditions are identical to Réville et al. (2020b). We use a spherical grid with a radial resolution that goes from 3 × 10 − 5 R ⊙ to 0.3 R ⊙ from the photosphere to 20 R ⊙ . There are 544 points in the radial and 512 in the latitudinal direction. The location of the transition region (TR) varies with the latitude and the temperature profiles. It is relatively constant in coronal holes at around 1.003 R ⊙ , i.e., 2000 km above the photosphere. The minimum height corresponds to the highest coronal temperatures at the very edge of the helmet streamer. The height of the TR then increases within the core of the streamer up to 1.007 R ⊙ , accordingly with lower coronal temperatures. These lower temperatures are likely due to increased density that result in increased cooling inside the loop (radiative losses are proportional to the squared density in the model, see Réville et al., 2020b). This structure has already been described in the work of Lionello et al. (2001) using a similar heating function.

The shock capturing ability of PLUTO, based on Riemann solvers, is essential to describe strong density and temperature gradients defining the transition region. Réville et al. (2020b) have developed an additional module to propagate Alfvén waves in the corona. Yet, since this module is based on a Wentzel-Krimers-Brillouin approximation, it is not suited for a precise description of the transition region. To further study the role of turbulence in the composition of the solar wind, we thus rely on the Shell-Atm code, which solves the non-linear incompressible Alfvén wave equations on given profiles along magnetic fields lines, neglecting magneto-acoustic modes. We extracted four different solar wind profiles from our 2.5 D simulation, shown in Figure 1. In the following section, we describe the Shell-Atm simulations made from these profiles.

Shell-Atm is a low dimension fluid turbulence code based on the shell technique (see Giuliani and Carbone, 1998; Nigro et al., 2004), that can model coronal loops (Nigro et al., 2004; Buchlin and Velli, 2007) and open solar wind solutions (Verdini et al., 2009). As fixed background profiles of flow speed, density and temperature, we use the four profiles extracted from our global MHD simulation. They are shown in Figure 2. Cases identified with CH correspond to open field regions in Figure 1. The CH1 profile is well in the coronal hole of the northern hemisphere in the fast wind region. CH2 is located close to the streamer and thus in the slow wind region of the simulation. LO cases are coronal loops, the LO1 case being at the open/close boundary with the minimum height of the transition region and LO2 being the smallest loop, well inside the helmet streamer.

Shell-Atm solves the reduced (incompressible) evolution of Alfvénic fluctuations, i.e. two coupled equations for the evolution of the parallel and anti-parallel wave populations. The transverse components of the fluctuations are discretized in the spectral space, starting from a scale k 0 , and then into 21 other shells at scales k n = 2 n k 0 . We define z n ± = δ v n ∓ δ b n / μ 0 ρ the Elsässer perturbations at scale λ n = λ 0 2 − n . The equations controlling the transport and dissipation of the z ± fields are described in details in Buchlin and Velli (2007) and Verdini et al. (2009). Non-linear interactions are operated in triads following the model of Giuliani and Carbone (1998). Dissipation in Shell-Atm is obtained via explicit resistivity and viscosity. We set the magnetic Prandtl number ν / η to unity and the Lundquist number S = L v A / η lies between 10 6 in the chromosphere and 10 8 in the corona. We write the heating per unit mass created by the cascade: Q ν / ρ = 1 2 ∑ n ν k n 2 ( | z n + | 2 + | z n − | 2 ) (2) Transverse motions’ amplitude δ v ⊙ is forced at the base of the domain in the chromosphere over three shells, starting at the scale k 3 = 8 k 0 . The forcing uses random fluctuation phases that are reset over a correlation time T ⋆ (see Buchlin and Velli, 2007). This is akin to setting a parallel frequency for the Alfvénic perturbations. T ⋆ is chosen around a few hundreds seconds to follow observations of transverse motions in the corona (see the review of Nakariakov and Verwichte, 2005).

We chose to maintain the amplitude δ v ⊙ at the lower boundary, around 100 km above the photosphere, and on both boundaries for coronal loops configurations, using the following condition: z n ± ( t ) = 2 δ v ⊙ − z n ∓ ( t ) (3) This makes the inner boundary condition partially reflective, with a reflection coefficient being a function of time and of the balance of inward and outward wave populations. Imposing a fully reflecting inner boundary conditions leads indeed to an increase of the base velocity perturbations up to unrealistic values.

Next, we define F ± = 1 4 ρ ( v ± v A ) ∑ n [ z n ± | 2 (4) the energy flux of the perturbations at any given location in the domain. For both coronal holes and loops cases, we have the total energy injected per second and surface unit at a given time E tot = F in A 0 − F out A 1 (5) = ( F 0 + + F 0 − ) A 0 − ( F 1 + + F 1 − ) A 1 (6) where the subscripts 0 and 1 denote the bottom and top of the computational domain respectively. In a statistical steady state, we expect to have the losses in the volume compensating for the input energy, i.e.,: 〈 E tot 〉 ∼ 〈 H ν 〉 + 〈 W 〉 (7) where H ν = ∫ Q ν A exp   d s (8) is the integrated turbulent heating in the domain and W is the work of the force exerted by the perturbations on the solar wind flow: W = ∫ ρ ( u ⋅ a w ) A exp   d s (9) The component of the ponderomotive acceleration a w along field lines can be written (Litwin and Rosner, 1998; Laming, 2004; Laming, 2009; Laming, 2012; Dahlburg et al., 2016): a w , s = ∂ ∂ s [ 〈 δ E 2 〉 2 | B | 2 ] (10) where 〈 δ E 〉 is the time averaged electric field due to the perturbations. Equation (10) also assumes that the ion cyclotron frequency is much larger than the wave frequency, which is easily verified in the inner heliosphere for typical Alfvén wave spectra and makes the force mass independent. In the remainder of the text, we will refer to Eq. 10 by ponderomotive force (per unit mass and unit volume), or ponderomotive acceleration indifferently.

As we consider one dimensional profile, we use the normalized area A exp , which is unity at the base of all flux tubes, and we write Eq. 7 in erg.cm − 2 s − 1 .

Table 1 sums up all simulations made with the shell model for this work. For each profile geometry, we explore three sets of turbulence parameters referred as cases A, B, and C. In general, we chose the input parameters to create a turbulent heating close to what is observed in the (open) solar wind and as such used in the MHD described in Background Wind and Transition Region Profiles. We find that a transverse perturbation of δ v ⊙ = 7  km/s yields heating rate close to 10 5   erg.cm − 2 . s − 1 — the minimum value to power the solar wind—in both coronal holes, which defines our reference set of turbulent inputs: cases A. Cases B use a higher forcing amplitude δ v = 10  km/s, which logically results in a higher heating rate. Both cases A and B have the same injection scale L inj 0 = 35000  km, which is the largest scale introduced in the system and corresponds to the size of supergranules. Cases C are set with a injection scale ten times lower, close to the size of granules. To get the heating in the coronal holes at the right order of magnitude, we had to decrease the input velocity perturbations to δ v ⊙ = 3  km/s. All cases have the same correlation time T ⋆ = 600  s.

In Figure 3, we show some properties of the solutions of the shell model runs A for δ v ⊙ = 7  km/s, T ⋆ = 600  s, and L inj 0 = 2 π / ( 8 k 0 ) = 35000  km. As the transverse perturbations are forced over three shells, the actual rms amplitude of the forcing is δ v * = 3 × δ v ⊙ 2 ∼ 12  km/s. These parameters are close to the one used in Verdini et al. (2019) for a typical coronal hole solution. The results shown in Figure 3 and in Table 1 are obtained with a time average of the solution between 150,000 and 200,000 s, which represents around ten Alfvén crossing times for the coronal holes profiles and more for the coronal loops. All simulations have reached a pseudo-steady state during the period. The top panel shows the outward and inward rms fluctuations z ± as a function of the curvilinear abscissa s along the profile. For coronal loops, we show s / L in a logit scale, where L is the total length of the loop, to emphasize the behaviour in the chromosphere and the transition region at both boundaries.

In coronal holes simulations, one striking feature is the regime difference between a mostly balanced turbulence ( z + ∼ z − ) below the transition region and an imbalanced turbulence beyond. This can be seen also in the bottom panel of Figure 3, where the inward and outward spectra E ± ( k n ) = ( z n ± ) 2 are very comparable at s = 1.001 R ⊙ (in blue), with a cascade covering four orders of magnitude and a spectral slope close to the usual Kolmogorov index − 5 / 3 . Higher in the corona, the outward wave dominates clearly while keeping a − 5 / 3 slope, while the inward component has a flatter spectrum with slope close to − 1 . In coronal loops, both populations are everywhere well represented and spectra seems perhaps closer to − 1 slopes. In the middle panel of Figure 3, we show the dissipation profile computed by the shell model. The plain black line is the true heating, while the gray line corresponds to the so-called phenomenological heating, a proxy often used in large scale fluid models (see, e.g. Dmitruk et al., 2002; Verdini and Velli, 2007; Chandran and Hollweg, 2009; Shoda et al., 2018a; Réville et al., 2020b). It can be written: Q p / ρ = 1 2 | z + | 2 | z − | + | z + | | z − | 2 2 L inj 0 (11) As said earlier, the total heating obtained in coronal holes solution H ν is close to the input energy of the global MHD simulations, itself chosen to provide enough power in the solar wind (see Table 1). Nevertheless, the total heating increases as we go from a typical polar coronal hole, to a slow denser wind around the streamer and to coronal loops. This is expected, and increased heating in the low corona is one explanation for the denser slower wind, which originates somewhere close to the loops. Interestingly, we see that the phenomenological heating is in general an overestimation of the true heating in the shell model (as already noted in Verdini et al., 2019), except around the transition region. It does, however, seems to be a reasonable estimate in imbalanced turbulence region, i.e., in the open coronal wind regions.

The properties of cases B are very similar to the one of cases A, only with larger wave amplitudes and larger heating rates. For cases C, the wave amplitude is logically smaller, and the cascade extent is also shorter as we remain with a fixed value of the viscosity/resistivity ν = η = 10 11   cm 2 s − 1 . In all cases, the dissipation starts to be significant for. k n > 10 4 R ⊙ − 1 .

The main objective of this work is to compare the effect of the turbulence properties on the FIP fractionation for different background configurations. Among the key parameter of the current FIP models is the ponderomotive force, or the wave pressure exerted by the perturbations on the background flow. Figure 4 shows the ponderomotive force obtained in all Table 1 cases, close to the inner boundary, in the chromosphere and transition region. Let us first study the cases A of Table 1, shown in Figure 4 in the top panel. The profile of a w , s has in all cases a similar shape, with a slow decrease and a peak located around the transition region. The TR peak, which is responsible for most of the FIP fractionation in the model that follows, has an interesting ordering. The LO2 profile, has usually the highest averaged ponderomotive acceleration and a higher peak, but the maximum of the LO1 is usually of the order of the coronal holes configuration peaks. Moreover, it seems that the amplitude of the peak is a growing function of the height of the TR.

For cases B, we observe for the ponderomotive acceleration essentially a shift up of all the curves, conserving the hierarchy of the previous cases as a function of the geometry. For most cases A and B, the strength of the ponderomotive force is higher than the opposing gravity of the Sun, around ∼ 2 × 10 4   cm.s − 2 , shown in black in Figure 4. The amplitude of the ponderomotive force is thus significant, especially at the transition region and Alfvén waves can have an influence on the coronal abundances. For cases C, however, the first striking feature is the net decrease of the ponderomotive acceleration below the Sun’s gravity pull. This means that although the turbulent heating is perfectly compatible with what is necessary to power the solar wind (compare especially case CH1 A with CH1 C), this set of parameter will likely not create a low-FIP bias through the ponderomotive acceleration. In these last cases, the turbulence injection scale is ten times smaller, i.e., the size of large granules. The initial amplitude is logically lower as the injection is made later in the cascade, which can explain the lower gradient of magnetic fluctuations.

We now compute the FIP fractionation created by the ponderomotive force in our models. We use for this the derivation of Laming (2009), that reads: ρ j cor ρ j chr = exp ( ∫ z chr z cor ξ j a w , s ν eff ν j , i v j 2 d z ) (12) Here, ρ j is the density of a given species j, ξ j the ionization fraction and ν eff = ν j , i ν j , n ξ j ν j , n + ( 1 − ξ j ) ν j , i (13) where ν j , i and ν j , n are the collision frequencies of the species j with ions and neutrals respectively. We use the formulation of Schwadron et al. (1999), Marsch et al. (1995) for the collision frequencies. The ionization fraction ξ j of each heavy ion is computed through an interpolation of Saha equilibria (Saha, 1920; Saha, 1921) for T ≲ 10 4  K and the CHIANTI database and the ChiantiPy interface (Dere et al., 1997) for T ≳ 10 4 K. Finally, v j 2 = c s , j 2 + v T 2 is the (quadratic) sum of the thermal speed and of a mass independent turbulent velocity v T , which represent a wave turbulence heating and is essential to avoid a mass dependant fractionation (see Schwadron et al., 1999).

Figure 5 shows the resulting FIP biases obtained with Eq. 12 and a constant v T = 15  km/s. We represent the minor ion density ratio between the low corona (7000 km above the photosphere) and the chromosphere (700 km above the photosphere), taking the Oxygen as a reference. The top and middle panel of Figure 5 shows a clear low-FIP bias in cases A and B. Elements Fe, Mg, and Si, show bias up to a factor 10 or more compared with Oxygen, in the red case, i.e., the smallest coronal loop. The profile of the FIP bias follows the hierarchy of the ponderomotive force amplitude shown in Figure 4. For CH1 and LO1 cases, the low-FIP bias is weak and roughly similar, while it is generally higher for the CH2 configuration. Finally, as intuited, cases C do not show any clear FIP biases. This is a direct result of the much weaker ponderomotive force obtained with these simulation parameters, and the heavy ion densities are the same in coronal holes and coronal loops.

In Figure 6, we repeat the same exercise, but choosing a turbulent velocity v T = ∑ k n > k ν δ v n 2 (14) We only sum the shell components for k n > k ν = 10 4 R ⊙ − 1 , where dissipation is significant. The average amplitude of v T is around 8 km/s for cases A, 12 km/s for cases B, and 5 km/s for cases C. The choice of v T is crucial to get FIP bias comparable with observations. It needs to be large enough to avoid the gravitational settling, and thus the mass dependence of the abundances in the corona. However, if v T is too large, it will kill the low-FIP bias. The low-FIP bias of cases A is thus slightly enhanced, because the average turbulent velocity given by Eq. 14 is lower than the value used in Figure 5. The FIP bias of case B is quenched and cases C start to show enrichment in low-FIP elements. Moreover, the hierarchy between the different magnetic configurations is much less visible in Figure 6. CH2 cases have for instance, properties very similar to LO2 cases. This shows that a very careful treatment of the turbulence is important to build FIP fractionation models.

In this work, we have combined several models and tools to assess the efficiency of the Alfvénic turbulence in enriching abundances of low-FIP elements in the low corona. First, we used a global MHD model to get realistic profiles of the Alfvén speed and velocity gradients in both coronal holes and coronal loops, starting from the low chromosphere. They are indeed fundamental in the amplification, reflection and non-linear interaction of propagating Alfvén waves. Then, we took advantage of the Shell-Atm code to compute the interaction between counter-propagating Alfvén waves, the turbulent cascade and the dissipation. Our approach certainly lacks self-consistency, as, for instance, the heating obtained with the shell model is significantly larger in coronal loops than in coronal holes (at least for cases A and B), while the MHD simulation has a much more homogeneous heating. It is, however, the first time that both coronal loops and coronal holes are treated with a shell turbulence model. This model solves the fully non-linear incompressible Alfvén wave equations, including the cascading process and the frequency dependant reflections and interactions of Alfvén wave populations. This is an improvement in comparison with analytical non-WKB approach of previous works (see Laming, 2004; Laming, 2009; Laming, 2012). Moreover, the dissipation is treated physically at very high Lundquist numbers, which cannot be achieved in direct numerical simulations (see, e.g., Dahlburg et al., 2016).

A few comments on the applicability of incompressible MHD to this problem are in order. Reduced MHD equations are derived assuming small gradients in the guide field (or even constant B 0 , see Strauss, 1976; Zank and Matthaeus, 1992; Zank and Matthaeus, 1993; Oughton et al., 2017) as well as in the density, conditions that are not verified in the upper chromosphere and the transition region. The excluded nonlinear couplings in the parallel and perpendicular directions (essentially compressible interactions) are responsible for the coupling of Alfvén waves to slow and fast magneto-acoustic modes, leading to steepening, shock formation and dissipation of slow magneto-acoustic modes and fast modes, and the refraction of the latter downwards back to the chromosphere for large perpendicular wavenumbers. The parametric decay instability (PDI), involving compressible processes, is also excluded. The latter plays a role in the formation of the turbulent spectrum and the balance of inward and outward Alfvén wave population, but is effective most in the lower beta regions higher up in the corona (see Tenerani and Velli, 2013; Shoda et al., 2018b; Réville et al., 2018). As noted in Verdini et al. (2019), Shoda et al. (2019), density fluctuations created by the PDI increase the inward Alfvén waves amplitude and eventually increase the turbulent heating in the open wind regions.

Wave steepening, shocks and mode conversion, which do create a cascade in the direction parallel to the magnetic field, have typical spectra ∝ k − 2 and are thus likely to be a secondary dissipation channel in comparison with the dissipation in the perpendicular plane, with spectra close to the Kolmogorov and Iroshnikov-Kraichnan phenomenology (see Figure 3). Alfvén wave coupling with magneto-acoustic modes may accelerate the perpendicular cascade, but because we run the simulation until a pseudo steady-state is reached, this is likely not a strong limitation of our approach. Numerous other modelers have made similar approximations to our own, as for example the studies by van Ballegooijen et al. (2011), Perez and Chandran (2013), van Ballegooijen and Asgari-Targhi (2016), van Ballegooijen and Asgari-Targhi (2017), and Chandran and Perez (2019). This study has no intention to be exhaustive, but we do examine a significant input parameter space in order to extract meaningful insights on the problem.

Our study shows that, assuming that turbulence is the dominant factor in the coronal heating and solar wind acceleration, a ponderomotive force can appear in the chromosphere and the transition region, and can be strong enough to create a low-FIP bias. This depends however on the turbulence parameters. Injecting energy at the scales of super granules provide the wave amplitude necessary for a low-FIP bias comparable with observations (see, e.g. Feldman, 1992). The force is related to the a amplification of the waves in the corona and the strong gradient that appears at the transition region. Nevertheless, if the energy is injected at the scale of granules, the resulting ponderomotive acceleration seems too weak to explain the observations. Recent works have debated of the right injection scale parameters (van Ballegooijen and Asgari-Targhi, 2017; Chandran and Perez, 2019), and the amplitude of the ponderomotive force could be a way of constraining solar wind turbulence models.

A second result of this study is that the low-FIP bias is not exclusive to coronal loops. Interestingly, we obtain significant low-FIP bias for open field configurations along the streamer (CH2). We also have usually little FIP bias for the LO1 loop, i.e., the loop at the very edge of the helmet streamer. This result might seem in opposition with previous works, such as the one of Laming (2004), Laming (2015). Their result is indeed the consequence of resonant Alfvén waves in the loops, likely triggered by reconnection, which amplify significantly the ponderomotive acceleration. Injecting transverse motions from the inner boundary, we do not observe such resonances, and further studies are needed to compare the relative importance of different triggers. We do show however that resonances are not needed to obtain a low-FIP bias in coronal loops and open slow wind regions. This study thus questions, anyhow, the necessity of interchange reconnection to explain the composition of the slow solar wind.

Finally, it is important to stress that our modeling of the solar chromosphere is very simple. We assume, following Laming (2004), Laming (2009) and writing Eq. 12, that proton drag, collisions and turbulent mixing, are exactly compensating for the Sun’s gravity pull in the chromosphere. The actual balance of all these processes remains very hard to estimate. Shocks and compressible fluctuations are also believed to be significant in the chromosphere and even in the transition region (see Carlsson et al., 2019, for a recent review) and are not accounted for in the shell simulations. Ongoing works are dedicated to build a compressible kinetic-fluid model, including heavy ions populations with various charge state and their interactions with protons and waves.