Inward-propagating Plasma Parcels in the Solar Corona: Models with Aerodynamic Drag, Ablation, and Snowplow Accretion

The complete equation of motion for a small coronal parcel would require the inclusion of forces from gravity, the gas pressure gradient, the wave pressure gradient (Bretherton & Garrett 1968; Belcher 1971), magnetic mirror effects (if the velocity distribution is anisotropic with respect to the magnetic field), MHD Lorentz force terms, and particle collisions. For both the parcel and the surrounding solar wind, many of these forces are expected to be comparable to one another in magnitude. Thus, as in studies of CME kinematics (e.g., Cargill 2004), we will account only for relative differences in forcing between the parcel and its surroundings. By specifying the acceleration of the ambient solar wind as a kind of attractor for the parcel flow, we can neglect the forces that maintain it in its assumed steady state.

Hydrodynamic and MHD forces of interaction between the parcel and the surrounding solar wind are often referred to as drag forces. When studying the motions of CMEs, there has been some debate about whether this force scales linearly with the relative velocity between the two fluids (as in high-viscosity Stokes drag; see Vršnak 2001), or whether it scales quadratically (as in low-viscosity turbulent drag; see Borgazzi et al. 2009). Comparison with MHD simulations tends to favor quadratic drag (e.g., Cargill et al. 1996; Cargill 2004), so that formalism is used here. Consider a finite parcel with mass M_i , positioned at a heliocentric distance r_i and flowing radially with velocity u_i . Its equations of motion are

$$\begin{eqnarray}&&\displaystyle \frac{{{dr}}_{i}}{{dt}}={u}_{i},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&({M}_{i}+{M}_{\mathrm{add}})\frac{{{du}}_{i}}{{dt}}=-{\rho }_{e}{{AC}}_{{\rm{D}}}({u}_{i}-{u}_{e})| {u}_{i}-{u}_{e}| ,\end{eqnarray} \tag{ 2 }$$

where subscript i refers to the internal properties of the parcel and subscript e refers to the external properties of the surrounding corona (see also Chen 1989, 1996; Subramanian et al. 2012; Vršnak et al. 2013; Dolei et al. 2014; Dumbović et al. 2018; Kay et al. 2020). Thus, u_e and ρ_e are specified as the speed and mass density of the ambient, time-steady solar wind. From the form of Equation (2), it is clear that the asymptotic solution is for u_i to approach u_e in the limit of t → ∞ .

Figure 1 illustrates a simplified parcel geometry, in which the volume is given by the product of a radial length scale L and a cross-sectional area A. The mass interior to the parcel is written as

$$\begin{eqnarray}&&{M}_{i}={\rho }_{i}{AL},\end{eqnarray} \tag{ 3 }$$

where all three components—including the parcel's mean mass density ρ_i —can be allowed to vary with radial distance. For the models discussed in this section, it is assumed that M_i remains constant in time. Note that the left-hand side of Equation (2) also contains an added mass (or virtual mass) term

$$\begin{eqnarray}&&{M}_{\mathrm{add}}={C}_{{\rm{A}}}{M}_{e},\end{eqnarray} \tag{ 4 }$$

where M_e = ρ_e AL is the mass of adjacent solar wind plasma with the same volume as the moving parcel. The dimensionless coefficient C_A is often assumed to be equal to 1/2, and it is discussed in more detail below. The M_add term accounts for the fact that a finite-sized parcel induces motions in the surrounding plasma and thus carries along a bit of extra inertia. For additional applications of this concept to solar flux tubes, see Spruit (1981), Ryutova & Priest (1993), and Cranmer & van Ballegooijen (2005).

Zoom In Zoom Out Reset image size

Much of the physics describing momentum transfer between the parcel and its surroundings is encapsulated in the dimensionless drag coefficient C_D. Computational and observational studies of CMEs have tried to use the large-scale motions of flux ropes to determine likely values for C_D in the corona and inner heliosphere. For example, Cargill et al. (1996) found that varying the magnetic field geometry in MHD simulations caused C_D to vary between values of 0 and 3. In this paper, where parcels are tracked primarily between about 7 and 12 R_⊙ (i.e., only over a relatively small dynamic range in distance), C_D is treated as an arbitrary constant. In situations where the parcel is interpreted as a compact MHD plasmoid, an alternate approach could be to take account of the full set of diamagnetic forces, including a "melon-seed" effect that accelerates the object in the direction of the weaker external field (see, e.g., Parker 1957; Schlüter 1957; Pneuman & Cargill 1985; Mullan 1990; Lin et al. 2008).

