Investigating the Origin of the First Ionization Potential Effect With a Shell Turbulence Model

The enrichment of coronal loops and the slow solar wind with elements that have low First Ionization Potential, known as the FIP effect, has often been interpreted as the tracer of a common origin. A current explanation for this FIP fractionation rests on the influence of ponderomotive forces and turbulent mixing acting at the top of the chromosphere. The implied wave transport and turbulence mechanisms are also key to wave-driven coronal heating and solar wind acceleration models. This work makes use of a shell turbulence model run on open and closed magnetic field lines of the solar corona to investigate with a unified approach the influence of magnetic topology, turbulence amplitude and dissipation on the FIP fractionation. We try in particular to assess whether there is a clear distinction between the FIP effect on closed and open field regions.

Introduction

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 ($\ge 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.

Background Wind and Transition Region Profiles

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\left(\frac{R_{\odot }}{r}\right)}^2\text{exp}\left(-\frac{r-R_{\odot }}{H}\right)(1)$

where $H=1R_{\odot }$, and $F_h={10}^5$ ${\mathrm{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\times {10}^{-5}{R}_{\odot }$ to $0.3R_{\odot }$ from the photosphere to 20 $R_{\odot }$. 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.003R_{\odot }$, 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.007R_{\odot }$, 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.

FIGURE 1

FIGURE 1. Plasma temperature obtained with the global MHD model. On the left panel, we show the temperature profile as a function of latitude and height. The white line shows the contour $T={10}^5$ K, characterizing the location of the transition region. Blue lines and orange lines are open field lines (referred as CH for coronal holes), green and red lines are coronal loops (referred as LO). They are almost perfectly radial at this scale. On the right panel, we show an extended view of the solution and integrated field lines in the meridional plane up to $5R_{\odot }$. The black line on the right is the Alfvén surface.

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 Turbulence Model

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.

FIGURE 2

FIGURE 2. Profiles of the Alfvén speed, flow speed (in dashed lines) and temperature for the four field lines extracted from the simulation, along the curvilinear abscissa s. For loops, we only show one side up to the apex. These profiles will be the source of the amplification and reflection (particularly in open regions) of perturbations, the latter leading to the creation of counter-propagating perturbations. Non-linear interactions between counter-propagating waves create the cascade and the dissipation at small scales.

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^{\pm }=\delta v_n\mp \delta b_n/\sqrt{{\mu }_0\rho }$ the Elsässer perturbations at scale ${\lambda }_n={\lambda }_0{2}^{-n}$. The equations controlling the transport and dissipation of the $z^{\pm }$ 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 $\nu /\eta$ to unity and the Lundquist number $S=L{v}_A/\eta$ 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_{\nu }/\rho =\frac{1}{2}{\displaystyle \sum _n\nu k_n^2}\left({\left|z_n^{+}\right|}^2+{\left|z_n^{-}\right|}^2\right)(2)$

Transverse motions’ amplitude $\delta v_{\odot }$ is forced at the base of the domain in the chromosphere over three shells, starting at the scale $k_3=8k_0$. The forcing uses random fluctuation phases that are reset over a correlation time $T_{\text{⋆}}$ (see Buchlin and Velli, 2007). This is akin to setting a parallel frequency for the Alfvénic perturbations. $T_{\text{⋆}}$ 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 $\delta v_{\odot }$ at the lower boundary, around 100 km above the photosphere, and on both boundaries for coronal loops configurations, using the following condition:

$z_n^{\pm }\left(t\right)=2\delta v_{\odot }-z_n^{\mp }\left(t\right)(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^{\pm }=\frac{1}{4}\rho \left(v\pm v_A\right){\displaystyle \sum _n{[z_n^{\pm }|}^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_{\text{tot}}=F_{\text{in}}{A}_0-F_{\text{out}}{A}_1(5)$

$=\left(F_0^{+}+F_0^{-}\right)A_0-\left(F_1^{+}+F_1^{-}\right)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.,:

$\langle E_{\text{tot}}\rangle \sim \langle H_{\nu }\rangle +\langle W\rangle (7)$

where

$H_{\nu }={\displaystyle \int Q_{\nu }{A}_{\text{exp}}\text{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={\displaystyle \int \rho \left(\mathbf{u}\cdot {\mathbf{a}}_w\right)A_{\text{exp}}\text{d}s}(9)$

The component of the ponderomotive acceleration ${\mathbf{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}=\frac{\partial }{\partial s}\left[\frac{\langle \delta {\mathbf{E}}^2\rangle }{2{\left|B\right|}^2}\right](10)$

where $\langle \delta \mathbf{E}\rangle$ 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_{\text{exp}}$, which is unity at the base of all flux tubes, and we write Eq. 7 in ${\mathrm{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 $\delta v_{\odot }=7$ km/s yields heating rate close to ${10}^5$ ${\mathrm{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 $\delta v=10$ km/s, which logically results in a higher heating rate. Both cases A and B have the same injection scale $L_{\text{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 $\delta v_{\odot }=3$ km/s. All cases have the same correlation time $T_{\text{⋆}}=600$ s.

TABLE 1

TABLE 1. Parameters and outputs for the Shell-Atm simulations.

In Figure 3, we show some properties of the solutions of the shell model runs A for $\delta v_{\odot }=7$ km/s, $T_{\mathrm{\star }}=600$ s, and $L_{\text{inj}}^0=2\pi /\left(8k_0\right)=35000$ km. As the transverse perturbations are forced over three shells, the actual rms amplitude of the forcing is $\delta v_{\text{*}}=\sqrt{3\times \delta v_{\odot }^2}\sim 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^{\pm }$ 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.

FIGURE 3

FIGURE 3. Overview of shell simulations A with the same input parameters on the four geometries of coronal holes and loops. The top panel shows the time averaged rms fluctuations $z^{\pm }$. For coronal loop cases, the curvilinear abscissa s is normalized by L the length of the loop and both boundaries are shown in log scale. The middle panel shows the heating obtained by dissipation at small scales $Q_{\nu }$ in black and compares it with a phenomenological form $Q_p$. The last panel gives the spectra of the fluctuations ($+$ in plain lines, $-$ in dashed lines) at various distances of the Sun. For open regions simulations, the blue, orange, green and red lines corresponds to $r=1.001,1.01,2.5$ and $10R_{\odot }$ respectively. For coronal loops, the blue and orange remain the same, while the green lines are the spectra at the apex of the loop.

In coronal holes simulations, one striking feature is the regime difference between a mostly balanced turbulence ($z^{+}\sim 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^{\pm }\left(k_n\right)={\left(z_n^{\pm }\right)}^2$ are very comparable at $s=1.001R_{\odot }$ (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/\rho =\frac{1}{2}\frac{{\left|z^{+}\right|}^2\left|z^{-}\right|+\left|z^{+}\right|{\left|z^{-}\right|}^2}{2L_{\text{inj}}^0}(11)$

As said earlier, the total heating obtained in coronal holes solution $H_{\nu }$ 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 $\nu =\eta ={10}^{11}$ ${\mathrm{cm}}^2$$s^{-1}$. In all cases, the dissipation starts to be significant for. $k_n>{10}^4{R}_{\odot }^{-1}$.

Ponderomotive Force and First Ionization Potential Fractionation

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.

FIGURE 4

FIGURE 4. Ponderomotive acceleration for all cases of Table 1 in the chromosphere and transition region. Top panel corresponds to cases (A), middle panel to cases (B), with a higher input $\delta v_{\odot }$, and bottom panel to cases (C), with a lower amplitude and higher injection scale. The Sun’s gravity pull is shown in black. The colored dots indicates the position of the transition region in the profiles. They are usually associated with a peak in the ponderomotive acceleration.

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 $\sim 2\times {10}^4$ ${\mathrm{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:

$\frac{{\rho }_j^{\text{cor}}}{{\rho }_j^{\text{chr}}}=\text{exp}\left({\displaystyle {\int }_{z_{\text{chr}}}^{z_{\text{cor}}}\frac{{\xi }_j{a}_{w,s}{\nu }_{\text{eff}}}{{\nu }_{j,i}{v}_j^2}}dz\right)(12)$

Here, ${\rho }_j$ is the density of a given species j, ${\xi }_j$ the ionization fraction and

${\nu }_{\text{eff}}=\frac{{\nu }_{j,i}{\nu }_{j,n}}{{\xi }_j{\nu }_{j,n}+\left(1-{\xi }_j\right){\nu }_{j,i}}(13)$

where ${\nu }_{j,i}$ and ${\nu }_{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 ${\xi }_j$ of each heavy ion is computed through an interpolation of Saha equilibria (Saha, 1920; Saha, 1921) for $T\mathrm{\lesssim }{10}^4$ K and the CHIANTI database and the ChiantiPy interface (Dere et al., 1997) for $T\mathrm{\gtrsim }{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.

FIGURE 5

FIGURE 5. FIP biases according to Eq. 12 in the cases (A)(top panel), (B)(middle panel), and (C)(bottom panel) for usual minor ions. The integration is made between $z_{\text{chr}}=700$ km and $z_{\text{cor}}=7000$ km above the solar surface. The turbulent velocity entering Eq. 12$v_T=15$ km/s is constant.

In Figure 6, we repeat the same exercise, but choosing a turbulent velocity

$v_T=\sqrt{{\displaystyle \sum _{k_n>k_{\nu }}\delta v_n^2}}(14)$

We only sum the shell components for $k_n>k_{\nu }={10}^4{R}_{\odot }^{-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.

FIGURE 6

FIGURE 6. Same as Figure 5 but using $v_T$ as in Eq. 14.

Discussion

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 $\propto 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.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

VR is the main developer of the multi-D MHD model and performed all simulations. AV, ÉB, and MV are the main developers of the SHELL-ATM code. All authors have contributed to the data analysis, text editing and discussions.

Funding

This research was funded by the ERC SLOW SOURCE project (SLOW SOURCE–DLV-819189). MV was supported by the NASA HERMES Drive Science Center grant No. 80NSSC20K0604. We thank M. Laming for very useful discussions. Numerical computations were performed using GENCI grants numbers A0080410133 and A0070410293. This study has made use of the NASA Astrophysics Data System.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Alazraki, G., and Couturier, P. (1971). Solar wind accejeration caused by the gradient of alfven wave pressure. Astron. Astrophys. 13, 380–389.

Google Scholar

Antiochos, S. K., Linker, J. A., Lionello, R., Mikić, Z., Titov, V., and Zurbuchen, T. H. (2012). The structure and dynamics of the corona-heliosphere connection. Space Sci. Rev. 172, 169–185. doi:10.1007/s11214-011-9795-7

CrossRef Full Text | Google Scholar

Belcher, J. W. (1971). Alfvénic wave pressures and the solar wind. Astrophys. J. 168, 509. doi:10.1086/151105

CrossRef Full Text | Google Scholar

Belcher, J. W., and Davis, L. (1971). Large-amplitude Alfvén waves in the interplanetary medium, 2. J. Geophys. Res. 76 (16), 3534–3563. doi:10.1029/ja076i016p03534

CrossRef Full Text | Google Scholar

Bochsler, P., Geiss, J., and Kunz, S. (1986). Abundances of carbon, oxygen, and neon in the solar wind during the period from August 1978 to June 1982. Sol. Phys. 103 (1), 177–201. doi:10.1007/bf00154867

CrossRef Full Text | Google Scholar

Bruno, R., and Carbone, V. (2013). The Solar wind as a turbulence laboratory. Living Rev. Sol. Phys. 10, 2. doi:10.12942/lrsp-2013-2

CrossRef Full Text | Google Scholar

Buchlin, E., and Velli, M. (2007). Shell models of RMHD turbulence and the heating of solar coronal loops. Astrophys. J. 662 (1), 701–714. doi:10.1086/512765

CrossRef Full Text | Google Scholar

Carlsson, M., De Pontieu, B., and Hansteen, V. H. (2019). New view of the solar chromosphere. Annu. Rev. Astron. Astrophys. 57 (1), 189–226. doi:10.1146/annurev-astro-081817-052044

CrossRef Full Text | Google Scholar

Chandran, B. D. G., and Hollweg, J. V. (2009). Alfvén wave reflection and turbulent heating in the solar wind from 1 solar radius to 1 Au: an analytical treatment. Astrophys. J. 707 (2), 1659–1667. doi:10.1088/0004-637x/707/2/1659

CrossRef Full Text | Google Scholar

Chandran, B. D. G., and Perez, J. C. (2019). Reflection-driven magnetohydrodynamic turbulence in the solar atmosphere and solar wind. J. Plasma Phys. 85 (4), 905850409. doi:10.1017/s0022377819000540

CrossRef Full Text | Google Scholar

Dahlburg, R. B., Laming, J. M., Taylor, B. D., and Obenschain, K. (2016). Ponderomotive acceleration in coronal loops. Astrophys. J. 831 (2), 160. doi:10.3847/0004-637x/831/2/160

CrossRef Full Text | Google Scholar

Dere, K. P., Landi, E., Mason, H. E., Monsignori Fossi, B. C., and Young, P. R. (1997). CHIANTI—an atomic database for emission lines. Astron. Astrophys. Suppl. Ser. 125 (1), 149–173. doi:10.1051/aas:1997368

CrossRef Full Text | Google Scholar

DeRosa, M. L., Brun, A. S., and Hoeksema, J. T. (2012). Solar magnetic field reversals and the role of dynamo families. Astrophys. J. 757 (1), 96. doi:10.1088/0004-637x/757/1/96

CrossRef Full Text | Google Scholar

Dmitruk, P., Matthaeus, W. H., Milano, L. J., Oughton, S., Zank, G. P., and Mullan, D. J. (2002). Coronal heating distribution due to low‐frequency, wave‐driven turbulence. Astrophys. J. 575 (1), 571–577. doi:10.1086/341188

CrossRef Full Text | Google Scholar

Feldman, U. (1992). Elemental abundances in the upper solar atmosphere. Phys. Scripta 46 (3), 202–220. doi:10.1088/0031-8949/46/3/002

CrossRef Full Text | Google Scholar

Geiss, J., Gloeckler, G., and von Steiger, R. (1995). Origin of the solar wind from composition Data. Space Sci. Rev. 72 (1–2), 49–60. doi:10.1007/978-94-011-0167-7_8

CrossRef Full Text | Google Scholar

Giuliani, P., and Carbone, V. (1998). A note on shell models for MHD Turbulence. Europhys. Lett. 43 (5), 527–532. doi:10.1209/epl/i1998-00386-y

CrossRef Full Text | Google Scholar

Gloeckler, G., and Geiss, J. (1989). The abundances of elements and isotopes in the solar wind. AIP Conf. Proc. 183 (1), 49–71. doi:10.1063/1.37985

CrossRef Full Text | Google Scholar

Jacques, S. A. (1977). Momentum and energy transport by waves in the solar atmosphere and solar wind. Astrophys. J. 215, 942. doi:10.1086/155430

CrossRef Full Text | Google Scholar

Laming, J. M. (2004). A unified picture of the first ionization potential and inverse first ionization potential effects. Astrophys. J. 614 (2), 1063–1072. doi:10.1086/423780

CrossRef Full Text | Google Scholar

Laming, J. M. (2009). Non-WKB models of the first ionization potential effect: implications for solar coronal heating and the coronal helium and neon abundances. Astrophys. J. 695 (2), 954–969. doi:10.1088/0004-637x/695/2/954

CrossRef Full Text | Google Scholar

Laming, J. M. (2012). Non-WKB models of the first ionization potential effect: the role of slow mode waves. Astrophys. J. 744 (2), 115. doi:10.1088/0004-637x/744/2/115

CrossRef Full Text | Google Scholar

Laming, J. M. (2015). The FIP and Inverse FIP Effects in Solar and Stellar Coronae. Living Rev. Sol. Phys. 12, 2. doi:10.1007/lrsp-2015-2

CrossRef Full Text | Google Scholar

Leer, E., Holzer, T. E., and Flå, T. (1982). Acceleration of the solar wind. Space Sci. Rev. 33, 161–200. doi:10.1007/BF00213253

CrossRef Full Text | Google Scholar

Lionello, R., Linker, J. A., and Mikić, Z. (2001). Including the transition region in models of the large‐scale solar corona. Astrophys. J. 546 (1), 542–551. doi:10.1086/318254

CrossRef Full Text | Google Scholar

Litwin, C., and Rosner, R. (1998). Coronal scale-height enhancement by magnetohydrodynamic waves. Astrophys. J. 506 (2), L143–L146. doi:10.1086/311656

CrossRef Full Text | Google Scholar

Marsch, E., von Steiger, R., and Bochsler, P. (1995). Element fractionation by diffusion in the solar chromosphere. Astron. Astrophys. 301, 261–276.

Google Scholar

Meyer, J.-P. (1985a). The baseline composition of solar energetic particles. Astrophys. J. Suppl. Ser. 57, 173–204.

CrossRef Full Text | Google Scholar

Meyer, J.-P. (1985b). The baseline composition of solar energetic particles. Astrophys. J. Suppl. Ser. 57, 151–171. doi:10.1086/191000

CrossRef Full Text | Google Scholar

Mignone, A., Bodo, G., Massaglia, S., Matsakos, T., Tesileanu, O., Zanni, C., et al. (2007). PLUTO: a numerical code for computational Astrophysics. Astrophys. J. Suppl. Ser. 170 (1), 228–242. doi:10.1086/513316

CrossRef Full Text | Google Scholar

Nakariakov, V. M., and Verwichte, E. (2005). Coronal Waves and Oscillations. Living Rev. Sol. Phys. 2, 3. doi:10.12942/lrsp-2005-3

CrossRef Full Text | Google Scholar

Nigro, G., Malara, F., Carbone, V., and Veltri, P. (2004). Nanoflares and MHD turbulence in coronal loops: a hybrid shell model. Phys. Rev. Lett. 92 (19), 194501. doi:10.1103/physrevlett.92.194501

PubMed Abstract | CrossRef Full Text | Google Scholar

Oughton, S., Matthaeus, W. H., and Dmitruk, P. (2017). Reduced MHD in astrophysical applications: two-dimensional or three-dimensional? Astrophys. J. 839 (1), 2. doi:10.3847/1538-4357/aa67e2

CrossRef Full Text | Google Scholar

Perez, J. C., and Chandran, B. D. G. (2013). Direct numerical simulations of reflection-driven, reduced magnetohydrodynamic turbulence from the Sun to the Alfvén critical point. Astrophys. J. 776 (2), 124. doi:10.1088/0004-637x/776/2/124

CrossRef Full Text | Google Scholar

Réville, V., Tenerani, A., and Velli, M. (2018). Parametric decay and the origin of the low-frequency alfvénic spectrum of the solar wind. Astrophys. J. 866 (1), 38. doi:10.3847/1538-4357/aadb8f

CrossRef Full Text | Google Scholar

Réville, V., Velli, M., Panasenco, O., Tenerani, A., Shi, C., Badman, S. T., et al. (2020a). The role of Alfvén wave dynamics on the large-scale properties of the solar wind: comparing an MHD simulation with parker solar probe E1 Data. Astrophy. J. Suppl. Ser. 246 (2), 24. doi:10.3847/1538-4365/ab4fef

CrossRef Full Text | Google Scholar

Réville, V., Velli, M., Rouillard, A. P., Lavraud, B., Tenerani, A., Shi, C., et al. (2020b). Tearing instability and periodic density perturbations in the slow solar wind. Astrophy. J. Suppl. Ser. 895 (1), L20. doi:10.3847/2041-8213/ab911d

CrossRef Full Text | Google Scholar

Saha, M. N. (1920). Ionisation in the solar chromosphere. Nature 105 (2634), 232–233. doi:10.1038/105232b0

CrossRef Full Text | Google Scholar

Saha, M. N. (1921). On a physical theory of stellar spectra. Proc. Royal Soc. London Ser. A 99, 135–153. doi:10.1098/rspa.1921.0029

CrossRef Full Text | Google Scholar

Schwadron, N. A., Fisk, L. A., and Zurbuchen, T. H. (1999). Elemental fractionation in the slow solar wind. Astrophy. J. 521 (2), 859–867. doi:10.1086/307575

CrossRef Full Text | Google Scholar

Shoda, M., Suzuki, T. K., Asgari-Targhi, M., and Yokoyama, T. (2019). Three-dimensional simulation of the fast solar wind driven by compressible magnetohydrodynamic turbulence. Astrophy. J. 880 (1), L2. doi:10.3847/2041-8213/ab2b45

CrossRef Full Text | Google Scholar

Shoda, M., Yokoyama, T., and Suzuki, T. K. (2018a). A self-consistent model of the coronal heating and solar wind acceleration including compressible and incompressible heating processes. Astrophy. J. 853 (2), 190. doi:10.3847/1538-4357/aaa3e1

CrossRef Full Text | Google Scholar

Shoda, M., Yokoyama, T., and Suzuki, T. K. (2018b). Frequency-dependent Alfvén-wave propagation in the solar wind: onset and suppression of parametric decay instability. Astrophy. J. 860 (1), 17. doi:10.3847/1538-4357/aac218

CrossRef Full Text | Google Scholar

Strauss, H. R. (1976). Nonlinear, three-dimensional magnetohydrodynamics of noncircular tokamaks. Phys. Fluids 19 (1), 134. doi:10.1063/1.861310

CrossRef Full Text | Google Scholar

Tenerani, A., and Velli, M. (2013). Parametric decay of radial Alfvén waves in the expanding accelerating solar wind. J. Geophys. Res. Space Phys. 118 (12), 7507–7516. doi:10.1002/2013ja019293

CrossRef Full Text | Google Scholar

Tu, C.-Y., and Marsch, E. (1995). MHD structures, waves and turbulence in the solar wind: observations and theories. Space Sci. Rev. 73, 1–210. doi:10.1007/BF00748891

CrossRef Full Text | Google Scholar

van Ballegooijen, A. A., and Asgari-Targhi, M. (2017). Direct and inverse cascades in the acceleration region of the fast solar wind. Astrophys. J. 835 (1), 10. doi:10.3847/1538-4357/835/1/10

CrossRef Full Text | Google Scholar

van Ballegooijen, A. A., Asgari-Targhi, M., Cranmer, S. R., and DeLuca, E. E. (2011). Heating of the solar chromosphere and corona by Alfvén wave turbulence. Astrophy. J. 736 (1), 3. doi:10.1088/0004-637x/736/1/3

CrossRef Full Text | Google Scholar

van Ballegooijen, A. A., and Asgari-Targhi, M. (2016). Heating and acceleration of the fast solar wind by Alfvén wave turbulence. Astrophy. J. 821 (2), 106. doi:10.3847/0004-637x/821/2/106

CrossRef Full Text | Google Scholar

Velli, M., Grappin, R., and Mangeney, A. (1989). Turbulent cascade of incompressible unidirectional Alfvén waves in the interplanetary medium. Phys. Rev. Lett. 63 (17), 1807–1810. doi:10.1103/physrevlett.63.1807

PubMed Abstract | CrossRef Full Text | Google Scholar

Verdini, A., Grappin, R., and Montagud-Camps, V. (2019). Turbulent Heating in the Accelerating Region Using a Multishell Model. Solar Phy. 294, 65. doi:10.1007/s11207-019-1458-y

CrossRef Full Text | Google Scholar

Verdini, A., and Velli, M. (2007). Alfvén waves and turbulence in the solar atmosphere and solar wind. Astrophy. J. 662 (1), 669–676. doi:10.1086/510710

CrossRef Full Text | Google Scholar

Verdini, A., Velli, M., and Buchlin, E. (2009). Turbulence in the sub-alfvénic solar wind driven by reflection of low-frequency Alfvén waves. Astrophy. J. 700 (1), L39–L42. doi:10.1088/0004-637x/700/1/l39

CrossRef Full Text | Google Scholar

Wood, B. E., and Linsky, J. L. (2010). Resolving the ξ BOO binary withchandra, and revealing the spectral type dependence OF the coronal “FIP effect”. Astrophy. J. 717 (2), 1279–1290. doi:10.1088/0004-637x/717/2/1279

CrossRef Full Text | Google Scholar

Zank, G. P., and Matthaeus, W. H. (1992). The equations of reduced magnetohydrodynamics. J. Plasma Phys. 48 (1), 85–100. doi:10.1017/s002237780001638x

CrossRef Full Text | Google Scholar

Zank, G. P., and Matthaeus, W. H. (1993). Nearly incompressible fluids. II: magnetohydrodynamics, turbulence, and waves. Phys. Fluid. Fluid Dynam. 5 (1), 257–273. doi:10.1063/1.858780

CrossRef Full Text | Google Scholar

Zhou, Y., and Matthaeus, W. H. (1989). Non-WKB evolution of solar wind fluctuations: a turbulence modeling approach. Geophys. Res. Lett. 16 (7), 755–758. doi:10.1029/gl016i007p00755

CrossRef Full Text | Google Scholar