Realistic values for the added mass coefficient C_A can be obtained from hydrodynamics. The standard value of 1/2, used in many drag-based CME models, is appropriate for a spherical object embedded in an incompressible fluid. However, the parcels considered here may become highly anisotropic, i.e., either prolate (longer in the radial direction than in the transverse direction) or oblate (shorter in the radial direction than in the transverse direction). Thus, to take one step beyond the assumption of sphericity, we utilize results for spheroidal objects in incompressible flows (e.g., Lamb 1932; Morel 2015; Fitzpatrick 2017). For a spheroid moving along the direction of its primary axis of symmetry, the total length along its trajectory is specified as L_∥, and its transverse diameter (i.e., the diameter of the circle one would see by looking down its axis) is L_⊥. Defining the aspect ratio α = L_∥/L_⊥, the spheroid is oblate if α < 1 and prolate if α > 1.

Figure 2 shows the result of an ideal hydrodynamics calculation for C_A as a function of α. The exact functions are lengthy but are given in full in Equations (7.139)–(7.140) and (7.162)–(7.163) of Fitzpatrick (2017). In the limit of α ≪ 1, the analytic result reduces to C_A ≈ 2/(π α). Note that C_A → ∞ in this case because the parcel becomes equivalent to a flat disk whose volume approaches zero. The added mass M_add remains finite and is given by 2ρ_e /π times the volume of an equivalent sphere that circumscribes the disk exactly. In the limit of α ≫ 1, the analytic result reduces to ${C}_{{\rm{A}}}\approx \mathrm{ln}(4{\alpha }^{2})/(2{\alpha }^{2})$, which approaches zero as α → ∞ . For the parameters shown in Figure 2, the following fitting formula always agrees with the exact calculation to within 5%:

$$\begin{eqnarray}&&{C}_{{\rm{A}}}\approx \frac{2}{\pi \alpha (1+0.287105{\alpha }^{0.843474})};\end{eqnarray} \tag{ 5 }$$

this expression is used in the numerical code that solves for the parcel's evolution, with L = L_∥ and $A=\pi {L}_{\perp }^{2}/4$.

Zoom In Zoom Out Reset image size

In order to determine all parameters in the model, there still needs to be a specification of the parcel's time-dependent shape and anisotropy. The cross-sectional area A is assumed to be identical to that of the large-scale background magnetic field, i.e., we assume that the parcel expands or contracts laterally to fill its local magnetic flux tube. As long as the radial magnetic field strength B_r is known, the principle of magnetic flux conservation provides the relative radial variation of $A\propto {B}_{r}^{-1}$. Thus, a given model of the solar wind has unique values for u_e (r), ρ_e (r), and A(r). The evolution in radial length of the parcel is specified with a scaling exponent σ as

$$\begin{eqnarray}&&L={L}_{0}{\left(\frac{A}{{A}_{0}}\right)}^{\sigma },\end{eqnarray} \tag{ 6 }$$

where quantities with subscript 0 are the initial conditions at t = 0. The exponent σ will be varied as a free parameter, and it is useful to note that there are several natural values that may pertain in different circumstances:

• 1.  
  A value of σ = − 1 would be appropriate if the volume of the parcel remains constant over time. This may occur for a parcel that retains a constant mass (M_i = constant) and is close to incompressible (ρ_i = constant).
• 2.  
  A value of σ = 0 constrains the length L to remain constant. This may be the case if the forces on the front and back faces of the parcel always remain roughly equal to one another. It may also be occurring in the outer corona for observed visible-light "flocculations" that appear to preserve their radial length scale from about 30 to 90 R_⊙ (see Figure 9 of DeForest et al. 2016).
• 3.  
  The specific value σ ≈ 0.14 would occur if the parcel's mass is conserved (M_i = constant) and its internal density remains proportional to the surrounding solar wind density (ρ_i ∝ ρ_e ). This value comes from the ZEPHYR coronal hole model discussed below, in which ρ_e ∝ r^−2.45 and A ∝ r^2.15 over the relevant range of heights to be compared with the observational data (7–12 R_⊙).
• 4.  
  Lastly, the value σ = 0.5 would occur if the aspect ratio α of the parcel remains constant over time. An initially spherical parcel may remain so if the local pressure gradient forces on it provide a roughly isotropic confinement.

It turns out that the models that successfully explain the observed COR2 parcel deceleration do not depend sensitively on the exponent σ. Thus, treating it as a free parameter probably does not limit the usefulness of these models. For another approach to simulating the evolution of parcel prolateness in the accelerating solar wind, see Pneuman & Cargill (1985).

For a specified background solar wind, the above model has two free parameters (C_D and σ) and five initial conditions (r_i,0, u_i,0, M_i,0, L₀, and A₀). The properties that are specified explicitly are r_i,0, u_i,0, L₀, and the initial value of the density ratio ρ_i /ρ_e . The initial parcel shape is assumed to be spherical, with ${A}_{0}=\pi {L}_{0}^{2}/4$, and M_i,0 is computed from the other quantities. For now, the parcel mass is presumed to remain constant over time, with M_i (t) = M_i,0.

Equations (1) and (2) will eventually be solved numerically, but it is also helpful to explore closed-form solutions when possible. Thus, we make the simplifying assumption that the wind speed u_e , background density ρ_e , and cross-sectional area A are all constants as a function of radial distance. This also constrains ρ_i and L to be constants as well. The assumption of spherical symmetry for the parcel gives the standard hydrodynamic value for the added mass coefficient, C_A = 1/2 (Lamb 1932).

Because essentially nothing in the simplified system depends on radial distance (other than the primary dependent variable u_i ), Equation (2) can be solved first for u_i (t), and then the radial trajectory r_i (t) can be determined subsequently. It is convenient to define two new dependent variables

$$\begin{eqnarray}&&U=| {u}_{i}-{u}_{e}| \,{\rm{and}}\,Q=\displaystyle \frac{{\rho }_{i}}{{\rho }_{e}}+\displaystyle \frac{1}{2},\end{eqnarray} \tag{ 7 }$$

where the two terms in Q correspond to the M_i and M_add terms in Equation (2). For now, the assumption is that Q remains fixed at its initial value Q₀. Equation (2) is thus rewritten as

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=-\displaystyle \frac{{C}_{{\rm{D}}}\,{U}^{2}}{{Q}_{0}L}.\end{eqnarray} \tag{ 8 }$$

If all quantities on the right-hand side are positive, then U will get smaller as time increases—i.e., asymptotically approaching zero as the parcel's motion becomes entrained into the surrounding solar wind. However, if the initial condition obeys u_i < u_e , then the time evolution of u_i can go from negative (inward flow) to positive (entrained outward flow). For a known initial condition U₀, this equation has an explicit solution

$$\begin{eqnarray}&&U(t)={U}_{0}{\left(1+{\omega }_{0}t\right)}^{-1}\end{eqnarray} \tag{ 9 }$$

with

$$\begin{eqnarray}&&{\omega }_{0}=\displaystyle \frac{{C}_{{\rm{D}}}\,{U}_{0}}{{Q}_{0}L}.\end{eqnarray} \tag{ 10 }$$

Note that significant deceleration occurs with a characteristic timescale of t ∼ 1/ω₀. The radial position of the parcel is given by integrating Equation (1),

$$\begin{eqnarray}&&{r}_{i}(t)={r}_{i,0}+{u}_{e}t-\frac{{U}_{0}}{{\omega }_{0}}\mathrm{ln}(1+{\omega }_{0}t),\end{eqnarray} \tag{ 11 }$$

where a sign convention was adopted that assumes u_i < u_e for all t ≥ 0.

Figure 3 shows a range of illustrative solutions of Equations (9) and (11). The initial conditions were chosen to agree with the observations of DeForest et al. (2014) for the largest height at which inflowing parcels were detected: r_i,0 = 12.4 R_⊙ and u_i,0 = −83.6 km s⁻¹. Here we also chose a fast solar wind speed of u_e = 500 km s⁻¹, typical for coronal holes at this distance.

Zoom In Zoom Out Reset image size

The only other parameter required to implement the above solutions is ω₀, and we can estimate likely values by examining the components of Equation (10). From the CME literature, a realistic range of values for C_D seems to be 0.001–1. The modified overdensity ratio Q₀ may be as small as 0.5 (if ρ_i < ρ_e ) and possibly as high as 100. The observations of DeForest et al. (2014) found that the parcel length scale L is probably no larger than about 0.5 R_⊙. Although it is unclear what the smallest value of L could be, the fact that these features are observable with STEREO/COR2 probably means that they are not much smaller than 1 pixel in size. Thus, we adopt a single-pixel lower limit for L of about 15'', or 0.016 R_⊙ (Howard et al. 2008). Combining the ranges for C_D, Q₀, and L with the initial speeds discussed above (i.e., U₀ = 583.6 km s⁻¹) gives a possible spread of about eight orders of magnitude for ω₀. Figure 3 shows the dynamical solutions for values between 3 × 10⁻⁸ s⁻¹ (green curve) and 1 s⁻¹ (black curve).

The green curves in Figure 3 correspond to very low values of ω₀. In this case, the drag force on the parcels is extremely weak; they just continue to flow at their initial downward velocity and eventually reach the solar surface. The violet and black curves in Figure 3 correspond to extremely large values of ω₀. In that case, the parcels become entrained very quickly into the ambient solar wind and flow out with speeds approaching u_e = 500 km s⁻¹. Note that all of the entrained solutions pass through the initial radius r_i,0 for a second time at the same velocity, u_i = 68.4 km s⁻¹. It is possible to write an analytic solution for this velocity, but it requires the use of the Lambert W function (Corless et al. 1996), with

$$\begin{eqnarray}&&{u}_{i}={u}_{e}\left\{1+{\left[{W}_{-1}\left(-\displaystyle \frac{{u}_{e}}{{U}_{0}}{e}^{-{u}_{e}/{U}_{0}}\right)\right]}^{-1}\right\},\end{eqnarray} \tag{ 12 }$$

and these solutions all exhibit ω₀ t = 0.35207.

Figure 3 also shows a fit to the "inflow ridge" (yellow region) that indicates the observed trend of parcel deceleration (DeForest et al. 2014) and is given by

$$\begin{eqnarray}&&{u}_{i,\mathrm{fit}}=118.39-46.739{r}_{i}+5.5295{r}_{i}^{2}-0.24790{r}_{i}^{3},\end{eqnarray} \tag{ 13 }$$

where r_i is in units of R_⊙ and u_i,fit is in units of km s⁻¹. This fitting formula is valid only between 7 and 12.4 R_⊙ and should not be applied outside this range. Note, however, that no choice of ω₀ matches the data well. Like the linear MHD-wave modes studied by Tenerani et al. (2016), the model curves have the opposite concavity to the observed trend.

It is possible that the analytic model described in Section 2.2 fails to match the data because its assumptions about constant values for the model parameters (e.g., u_e , ρ_e , A, L) were too simplistic. Thus, more realistic radial dependences for these plasma properties can be adopted, and the original drag equations can be solved numerically. Figure 4 shows coronal hole parameters from ZEPHYR, a one-fluid model of turbulent heating and acceleration along open magnetic flux tubes (Cranmer et al. 2007). Along with u_e (r), A(r), and ρ_e (r), we also show the radial Alfvén speed

$$\begin{eqnarray}&&{V}_{{\rm{A}}r}=\displaystyle \frac{{B}_{r}}{\sqrt{4\pi {\rho }_{e}}}\end{eqnarray} \tag{ 14 }$$

in the vicinity of relevant radii for the COR2 parcels. For this model, the Alfvén radius is located at r_A = 10.73 R_⊙, which is notably below the initial-condition parcel radius of 12.4 R_⊙. This means that a parcel that flows back down to the Sun from 12.4 R_⊙ must have an initial speed that is locally faster than the Alfvén speed. Some implications of this situation are discussed in Section 5, but for now we will assume that parcels can flow back toward the Sun at all radii ≤ 12.4 R_⊙.

Zoom In Zoom Out Reset image size

A numerical code was developed to integrate Equations (1) and (2) forward in time using first-order Euler steps of size Δt = 10⁻⁴ R_⊙/∣u_i,0∣. At each step, the current value of r_i was used to interpolate the local values of u_e , ρ_e , and A from the ZEPHYR model grids and to compute the local values of L and ρ_i from the expressions given in Section 2.1. At least 10,000 trial models were constructed with randomized choices for four key parameters: (1) C_D was sampled from a uniform grid (in logarithm) between 10⁻⁶ and 10⁺¹, (2) the shape exponent σ was sampled from a uniform grid of values between −6 and +6, (3) the initial parcel length L₀ was sampled from a uniform grid (in logarithm) between 0.001 and 0.5 R_⊙, and (4) the initial parcel overdensity ratio ${\left({\rho }_{i}/{\rho }_{e}\right)}_{0}$ was sampled from a uniform grid (in logarithm) between 10⁻³ and 10⁺³.

For each numerical model, the values of u_i at three specified radial distances were saved, compared with the observed inflow speeds, and used to compute a reduced χ² goodness-of-fit parameter,

$$\begin{eqnarray}&&{\chi }^{2}=\displaystyle \frac{1}{N}\displaystyle \sum _{j=1}^{N}{\left[\displaystyle \frac{{u}_{i}({r}_{j})-{u}_{i,\mathrm{obs}}({r}_{j})}{\delta {u}_{i}}\right]}^{2},\end{eqnarray} \tag{ 15 }$$

where N = 3. The sum is taken over the observations at r_j = {7, 9, 11} R_⊙, with corresponding velocities u_i,obs = {−22.86, −35.08, −56.62} km s⁻¹. A somewhat arbitrary uncertainty width of δ u_i = 5 km s⁻¹ was adopted as well; see Figure 7 of DeForest et al. (2014). If a given trial model never made it down to 7 R_⊙ (i.e., if it was rapidly entrained into the outflowing solar wind like the high-ω₀ models in Figure 3), it was assigned an artificially large value of χ² as a flag to neglect that set of input parameters.

After randomly sampling the four-dimensional parameter space discussed above, we found no deceleration trajectories that ended up matching the observed trend very well. The "best" models (which exhibited a minimum χ² ≈ 6) were very similar in appearance to the analytic solutions shown in Figure 3 for values of $\mathrm{log}{\omega }_{0}$ between about −6 and −5.5. The model with the lowest χ² had a value of ω₀ ≈ 1.7 × 10⁻⁶ s⁻¹, and the parcel reached a minimum radius of about 6.5 R_⊙ before turning around and flowing out with the solar wind. Thus, even when taking into account the radial variation of background solar wind parameters (which the analytic model above did not do), it seems that the assumption of a constant-M_i parcel cannot explain the inflow pattern observed by DeForest et al. (2014).

In order to better match the observed deceleration trend, it is clear that ω₀ should be allowed to vary even more than it does in the numerical models discussed in this section. Figure 3 indicates that ω₀ should be rather large when the parcel is first released and then become smaller as the parcel reaches the lowest observed height. In other words, the magnitude of the drag acceleration du_i /dt must decrease over time.

This model now consists of three coupled differential equations that must be solved simultaneously: Equations (1) and (2), and a new equation that describes the rate of change of parcel mass. We generalize from earlier versions of this kind of snowplow (or bulldozer) equation, such as that of Tappin (2006), with

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{i}}{{dt}}={\rho }_{e}{{AC}}_{{\rm{S}}}| {u}_{i}-{u}_{e}| \end{eqnarray} \tag{ 16 }$$

(see also Howard et al. 2007; Feng et al. 2015; Takahashi & Shibata 2017). The constant C_S is introduced as an effective snowplow efficiency that describes what fraction of the plasma in front of the parcel is actually incorporated into it. Note that Tappin (2006) assumed C_S = 1, which describes the parcel gaining an incremental mass dM_i as it traverses a radial length dL = ∣u_i − u_e ∣dt in the rest frame of the fluid, thus sweeping out an incremental volume dV = AdL containing particles with density ρ_e . Other values of C_S < 1 were considered by Takahashi & Shibata (2017). The right-hand side of Equation (16) is always positive, so M_i grows monotonically over time. Thus, Equation (2) says that the magnitude of the acceleration du_i /dt decreases over time, at least in comparison to the constant-M_i case.

As in Section 2.2, it is possible to first explore an analytic solution of the system in question (for the limiting case of constant values of u_e , ρ_e , A, and L) before solving the full equations numerically. In this case there are two coupled differential equations to solve: Equations (2) and (16). The same substitution of variables from Equation (7) allows the two equations to be written as

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=-\displaystyle \frac{{C}_{{\rm{D}}}\,{U}^{2}}{{QL}}\,{\rm{and}}\,\displaystyle \frac{{dQ}}{{dt}}=\displaystyle \frac{{C}_{{\rm{S}}}\,U}{L}.\end{eqnarray} \tag{ 17 }$$

As before, U gets smaller as time increases, and now the parcel mass proxy variable Q grows larger as time increases. The two coupled first-order equations can be combined into a single second-order equation,

$$\begin{eqnarray}&&U\,\displaystyle \frac{{d}^{2}U}{{{dt}}^{2}}=\left(\displaystyle \frac{{C}_{{\rm{S}}}}{{C}_{{\rm{D}}}}+2\right){\left(\displaystyle \frac{{dU}}{{dt}}\right)}^{2},\end{eqnarray} \tag{ 18 }$$

which has a closed-form solution of the form

$$\begin{eqnarray}&&U(t)={U}_{0}{\left(1+\omega t\right)}^{-1/\gamma }\end{eqnarray} \tag{ 19 }$$

with

$$\begin{eqnarray}&&\gamma =\displaystyle \frac{{C}_{{\rm{S}}}}{{C}_{{\rm{D}}}}+1\,{\rm{and}}\,\omega =\displaystyle \frac{\gamma {C}_{{\rm{D}}}\,{U}_{0}}{{Q}_{0}L}.\end{eqnarray} \tag{ 20 }$$

The full solution is completed by

$$\begin{eqnarray}&&Q(t)={Q}_{0}{\left(1+\omega t\right)}^{(\gamma -1)/\gamma }.\end{eqnarray} \tag{ 21 }$$

Note that both initial conditions U₀ and Q₀ are positive, and γ ≥ 1. The radial position of the parcel is given by integrating Equation (1),

$$\begin{eqnarray}&&{r}_{i}(t)={r}_{i,0}+{u}_{e}t-\displaystyle \frac{{Q}_{0}L}{{C}_{{\rm{S}}}}\left[{\left(1+\omega t\right)}^{(\gamma -1)/\gamma }-1\right],\end{eqnarray} \tag{ 22 }$$

where the same sign convention as in Equation (11) was used. The solutions from Section 2.2 are recovered exactly in the limit of γ = 1, which is equivalent to C_S = 0 (see also Mahajan 2020).

Figure 5(a) shows a range of illustrative solutions. The choices for r_i,0, u_i,0, and u_e were the same as in Section 2.2. For simplicity, we also chose Q₀ = 1.5 (i.e., initial density equipartition between the interior and exterior of the parcel) and L = 0.1 R_⊙. In this case, it is possible to find combinations of parameters that match the DeForest et al. (2014) data very well. It was found that one specific value for the ratio of snowplow to drag coefficients (i.e., γ ≈ 26) produced excellent agreement with the data. The best-fitting (central) curve in Figure 5(a) is for a model with C_S = 0.03 and C_D = 0.0012. The other curves were computed with C_S held fixed and C_D varied up and down by factors of (up to) two. Lower values of C_D imply weaker drag, so an inwardly moving parcel just keeps falling in toward the Sun. Higher values of C_D imply stronger drag, so the parcels are eventually swept into the outward-flowing solar wind.

Zoom In Zoom Out Reset image size

For the best-fitting analytic model in Figure 5(a), the parcel flows in from 12.4 to 7 R_⊙ over a time of about 26.1 hr. However, this model exhibits a value of ω = 1.74 × 10⁻⁴ s⁻¹, or a characteristic deceleration timescale of only 1/ω ≈ 1.6 hr. This is not an inconsistency because the definition of ω involves only the initial conditions of the parcel. If a new inverse-timescale variable was defined using the time-variable quantities (i.e., $\omega ^{\prime} =\gamma {C}_{{\rm{D}}}U(t)/Q(t)L$), the combined effects of mass gain and deceleration reduce its value over time. By the time the parcel reaches 7 R_⊙, this quantity has decreased to $\omega ^{\prime} \approx {10}^{-5}$ s⁻¹, driven mainly by an increase in Q(t) from 1.5 to 23. The corresponding timescale $1/\omega ^{\prime}$ is thus 27.7 hr, comparable to the total transit time.

The numerical integration code described in Section 2.3 was extended to solve Equations (1), (2), and (16). These models have radially varying trends in ρ_e , u_e , A, and L, given by the ZEPHYR polar coronal hole model shown in Figure 4. Another large set of randomized trial models was constructed for the parcel motion and mass gain, this time with different seed constants for the pseudorandom number generators than the ones used above. The parameters σ, L₀, and ${\left({\rho }_{i}/{\rho }_{e}\right)}_{0}$ were sampled as described in Section 2.3, and in this case both coefficients C_D and C_S were sampled randomly from a uniform grid (in logarithm) between 10⁻⁶ and 10⁺¹.

For this set of models with parcel mass gain, there were many solutions that ended up matching the observed deceleration trend of DeForest et al. (2014). More than 10,000 total models were generated, and the code was run until a sample size of 500 models with χ² ≤ 1 was accumulated. The smallest value of χ² in that sample was 0.059. Figure 5(b) shows nine sample trajectories of solutions with χ² ≤ 1. Note the large variety in the possible trajectories—from parcels that turn around just below the minimum observed height of 7 R_⊙ to ones that keep flowing monotonically down to the solar surface—that are all reasonably consistent with the observed velocities.

Figure 6 shows histograms of the input parameters for the set of 500 solutions with χ² ≤ 1. The ranges of good-fitting values for γ (which is essentially the ratio of snowplow to drag coefficients) and ω (the inverse deceleration timescale) were found to be rather narrow, with median values of γ = 27.3 and ω = 6.62 × 10⁻⁶ s⁻¹. The former value is almost identical to the best-fitting value of γ found for the analytic models described in Section 3.2. For the numerical models, the drag coefficient C_D also seemed to prefer intermediate values (with a median of C_D = 5.56 × 10⁻⁴), though the distribution of values matching the observations was broader than in the cases of γ and ω. The distribution of input values for C_S is not shown, but it had a median value of 0.0135. The combination of (more tightly constrained) medians for γ and C_D points to a slightly higher most likely value of C_S ≈ 0.0146.

Zoom In Zoom Out Reset image size

It is evident from Figure 6 that there do not seem to be any preferred values for σ, L₀, or ρ_i /ρ_e . Models that agree with the observed pattern of deceleration appear to be possible for nearly any of these input parameters. However, if we apply an independent estimate for a likely range of values for L₀ (say, 0.01–0.1 R_⊙), it is then possible to use the medians for C_D and ω to obtain a similarly likely range of values for the overdensity parameter Q₀ between about 20 and 200.

Additional information about how much snowplowing occurs in these models can be found in Figure 7. Figure 7(a) shows the monotonic gains in parcel mass M_i for a random subset of 50 models with χ² ≤ 1. At the lowest COR2 height of 7 R_⊙, these models exhibit a large range of final-to-initial mass ratios, i.e., between 5 and 7000, with a median ratio of about 30. Note that for the larger set of "bad" models (i.e., with χ² > 1; not shown) in which the parcels made it down to 7 R_⊙, the distribution of mass-gain ratios was much broader: sometimes as low as 1.0001, sometimes as high as 10⁵. Figure 7(b) shows how ρ_i /ρ_e evolves for the same subset of models from Figure 7(a). Note that the initial range of 10⁻³–10⁺³ expands gradually as the parcels evolve from 12.4 to 7 R_⊙. Due to the broad range of randomly chosen σ exponents (which controls the evolution of parcel volume), there is no strong correlation between the mass-gain ratio and the local overdensity ratio.

Zoom In Zoom Out Reset image size

It remains to be determined why the best-fitting values of both C_D and C_S fall so far below the order-unity expectations for these dimensionless coefficients. For the drag coefficient, it is possible that its small value can be understood as an effect of magnetic tension. Cargill et al. (1996) found that a horizontal flux rope, with an axial field pointing along the x-axis and rising along the z-axis, experiences strong turbulent drag with C_D ≈ 1. However, if the z-component of the field becomes strong enough to exceed a significant fraction of the flux rope's axial field, the vertical magnetic tension suppresses turbulent eddy-like vortices from forming behind the object and reduces C_D to values near zero. Conceptually, this may be similar to the reduction in drag when a blunt solid object (such as a sphere) is streamlined by lengthening and tapering its surface in the downwind direction. In both cases, trailing streamlines become stretched out more along the direction of motion, and the volume of the turbulent wake behind the object is reduced drastically. For high Reynolds number hydrodynamic flows, it is common to see reductions in C_D from values of order 1 to values of order 0.001 (Hoerner 1965; Munson et al. 2009). Thus, a similar reduction in drag may occur for the parcels considered here. Although we have no firm constraints on their internal magnetic geometry, they do seem to be moving parallel to a strong background field as in the high-B_z cases of Cargill et al. (1996).

Because the mass-gain mechanism assumed in Section 3 was just a local snowplow, the rate of gain was straightforwardly proportional to the relative velocity u_i − u_e . However, for ablative mass loss, it is usually assumed that the ram pressure (i.e., kinetic energy density) of the relative flow is the primary driver of the interaction. Thus, the mass-loss rate is proportional to ${\left({u}_{i}-{u}_{e}\right)}^{3}$. We use a similar formula to that given by Bronshten (1983),

$$\begin{eqnarray}&&\displaystyle \frac{{{dM}}_{i}}{{dt}}=-{\rho }_{e}A\,\displaystyle \frac{| {u}_{i}-{u}_{e}{| }^{3}}{{V}_{\mathrm{eff}}^{2}}\,{ \mathcal H }({M}_{i}-{M}_{e}),\end{eqnarray} \tag{ 23 }$$

where V_eff is a velocity-like quantity that describes the efficiency of the ablation. Note the presence of a Heaviside step function ${ \mathcal H }({M}_{i}-{M}_{e})$ on the right-hand side. This constrains the ablative mass loss to occur only when M_i > M_e . The process is assumed to shut off completely when M_i < M_e , i.e., when the parcel's inertia is totally dominated by that of the surrounding solar wind.

Because there does not seem to be any well-developed theory for the ablation of plasma parcels in the solar wind, we treat the ablation efficiency V_eff as a constant free parameter. We also speculate about a variety of approaches to estimating its value:

• 1.  
  For solid meteors entering a planetary atmosphere, ${V}_{\mathrm{eff}}^{2}$ is specified as a latent heat of ablation (with its velocity-squared units more typically reported as energy per unit mass). For rocky and metallic meteoric material, V_eff ≈ 1–3 km s⁻¹ (Baldwin & Sheaffer 1971; Chyba et al. 1993). These particular values are not likely to be relevant for plasma parcels in the solar wind, but we find it interesting to note the parallels with this other application.
• 2.  
  For astrophysical ram pressure stripping—which efficiently removes mass from regions with gas pressure less than the large-scale flow's ram pressure—it may be appropriate to assume that V_eff ≈ c_s , the adiabatic sound speed of the gas. For the solar corona, this hints at values around 150 km s⁻¹.
• 3.  
  Another way to estimate V_eff is to examine hydrodynamic models of large interstellar clouds that are destroyed when they encounter shocks from nearby supernovae (Cowie & McKee 1977; Klein et al. 1994). Klein et al. (1994) estimated that an initially spherical parcel would be mostly destroyed over a cloud-crushing timescale

  $$\begin{eqnarray}&&{\tau }_{{cc}}\approx \sqrt{\displaystyle \frac{{\rho }_{i}}{{\rho }_{e}}}\displaystyle \frac{L}{| {u}_{i}-{u}_{e}| },\end{eqnarray} \tag{ 24 }$$

  where we have converted the original expressions into the notation of this paper. However, another way to estimate the timescale over which a parcel undergoes substantial mass loss would be

  $$\begin{eqnarray}&&{\tau }_{m{\ell }}=\displaystyle \frac{{M}_{i}}{| {{dM}}_{i}/{dt}| }=\left(\displaystyle \frac{{\rho }_{i}}{{\rho }_{e}}\right)\displaystyle \frac{{{LV}}_{\mathrm{eff}}^{2}}{| {u}_{i}-{u}_{e}{| }^{3}}\end{eqnarray} \tag{ 25 }$$

  using Equation (23). Setting τ_cc = τ_{m ℓ} , we can solve for the ablation efficiency parameter

  $$\begin{eqnarray}&&{V}_{\mathrm{eff}}\approx | {u}_{i}-{u}_{e}| {\left(\displaystyle \frac{{\rho }_{i}}{{\rho }_{e}}\right)}^{-1/4},\end{eqnarray} \tag{ 26 }$$

  which is never far from just ∣u_i − u_e ∣ itself, i.e., a few hundred kilometers per second for our coronal parcels.

Following Sections 2.2 and 3.2, we describe a simplified system in which u_e , ρ_e , A, and L are assumed to be constants. Using the same dependent variables as above, there are two coupled differential equations,

$$\begin{eqnarray}&&\displaystyle \frac{{dU}}{{dt}}=-\displaystyle \frac{{C}_{{\rm{D}}}\,{U}^{2}}{{QL}}\,{\rm{and}}\,\displaystyle \frac{{dQ}}{{dt}}=-\displaystyle \frac{{U}^{3}}{{{LV}}_{\mathrm{eff}}^{2}}.\end{eqnarray} \tag{ 27 }$$

Inspired by the mathematical technique described by Bronshten (1983), we divide the dQ/dt equation by the dU/dt equation. This gives a first-order separable equation for Q as a function of U, which integrates to

$$\begin{eqnarray}&&Q={Q}_{0}\exp \left(\displaystyle \frac{{U}^{2}-{U}_{0}^{2}}{2{C}_{{\rm{D}}}{V}_{\mathrm{eff}}^{2}}\right).\end{eqnarray} \tag{ 28 }$$

As time increases, U decreases, so Q also decreases from its initial value of Q₀ to a final asymptotic value

$$\begin{eqnarray}&&{Q}_{\infty }={Q}_{0}\exp \left(\displaystyle \frac{-{U}_{0}^{2}}{2{C}_{{\rm{D}}}{V}_{\mathrm{eff}}^{2}}\right).\end{eqnarray} \tag{ 29 }$$

However, the above solution may not apply for all time. Under the present set of approximations, when Q first becomes equal to 3/2, it means that M_i = M_e and the Heaviside step function in Equation (23) becomes equal to zero. Thus, once Q reaches a value of 3/2, it remains fixed at this value for all future times. In that case, U(t) behaves like the mass-conserving model (Equation (9)) during this later phase of parcel evolution.

During the times when M_i > M_e and ablation occurs, Equation (28) can be inserted back into Equation (27) and integrated again. There does not appear to be a closed-form solution for U(t), but an inverted solution for t(U) can be written. If we define a new dimensionless variable ${y}^{2}={U}^{2}/(2{C}_{{\rm{D}}}{V}_{\mathrm{eff}}^{2})$, then

$$\begin{eqnarray}&&{\rm{\Omega }}t=\sqrt{\pi }\left[{\rm{erfi}}({y}_{0})-{\rm{erfi}}(y)\right]+\displaystyle \frac{{e}^{{y}^{2}}}{y}-\displaystyle \frac{{e}^{{y}_{0}^{2}}}{{y}_{0}},\end{eqnarray} \tag{ 30 }$$

where

$$\begin{eqnarray}&&{\rm{\Omega }}=\displaystyle \frac{{2}^{1/2}{C}_{{\rm{D}}}^{3/2}{V}_{\mathrm{eff}}}{{Q}_{\infty }L}\end{eqnarray} \tag{ 31 }$$

and the imaginary error function ${\rm{erfi}}(z)=-i\,{\rm{erf}}({iz})$ is real-valued for a real argument z. In the weak ablation limit of V_eff → ∞ , Equation (30) reduces to Equation (9) as it should.

As in Sections 2.3 and 3.3, a large set of numerical models was constructed to determine whether any combinations of parameters could reproduce the observed deceleration trend of DeForest et al. (2014). In this case, similar ranges of random values for C_D, σ, and L₀ were explored. The initial parcel density ratio ${\left({\rho }_{i}/{\rho }_{e}\right)}_{0}$ was sampled from a smaller range of 1.001–10⁺³, i.e., excluding values less than unity that would not exhibit any ablation. The efficiency parameter V_eff was sampled randomly from a uniform grid (in logarithm) over 12 orders of magnitude: from 10⁻⁶ c_s to 10⁺⁶ c_s , with a fiducial value of the coronal sound speed c_s = 150 km s⁻¹.

Unfortunately, there were no combinations of input parameters that agreed with the observed deceleration trend. The "best" models were ones with extremely weak ablation (i.e., the largest values of V_eff), and these were virtually identical to the mass-conserving models discussed in Sections 2.2–2.3. Minimum χ² values were similarly around 6, for minimum radii of about 6.5 R_⊙ and typical values of ω₀ ≈ 1.7 × 10⁻⁶ s⁻¹. Many of the ablated solutions ended up with overdensity ratios ρ_i /ρ_e ≲ 1 for r_i between 7 and 12 R_⊙. Parcels of this kind may not exhibit enough scattered-light contrast to ever be visible with instruments like COR2, so it is probably not surprising that no parcels with these kinematic properties were observed.