Evolving waves and turbulence in the outer corona and inner heliosphere: the accelerating expanding box

TL;DR

This paper introduces the accelerating expanding box (AEB) model, a numerical simulation tool that captures the evolution of solar wind turbulence from the corona to 1 AU, including effects of acceleration and inhomogeneity.

Contribution

The AEB model extends previous models by incorporating solar wind acceleration and background inhomogeneity, enabling more realistic simulations of solar wind turbulence evolution.

Findings

Successfully simulates Alfvén wave propagation across the Alfvén critical point.

Validates the model by reproducing expected wave behaviors in the solar wind.

Provides a framework for future detailed turbulence studies in the heliosphere.

Abstract

Alfv\'enic fluctuations in the solar wind display many properties reflecting an ongoing nonlinear cascade, e.g. a well-defined spectrum in frequency, together with some characteristics more commonly associated with the linear propagation of waves from the Sun, such as the variation of fluctuation amplitude with distance, dominated by solar wind expansion effects. Therefore both nonlinearities and expansion must be included simultaneously in any successful model of solar wind turbulence evolution. Because of the disparate spatial scales involved, direct numerical simulations of turbulence in the solar wind represent an arduous task, especially if one wants to go beyond the incompressible approximation. Indeed, most simulations neglect solar wind expansion effects entirely. Here we develop a numerical model to simulate turbulent fluctuations from the outer corona to 1 AU and beyond, including the sub-Alfv\'enic corona. The accelerating expanding box (AEB) extends the validity of previous expanding box models by taking into account both the acceleration of the solar wind and the inhomogeneity of background density and magnetic field. Our method incorporates a background accelerating wind within a magnetic field that naturally follows the Parker spiral evolution using a two-scale analysis in which the macroscopic spatial effect coupling fluctuations with background gradients becomes a time-dependent coupling term in a homogeneous box. In this paper we describe the AEB model in detail and discuss its main properties, illustrating its validity by studying Alfv\'en wave propagation across the Alfv\'en critical point.

Equations

$$\begin{array}{rl}& x=r\sin \theta \cos \phi ,\\ & y=r\sin \theta \sin \phi ,\\ & z=r\cos \theta ,\end{array}$$

$$\begin{array}{rl}\frac{\partial \mathbf{V}}{\partial t}+\mathbf{V}\cdot \nabla \mathbf{V}=& -\frac{1}{\rho }\nabla \left(p+\frac{1}{8\pi }{B}^2\right)\\ & -\frac{G{M}_{\odot }}{r^2}\hat{\mathbf{r}}+\frac{1}{4\pi \rho }\mathbf{B}\cdot \nabla \mathbf{B},\end{array}$$

$$\frac{\partial \mathbf{B}}{\partial t}+\mathbf{V}\cdot \nabla \mathbf{B}=-\mathbf{B}(\nabla \cdot \mathbf{V})+\mathbf{B}\cdot \nabla \mathbf{V},$$

$$\frac{\partial p}{\partial t}+\mathbf{V}\cdot \nabla p=-\gamma p(\nabla \cdot \mathbf{V}),$$

$$\frac{\partial \rho }{\partial t}+\nabla \cdot (\rho \mathbf{V})=0.$$

$$\mathbf{V}(R,\mathbf{r},t)=\mathbf{U}(R)+\mathbf{u}(R,\mathbf{r},t),$$

$$\mathbf{B}(R,\mathbf{r},t)={\mathbf{B}}_0(R)+\mathbf{b}(R,\mathbf{r},t)$$

$$p(R,\mathbf{r},t)=p_0+\delta p(R,\mathbf{r},t),$$

$$\rho (R,\mathbf{r},t)={\rho }_0(R)+\delta \rho (R,\mathbf{r},t).$$

$$\frac{d}{dr}=\frac{\partial }{\partial r}+\frac{\ell }{R_0}\frac{\partial }{\partial R}\equiv \frac{\partial }{\partial r}+\epsilon \frac{\partial }{\partial R}.$$

$$\begin{array}{rl}{\rho }_0& \mathbf{U}\cdot \nabla \mathbf{U}=-\nabla p_0-{\rho }_0\frac{M_{\odot }G}{r^2}\hat{\mathbf{r}}+\frac{1}{4\pi }\left(\nabla \times {\mathbf{B}}_0\right)\times {\mathbf{B}}_0\\ & +\nabla :\left(\frac{1}{4\pi }<\mathbf{bb}>-<\rho \mathbf{uu}>-\frac{1}{8\pi }<b^2>\mathcal{I}\right)\\ & -<\delta \rho \frac{\partial \mathbf{u}}{\partial t}>-<\delta \rho (\mathbf{u}\cdot \nabla \mathbf{U}+\mathbf{U}\cdot \nabla \mathbf{u})>,\end{array}$$

$$\begin{array}{rl}& \mathbf{U}\cdot \nabla {\mathbf{B}}_0=-{\mathbf{B}}_0(\nabla \cdot \mathbf{U})+{\mathbf{B}}_0\cdot \nabla \mathbf{U}\\ & -\nabla :\left(<\mathbf{ub}>-<\mathbf{bu}>\right),\end{array}$$

$$\mathbf{U}\cdot \nabla p_0=-\gamma p_0(\nabla \cdot \mathbf{U})-\gamma <\nabla \cdot (\delta p\mathbf{u})>,$$

$$\mathbf{U}\cdot \nabla {\rho }_0=-{\rho }_0(\nabla \cdot \mathbf{U})-<\nabla \cdot (\delta \rho \mathbf{u})>,$$

$$\tilde{\mathbf{b}}=\frac{\mathbf{b}}{\sqrt{4\pi {\rho }_0}},\tilde{p}=\frac{p}{c_s{\rho }_0},$$

$$\begin{array}{rl}\frac{\partial \mathbf{u}}{\partial t}+& \mathbf{U}\cdot \nabla \mathbf{u}+\mathbf{u}\cdot \nabla \mathbf{U}+\mathbf{u}\cdot \nabla \mathbf{u}-{\mathbf{V}}_{\mathbf{a}}\cdot \nabla \tilde{\mathbf{b}}\\ & -\tilde{\mathbf{b}}\cdot \nabla {\mathbf{V}}_{\mathbf{a}}-\tilde{\mathbf{b}}\cdot \nabla \tilde{\mathbf{b}}=\frac{1}{2}\tilde{\mathbf{b}}(V_{a,r}+{\tilde{b}}_r)(\ln {\rho }_0{)}^{\prime }\\ & +\frac{1}{2}{\mathbf{V}}_{\mathbf{a}}{\tilde{b}}_r(\ln {\rho }_0{)}^{\prime }-\nabla {\tilde{P}}_T\\ & -\frac{\delta }{1+\delta }\left[\tilde{\mathbf{b}}\cdot \nabla \tilde{\mathbf{b}}+\tilde{\mathbf{b}}\cdot \nabla {\mathbf{V}}_{\mathbf{a}}+{\mathbf{V}}_{\mathbf{a}}\cdot \nabla \tilde{\mathbf{b}}+{\mathbf{V}}_{\mathbf{a}}\cdot \nabla {\mathbf{V}}_{\mathbf{a}}\right]\\ & -\frac{\delta }{1+\delta }\left[\frac{1}{2}(\tilde{\mathbf{b}}+{\mathbf{V}}_{\mathbf{a}})(V_{a,r}+{\tilde{b}}_r)(\ln {\rho }_0{)}^{\prime }\right]-⟨\dots ⟩,\end{array}$$

$$\begin{array}{rl}\frac{\partial \tilde{\mathbf{b}}}{\partial t}+& \mathbf{U}\cdot \nabla \tilde{\mathbf{b}}+\mathbf{u}\cdot \nabla \tilde{\mathbf{b}}\\ & +\mathbf{u}\cdot \nabla {\mathbf{V}}_{\mathbf{a}}-\tilde{\mathbf{b}}\cdot \nabla \mathbf{U}-\tilde{\mathbf{b}}\cdot \nabla \mathbf{u}-{\mathbf{V}}_{\mathbf{a}}\cdot \nabla \mathbf{u}=\\ & -\frac{1}{2}\tilde{\mathbf{b}}(U+u_r)(\ln {\rho }_0{)}^{\prime }-\frac{1}{2}{\mathbf{V}}_{\mathbf{a}}{u}_r(\ln {\rho }_0{)}^{\prime }\\ & -\tilde{\mathbf{b}}(\nabla \cdot \mathbf{U})-\tilde{\mathbf{b}}(\nabla \cdot \mathbf{u})-{\mathbf{V}}_{\mathbf{a}}(\nabla \cdot \mathbf{u})-⟨\dots ⟩,\end{array}$$

$$\begin{array}{rl}\frac{\partial \delta \tilde{p}}{\partial t}& +\mathbf{U}\cdot \nabla \delta \tilde{p}+\mathbf{u}\cdot \nabla \delta \tilde{p}+u_r\frac{1}{\gamma }{c}_s^{\prime }\\ & +\frac{c_s}{\gamma }{u}_r(\ln {\rho }_0{c}_s{)}^{\prime }+\delta \tilde{p}(U+u_r)(\ln {\rho }_0{c}_s{)}^{\prime }=\\ & -(c_s+\gamma \delta \tilde{p})(\nabla \cdot \mathbf{u})-\gamma \delta \tilde{p}(\nabla \cdot \mathbf{U})-⟨\dots ⟩,\end{array}$$

$$\begin{array}{rl}& \frac{\partial }{\partial t}\left(\frac{\delta \rho }{{\rho }_0}\right)+\mathbf{U}\cdot \nabla \left(\frac{\delta \rho }{{\rho }_0}\right)+\mathbf{u}\cdot \nabla \left(\frac{\delta \rho }{{\rho }_0}\right)\\ & +U\frac{\delta \rho }{{\rho }_0}(\ln {\rho }_0{)}^{\prime }+u_r\left(\frac{\delta \rho }{{\rho }_0}+1\right)(\ln {\rho }_0{)}^{\prime }=-\left(\frac{\delta \rho }{{\rho }_0}\right)(\nabla \cdot \mathbf{U})\\ & -\left(\frac{\delta \rho }{{\rho }_0}+1\right)(\nabla \cdot \mathbf{u})-⟨\dots ⟩.\end{array}$$

$$\begin{array}{rl}\nabla {\tilde{P}}_T=\frac{1}{\rho }\nabla & \left[\delta \tilde{p}({\rho }_0{c}_s)+\frac{1}{2}{\rho }_0({\tilde{b}}^2+2\tilde{\mathbf{b}}\cdot {\mathbf{V}}_{\mathbf{a}})\right]\\ & -\frac{1}{{\rho }_0}\frac{\delta }{1+\delta }\nabla \left(p_0+\frac{1}{2}{V}_a^2\right).\end{array}$$

$$\frac{dR}{dt}=U(R).$$

$$\delta \mathbf{r}=(\delta x,y,z),\delta x=x-R(t).$$

$$U(r)=U(R)+U^{\prime }\delta x+\mathcal{O}({ϵ}^2,{\alpha }^2),$$

$$\mathbf{U}(r)=U(r)\left[\frac{\delta x+R}{r}\right]\hat{\mathbf{x}}+U(r)\frac{y}{r}\hat{\mathbf{y}}+U(r)\frac{z}{r}\hat{\mathbf{z}},$$

$$\mathbf{U}(R)=\left[U(R)+U^{\prime }\delta x\right]\hat{\mathbf{x}}+U(R)\frac{y}{R}\hat{\mathbf{y}}+U(R)\frac{z}{R}\hat{\mathbf{z}}.$$

$$\begin{array}{rl}& \bar{x}=[x-R(\bar{t})]\left(L_{0x}/L_x(t)\right),\\ & \bar{y}=y\left(R_0/R(t)\right),\\ & \bar{z}=z\left(R_0/R(t)\right),\\ & \bar{t}=t,\end{array}$$

$$\frac{dR}{d\bar{t}}=U,\frac{d{L}_x}{d\bar{t}}=L_x{U}^{\prime }.$$

$$\nabla =\frac{L_{0x}}{L_x}\frac{\partial }{\partial \bar{x}}\hat{\mathbf{x}}+\frac{R_0}{R(t)}\frac{\partial }{\partial \bar{y}}\hat{\mathbf{y}}+\frac{R_0}{R(t)}\frac{\partial }{\partial \bar{z}}\hat{\mathbf{z}}\equiv \bar{\nabla },$$

$$\epsilon \frac{\partial }{\partial R}=\epsilon \frac{1}{U}\frac{\partial }{\partial \tau }.$$

$$\begin{array}{r}\bar{\mathbf{U}}\left[U(R)+U^{\prime }\left(\frac{L_x}{L_{0x}}\right)\bar{x}\right]\hat{\mathbf{x}}+U(R)\frac{\bar{y}}{R_0}\hat{\mathbf{y}}+U(R)\frac{\bar{z}}{R_0}\hat{\mathbf{z}},\end{array}$$

$$\frac{\partial }{\partial t}=\frac{\partial }{\partial \bar{t}}-\bar{\mathbf{U}}\cdot \bar{\nabla }.$$

$$\nabla \cdot \mathbf{U}=\frac{2U(R)}{R(t)}+U^{\prime }(R),$$

$$\mathbf{A}\cdot \nabla \mathbf{U}(R)=U^{\prime }(R){\bar{A}}_x\hat{\mathbf{x}}+\frac{U(R)}{R(t)}{\bar{A}}_y\hat{\mathbf{y}}+\frac{U(R)}{R(t)}{\bar{A}}_z\hat{\mathbf{z}}.$$

$$\frac{1}{{\bar{\rho }}_0}\frac{d{\bar{\rho }}_0}{dR}=-\left(\frac{2}{R}+\frac{U^{\prime }}{U}\right)\Rightarrow {\bar{\rho }}_0={\rho }_{\ast }{\left(\frac{R_0}{R}\right)}^2\frac{U_{\ast }}{U},$$

Full text

\pacs

1. Introduction

In this paper we propose and discuss the accelerating expanding box (AEB). The AEB is a three-dimensional numerical model based on the compressible Magnetohydrodynamic (MHD) description of the plasma, that allows numerical studies of the evolution of turbulence and structures in the accelerating, expanding solar wind. The model presented here generalizes the expanding box approximation (Grappin et al., 1993) to regions close to the sun, overcoming the limitations of a previous version of the model introduced in Tenerani & Velli (2013) to guarantee the exact conservation of the proper non-WKB invariants for counter-propagating Alfvén waves in the limit of small amplitude (Heinemann & Olbert, 1980).

Understanding how the wind expansion affects plasma dynamics and the evolution of fluctuations is the starting point to understand the more general problem of solar wind heating and acceleration via turbulent dissipation and Reynolds stresses (wave-driven winds). In-situ measurements show that the solar wind is indeed in a turbulent state: while in the slow wind the magnetic energy spectrum displays in general a power-law index close to the standard $-5/3$ predicted by the Kolmogorov phenomenology, hot fast streams originating from coronal holes appear to be permeated by almost incompressible, large amplitude Alfvénic fluctuations (Coleman, 1966; Belcher & Davis, 1971) characterized by a “younger” (flatter) spectrum which evolves with heliocentric distance (Bavassano et al., 1982). Despite the wealth of observations dating back from the ’60s, and presently covering distances from 0.3 out to many AU over a wide range of heliographic latitudes (see for instance Horbury et al. (2005) and Bruno & Carbone (2013)), the mechanisms leading to this Alfvénic turbulent state and the possible role of such large amplitude fluctuations in the acceleration and heating, specifically of fast streams (Cranmer et al., 2015), remain one of the main questions concerning the solar wind. In addition, there is no present comprehensive understanding of the variation of turbulence characteristics with solar wind state (slow, fast, et.c.) throughout the solar activity cycle. In particular, a high degree of Alfvénicity has been found also in specific slow wind periods, most probably associated with streams originating close to the boundary of coronal holes (D’Amicis et al., 2016).

Alfvénic turbulence in fast streams displays velocity-magnetic field correlations corresponding to waves largely propagating away from the sun. The Alfvénic spectrum dominates the frequency range $f\simeq 10^{-4}-10^{-2}$ Hz and shows a low frequency part where the fluctuation energy density scales as $E(f)\sim f^{-1}$, and a high frequency one where $E(f)\sim f^{-5/3}$. Observations indicate that the spectrum evolves towards the same “standard” state found in the slow wind at increasing heliocentric distances, as the knee separating the steeper from the flatter parts of the spectrum gradually shifts towards the lower frequencies. This is evidence that a nonlinear energy cascade must be at play, as was originally suggested by Coleman (1968). At the same time the imbalance between forward and backward propagating perturbations (cross helicity) is observed to decrease with heliocentric distance (Roberts et al., 1987). This observation is at odds with theory and simulations of homogeneous, incompressible turbulence showing that an initial imbalance between counter-propagating waves would be reinforced by nonlinear interactions via a process called dynamical alignment (Dobrowolny et al., 1980).

In this respect, neither phenomenological theories nor direct numerical simulations of homogeneous turbulence seem adequate to understand the evolution and energy cascade of Alfvénic fluctuations in the solar wind. The major difficulty comes from the fact that the solar wind is an inhomogeneous medium, complicating both theoretical and numerical approaches to the problem. The solar wind inhomogeneity comprises many different aspects.

First and foremost the overall expansion of the underlying solar atmosphere affects the evolution of fluctuations through the gradients in the average fields that introduce reflection and coupling among different wave modes, especially at frequencies low with respect to the expansion rate (Heinemann & Olbert, 1980; Lou, 1993a, b, c; Zhou & Matthaeus, 1989; Velli et al., 1989, 1991; Matthaeus et al., 1999a). These couplings exist in the equations quite independently of any normal mode analysis and have an important effect on the energy density in the fluctuations, causing a slow radial evolution. This, in linear theory and in the absence of coupling to compressible modes and transverse gradients, follows the conservation law of wave action (Bretherton & Garrett, 1969; Heinemann & Olbert, 1980). Not so surprisingly an overall radial decrease of the root-mean-square energies $<\delta b^{2}>$ is observed in the low frequency part of the Alfvénic spectrum, in good agreement with the theory of wave propagation in the expanding solar wind that, for frequencies larger than the expansion rate, predicts $<\delta b^{2}>\sim r^{-3}$ (Bavassano et al., 1982; Roberts et al., 1990). Nevertheless, this approach based on noninteracting waves fails in accounting for many other features of the solar wind, such as the observed radial decrease of cross helicity at high frequencies (Velli, 1993) or the evolution of the Alfvén ratio (ratio between the kinetic and magnetic energy of the fluctuations) with heliocentric distance.

Other aspects of solar wind inhomogeneity might therefore play an important role. Among these are the transverse gradients (shear) between fast and slow solar wind streams, a gradient which also steepens with distance from the Sun evolving into corotating interaction regions. The interaction of fluctuations with such shears was proposed to provide a mechanism for driving the turbulent cascade, even when the solar wind is initially dominated by outwards-type fluctuations (Roberts et al., 1990) and the stream shear interface is stable. Also, results from incompressible nonlinear simulations (Roberts et al., 1992) and transport turbulence models (Zank et al., 1996; Matthaeus et al., 1999b; Breech et al., 2005) are in good agreement with observations in the outer heliosphere pointing to the fact that the regions of strong shears found in the ecliptic plane display fluctuations evolving from an almost purely Alfvénic state towards a state with small cross-helicity (“mixed” state) more rapidly than in the polar regions, where the flow shear is weaker. Note that (Grappin & Velli, 1996) also pointed to the importance of such effects in their expanding box simulations with stream shears. However, as shown in  (Velli et al., 1990; Velli, 1993) the analytical closure models face important difficulties when approaching the Alfvén point.

It becomes clear that any study of wave and turbulence evolution in the solar wind needs to include nonlinearities, expansion effects and stream structures. Unfortunately, fully nonlinear eulerian simulations spanning a wide range of spatial scales which could include gradients of the background all the way down to the smallest scales reached through the turbulent cascade (say, from one AU all the way down to $10^{-4}-{10^{-6}}$ AU), are still prohibitively expensive in terms of numerical costs, especially if one wants to go beyond the incompressible approximation.

Turbulence driven by Alfvén wave reflection triggered by the strong inhomogeneities in the lower corona has been investigated numerically, although in most cases using simplified nonlinear models to describe the fluctuations. Shell model calculations, which approximate incompressible nonlinear interactions in directions orthogonal to the radial to form a set of coupled first order nonlinear PDEs for each perpendicular scale, have shed some light on the evolution of Alfvén waves propagating away from the Sun (Verdini et al. (2012) and references therein). For example, Verdini et al. (2012) found that the double power-law spectrum, $-1$ and $-5/3$, measured in situ might arise from the combination of strong and weak cascades respectively for reflected (inward propagating) and outward modes at the base of the solar corona. Semi-phenomenological models which adopt ad-hoc terms for turbulent spectral transfer and heating (Hossain et al., 1995) have also been useful to investigate the role of reflection-driven turbulence in coronal heating and solar wind acceleration (Matthaeus et al., 1999a; Dmitruk et al., 2001; Verdini et al., 2010).

Recently, direct three-dimensional simulations of turbulence in a given solar wind flux tube extending up to the Alfvén critical point (i.e. where the solar wind speed equals the Alfvén speed) have been carried out by Perez & Chandran (2013), yet within the reduced MHD framework. These simulations showed that the anomalous inwards Alfvén modes induced by the inhomogeneity trigger a nonlinear cascade in the sub-Alfvénic region via wave-reflection-driven turbulence (Velli et al., 1989), eventually developing a $k_{\bot}^{-1}$ spectrum in wave-number space. The simplifying assumption of incompressibility, an assumption commonly accepted in the studies of Alfvénic turbulence, may be justified on the basis of observational evidence: incompressible fluctuations seem to be the only ones to survive at large heliocentric distances, at least in the fast streams. Slow modes, which tend naturally to steepen into shocks, are expected to be rapidly dissipated in the chromosphere and transition region, while oblique fast modes should suffer strong refraction because of the large gradients in those regions. Nevertheless, compressible effects may arise as the outcome of transient processes which should not be neglected, an interesting example being the parametric decay of large amplitude Alfvén waves that may take place in the outer corona and inner heliosphere (Malara & Velli, 1996; Del Zanna et al., 2001; Tenerani & Velli, 2013) and nonlinear mode conversion from Alfvén to slow mode that may contribute to coronal heating and solar wind acceleration via shock formation and dissipation (Suzuki & Inutsuka, 2005; Matsumoto & Suzuki, 2012).

The expanding box model (EBM) was introduced a few years ago for the purpose of taking into account in a simple way the basic effects of the expansion on a compressible MHD plasma, in the limit of a radially uniform and super-Alfvénic solar wind (Velli et al., 1992). The expanding box approximation is based on a separation of spatial scales which allows for a simultaneous description of the large scale evolution due to the expansion of the solar wind and of the more rapid dynamics at smaller scales via a semi-lagrangian approach. Expansion effects can be taken into account in the traditional EBM easily, by rescaling variables according to the radial expansion and by including first order effects of the solar wind velocity gradients, provided the dimensions of the box are small with respect to the average heliocentric distance so that a local cartesian geometry is still appropriate (Grappin et al., 1993; Del Zanna et al., 2015). The EBM has proved a useful tool, and it has been extensively employed to investigate MHD instabilities and turbulence (Grappin et al., 1993; Rappazzo, 2005; Dong et al., 2014; Verdini & Grappin, 2016), even within stream interaction regions (Grappin & Velli, 1996), as well as the onset and development of kinetic instabilities driven by solar wind expansion with the hybrid expanding box (Liewer, 2001; Matteini et al., 2006; Hellinger and Travnicek, 2013, 2015; Maneva et al., 2015).

The AEB extends and improves the previous EBM to the regions close to the Sun, at a few solar radii, where it becomes necessary not only to include the radial acceleration of the solar wind speed directly, but also to extend the validity of the EBM close to the sonic and Alfvén critical points by taking into account explicitly the inhomogeneity of the background density and magnetic field. This model and its derivation therefore explicitly rely on a multiple scale approach that is similar in spirit to turbulence transport theories (Marsch & Tu, 1989; Zhou & Matthaeus, 1990; Zank et al., 1996; Breech et al., 2005, 2008; Zank et al., 2017). However, while the latter are designed to describe the radial evolution of a limited number of moments of the fluctuations (e.g. energy and cross helicity) and require closure assumptions, the resulting nonlinear equations for the fluctuations in the AEB can be solved via direct numerical integration in the frame co-moving with the solar wind. In this sense our model differs from the one adopted by Nariyuki (2015), who uses an accelerating expanding box for studying nonlinear Alfvén wave evolution in the Derivative Nonlinear Schrödinger approximation: in analogy with the traditional EBM, Nariyuki (2015) derives a set of MHD equations that includes also the effects of the radial acceleration of the wind. However, this description does not include the expansion effects that arise from the gradients of background fields other than the wind speed — in particular gradients of the Alfvén velocity which are non negligible in the acceleration region; second, contrary to the present model, Nariyuki (2015) considers a uniformly accelerating wind.

The inclusion of the acceleration region into the model, and in particular of the sub-Alfvénic corona, is important because in that region the strongest wave reflections occur and fluctuations’ amplitudes are expected to reach their maximum – both central elements for the development of a turbulent cascade (Velli et al., 1989; Matthaeus et al., 1999a; Verdini et al., 2009; Perez & Chandran, 2013). Simulations of evolving fluctuations from the acceleration region out to the inner heliosphere is of further interest in view of the upcoming missions Solar Probe Plus (Fox et al., 2015) and Solar Orbiter (Müller et al., 2013). In particular, Solar Probe Plus will arrive closer to the Sun than any previous spacecraft by crossing the Alfvén critical point, placed somewhere between $10-20\,R_{s}$ (DeForest et al., 2014), providing for the first time in-situ evidence of the state of the turbulence in the sub-Alfvénic corona.

The paper is organized as follows: in Section 2 we introduce the basic idea of the (accelerating) expanding box; in Section 3 we write the equations describing the average solar wind that represents our resulting solar wind model; in Section 4 we derive the AEB governing equations for fluctuations and in Section 5 we write down the governing equations of the AEB in terms of generalized Elsässer variables used for numerical integration; in Section 6 we discuss the range of validity of the AEB and a numerical test of Alfvén wave propagation, while the final discussion and summary is deferred to Section 7.

Here we lay emphasis mostly on propagation of MHD waves to validate the model. Nevertheless, the latter applies to a wider class of problems involving some localized structure, such as corotating stream interaction regions at sector boundaries of fast-slow wind (Grappin & Velli, 1996) and current sheets in the streamer belt (Rappazzo, 2005). Finally, the model also provides a general approach to describe the dynamics of a given parent system of equations in an expanding, accelerating medium.

2. The Accelerating Expanding Box

Start from MHD with a polytropic energy equation of index $\gamma$, within the absolute frame of reference $\{x,y,z\}$ with origin in the Sun defined by

[TABLE]

$r$, $\theta$, and $\varphi$ being the radial distance from the Sun, the latitude, and the longitude of a given point of the plasma, respectively. In this frame, the set of equations for velocity, magnetic field, pressure and density is

[TABLE]

From eqs. (2d) it is possible to obtain, through Reynolds average, two coupled set of equations describing respectively fluctuations and a spherically symmetric, statistically steady background (see e.g. Whang (1980); Zhou & Matthaeus (1990)). Within this standard framework, we employ a two-scale analysis by assuming that the average fields, labelled with a “zero”, are functions of a slow radial variable $R=(\ell/R_{0})\,r$, with $R_{0}$ a reference heliocentric distance and $\ell$ a dynamical fast scale length such that $\ell/R_{0}\sim\varepsilon\ll 1$. The parameter $\varepsilon$ will be de facto our power expansion parameter. Fluctuations, on the other hand, are in general functions of $R$, via couplings with the background, of a fast spatial variable ${\bf r}$ and of time $t$. Fields may therefore be separated into an average part, steady in time and slowly varying in space, and a fluctuating part as follows,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Note that introduction of $R$ is a useful artifice, but $R$ and $r$ are not two independent variables (see e.g. Bender & Orszag (1978)), so that the radial derivative can be inferred from the usual chain rule for differentiation,

[TABLE]

Hereafter it will be implied that gradients along the radial direction are given by both the small and the large scale variations as expressed by eq. (7).

Reynolds average over eqs. (2d) leads to the following solar wind equation,

[TABLE]

where brackets denote time-average, $\mathcal{I}$ is the unit dyadic, $M_{\odot}$ the solar mass and $G$ the gravitational constant. Equation (8) is complemented by the equations for the background magnetic field, pressure and density,

[TABLE]

and by the governing equations for the fluctuations, obtained by subtracting eqs. (8)–(9c) from eqs. (2d).

The equations for the fluctuations are best written in terms of the normalized magnetic field and pressure

[TABLE]

where $c_{s}=\sqrt{\gamma p_{0}/\rho_{0}}$ is the sound speed. We write them below by denoting with a “prime” the large scale radial derivative $\partial/\partial R$, and where $\delta=\delta\rho/\rho_{0}$ and ${\bf V_{a}}={\bf B_{0}}/\sqrt{4\pi\rho_{0}}$ is the Alfvén velocity,

[TABLE]

Again, brackets $\langle\dots\rangle$ stand for the nonlinear feedback on the average fields that must be subtracted, whereas the gradient of the total pressure $\tilde{P}_{T}$ in eq. (11a) is given by

[TABLE]

The aim of the AEB is to describe how inhomogeneity and expansion affect the radial evolution of fluctuations given in eqs. (11d). For this reason, the AEB does not take into account the nonlinear feedback of fluctuations on the average quantities in eqs. (9c), and the profile of ${\bf U}$ is given a priori in this model. The main idea is to change to a new frame of reference, the co-moving frame, that allows to follow in a semi-lagrangian fashion the evolution of a small angular sector of plasma (the simulation box) which is in accelerating, spherical expansion because of the drag of the average outflow ${\bf U}=U(R){\bf\hat{r}}$. A schematic representation of the AEB is shown in Fig. 1. In this model the expansion is taken into account by introducing an explicit dependence on the variable $R(t)$, which is defined here as the average heliocentric distance of the plasma sector, taken along the $x$ axis. The average distance $R(t)$ is parameterized in time and can be considered the lagrangian variable of the background, obeying to the following equation,

[TABLE]

At any given time, a point of the plasma parcel has coordinates defined by $R(t)$ and by the local position vector ${\delta\bf{r}}$,

[TABLE]

As can be seen from eqs. (9c) and (11d), effects of the solar wind expansion on the temporal evolution of scalar and vectorial fields enter explicitly through the convective derivative associated with the flow ${\bf U\cdot\boldsymbol{\nabla}}$, its divergence ${\bf\boldsymbol{\nabla}\cdot U}$, and, for a generic vector ${\bf A}$, the material derivative ${\bf A\cdot\boldsymbol{\nabla}U}$. If the plasma volume is small, i.e., if the length of the box is of order $\delta x/R\sim\epsilon\ll 1$ and the transverse dimensions are of order $y/R,\,z/R=\alpha\ll 1$, at any given time one can expand the spherical coordinates around $R$ to first order in the small parameters $\alpha$ and $\epsilon$, so as to obtain an explicit dependence on $R$. This procedure allows us to include the expansion effects on the local dynamics. The difference with the EBM consists in retaining the first correction due to the acceleration of $U(r)$, that gives an additional contribution to the expansion rates, as we will see.

In practice, if the radial inhomogeneity of the solar wind speed $U(r)$ is retained up to first order derivative, $U(r)$ can be written locally as

[TABLE]

where $U^{\prime}=({dU}/{dr})|_{R}$ has at most the order of magnitude $U^{\prime}\sim U/R$. The components of ${\bf U}$ are given by

[TABLE]

which, together with eq. (15), yields to first order in $\alpha,\epsilon$:

[TABLE]

According to the spherical expansion expressed by eq. (17), the plasma parcel undergoes a stretching in the transverse directions $y$ and $z$ at a rate $U/R$ because of the advection of the solar wind velocity ${\bf U}$. On the other hand, while considering a finite radial acceleration within the box, there is also an additional stretching along the parallel (radial) direction. Since we are assuming that the length of the box $L_{x}$ is small, $L_{x}/R\ll 1$, the latter increases at a rate $dL_{x}/dt=L_{x}U^{\prime}$.

The co-moving frame is thus defined by the following change of variables:

[TABLE]

[TABLE]

In the equations above, $R_{0}$ and $L_{0x}$ are the initial average heliocentric distance and parallel length of the box, respectively. It is easily verified that fast spatial derivatives transform as

[TABLE]

whereas the slow spatial variation along the radial direction, that we have introduced in eq. (7), becomes a slow time variation as follows

[TABLE]

We recall however that in reality fields are functions of $t$ alone, and $\tau$ has to be ultimately included within the variable $t$.

The solar wind velocity can be written, in the new frame,

[TABLE]

and consequently the time derivative is given by

[TABLE]

The time-dependent rescaling of coordinates defined in eqs. (18) simplifies the convective derivative associated with ${\bf U}$, while the slow radial gradient implicitly becomes a slow time variation, $U(\partial/\partial R)=\partial/\partial\tau$. In this way, any dependence on the local coordinates due to the advection is eliminated, and we can assume a background locally homogeneous, yet slowly varying in time, which allows for the use of periodic boundary conditions. Finally, the divergence of ${\bf U}$ and the material derivative of a generic vector ${\bf A}$ can be written in the co-moving frame, according to our ordering assumptions, in the following way:

[TABLE]

[TABLE]

In summary, the expansion effects in the co-moving frame $\{\bar{x},\,\bar{y},\,\bar{z}\}$ are represented by: $(i)$ the rescaling of spatial derivatives, described by eqs. (18); $(ii)$ the divergence of the solar wind velocity, eq. (24); $(iii)$ the material derivative of ${\bf U}$ with respect to a given vectorial field, eq. (25).

The quantities expressed by eqs. (24)–(25) are the expansion terms which appear as effective drag terms in the reduced set of equations to be discussed in the following sections. These terms introduce an expansion time scale $\tau_{e}=R_{0}/U_{\infty}$ which is assumed large with respect to the local dynamical time scale, where $U_{\infty}$ is a reference value for $U(R)$, e.g. its asymptotic value. The new expansion terms include effects of the wind acceleration and are an extension of those already present in the EBM.

However, within the acceleration region, where the solar wind crosses the critical sound and Alfvén points, the coupling of fluctuations with the inhomogeneity of background quantities other than ${\bf U}$, e.g. background density, sound and Alfvén speed, may not be neglected. Indeed, this coupling introduces corrections of order of $c_{s}/U$ and $V_{a}/U$ with respect to the expansion terms. For instance, take the induction equation (11b): the term which couples fluctuations and gradient of the background magnetic field is $-{\bf u\cdot\boldsymbol{\nabla}V_{a}}\sim u(V_{a}/R)$. It can be seen easily that the latter is of order $V_{a}/U$ with respect to the expansion term ${\bf\tilde{b}}(\boldsymbol{\nabla}\cdot{\bf U})\sim\tilde{b}(U/R)$, assuming an Alfvénic correlation of the type $u\sim\tilde{b}$. As a consequence, if one wants to describe correctly plasma dynamics in the acceleration region near the critical Alfvén point where $V_{a}=U$, this term must be retained. Similar considerations hold for the pressure, continuity, and momentum equations where the additional coupling of fluctuations with the large scale gradients of the sound speed introduces terms $\sim c_{s}/U$. As we will discuss, all these large scale gradient terms that we introduce in this model are the first necessary corrections to be kept in order to conserve the wave action flux of small amplitude waves. In particular, we will introduce the simplest form of such large scale gradients necessary to conserve the generalized wave action flux of Alfvén waves, and we will consider the case $c_{s}/U\ll 1$.

3. The average solar wind

We prescribe a profile of the solar wind speed $U(R)$, so that the background mass density, (non normalized) pressure and magnetic field evolve self-consistently according to eqs. (9c) together with eqs. (24)–(25). Taking $R_{0}$ as a reference heliocentric distance and indicating with a “star” quantities at $R_{0}$, we obtain the following radial evolution for density, pressure and magnetic field, that we label with an over-bar:

{widetext}

[TABLE]

By introducing the solar wind speed and $R(\bar{t})$, the former slow radial dependence in $R/R_{0}$ in the fixed frame of reference becomes an explicit slow time dependence in $\bar{t}/\tau_{e}$, $\tau_{e}=R_{0}/U_{\infty}$ being the expansion time scale. For example, consider the asymptotic regime of expansion, where $U=U_{\infty}$ and

[TABLE]

In this case, density and pressure decrease respectively as

[TABLE]

[TABLE]

according to the conservation of mass flux and adiabatic cooling, as expected in our model for the solar wind. At small heliocentric distances in the acceleration region the expansion rates are larger, and density and pressure decrease faster.

Within the same regime we get also the following expressions for the magnetic field:

[TABLE]

The scaling of the angle of the Parker spiral with heliocentric distance, which is $\propto R/R_{0}\propto\bar{t}/\tau_{e}$, is thus recovered. The set given by eqs. (26d) for given $U(R)$ will represent our average solar wind. Note that at our order of expansion of $U(r)$, the magnetic field component $\bar{B}_{0,x}$ and $\bar{B}_{0,y}$ respectively mimic the radial and the azimuthal component of the average magnetic field.

4. Derivation of the AEB governing equations

We want to approximate eqs. (11d) in the co-moving frame in order to describe some fundamental properties of fluctuations – in particular of Alfvén waves – related to solar wind inhomogeneity and expansion, while keeping a locally homogeneous and plane geometry.

In order to obtain the evolution equations in the AEB from eqs. (11d) it is therefore necessary not only to introduce explicitly the expansion terms (${\bf A\cdot\boldsymbol{\nabla}U}$ and $\boldsymbol{\nabla}\cdot{\bf U}$), but also to write, as a function of the known quantities $U$, $U^{\prime}$ and $R$, all spatial derivatives of the other background fields (${\bf V_{a}}$, $\rho_{0}$, $c_{s}$) before changing to the co-moving frame. If not, such important terms which couple fluctuations with the gradients of the average solar wind would be lost, as the background is, by construction, locally homogeneous in the (accelerating) expanding box approximation. We will however retain only the first large scale gradient corrections necessary to describe correctly Alfvén waves in the so-called non-WKB approximation, in which the generalized wave action flux, involving both “outwards” and “inwards” waves, must be conserved. To this end, we rely on the following approximations: $(i)$ we assume a radial background magnetic field, since, as explained, the large scale gradient corrections are non-negligible only close the Alfvénic point, where the spiral angle is very small, but equations can be generalized to include terms $\sim V_{a,\varphi}/U$; $(ii)$ we consider the supersonic region and neglect large scale gradient terms of order $c_{s}/U$ involving compressible fluctuations, as well as those large scale gradients in the momentum equation coupling with density perturbations – inclusion of these effects on compressible fluctuations may be deferred to a later work; $(iii)$ we retain the remaining large scale gradients up to first order in the fluctuations’ amplitude. We therefore neglect those terms in the $\boldsymbol{\nabla}\tilde{P}_{T}$ since the latter is essentially nonlinear for Alfvén waves.

Within these assumptions, we obtain the following equations valid in the Sun-centered fixed frame of reference – notice that the large scale advection at the solar wind and Alfvén speed, of the form $U(\partial/\partial R)$ and $V_{a}(\partial/\partial R)$, are implicitly included in the $\boldsymbol{\nabla}$ operator:

{widetext}

[TABLE]

The gradients of the Alfvén velocity can now be taken into account explicitly by expanding its components to first order in $\epsilon$ and $\alpha$, in the same way we expanded the solar wind velocity in eq. (17). The material derivative of ${\bf V_{a}}$, for a generic vector ${\bf A}$, therefore reads

[TABLE]

with

[TABLE]

We now change to the co-moving frame by using the coordinate transformation (18) and (23), eqs. (24)–(25) and eq. (32), and retain the simplest Alfvénic couplings introduced by inhomogeneity that reproduce the known properties of Alfvén waves, in a way which allows to conserve the condition $\boldsymbol{\nabla}\cdot{\bf b}=0$. We therefore impose that the coupling introduced by the Alfvén speed inhomogeneity be between the magnetic field and the incompressible part of the velocity. A more detailed discussion on the problem of the conservation of the zero-divergence of the magnetic field, arising from our simplification of the term ${\bf\boldsymbol{\nabla}\times(u\times B_{0})}$, can be found in Appendix A.

Since the background fields in the AEB are homogeneous in space, we can recast the system of equations in a more compact form. We obtain in this way the AEB model, complemented by eqs. (26d):

{widetext}

[TABLE]

For the sake of notation we have dropped the over-bar, and the projectors $\hat{P}_{u}$, $\hat{P}_{b}$ and the gradient are defined as follows:

[TABLE]

[TABLE]

[TABLE]

In summary, equations (34d) contain the same expansion terms of the previous expanding box model (Grappin & Velli, 1996) extended to include effects of the wind acceleration. In addition, terms of order $V_{a}/U$ with respect to the traditional expansion ones now arise in the induction and momentum equation, coupling the incompressible part of the velocity (${\bf u}^{\text{inc.}}$) and the magnetic field ${\bf\tilde{b}}$. Those terms are the large scale radial advection at the Alfvén speed, corresponding to the second term on the left-hand-side of either equations (this term originates from the implicit advection term $V_{a,r}(\partial/\partial R)$), and the last two terms on the right-hand-side (essentially originating from the inhomogeneity of the Alfvén speed): these are the first necessary corrections to be kept in order to conserve Alfvén wave action flux and to describe correctly the radial evolution of the intrinsic frequency of waves at small distances from the Sun, as we will discuss in Section 6.1.

However, as can be seen, the large scale advection at the Alfvén speed along the radial direction couples the time derivatives of the incompressible part of ${\bf u}$ and ${\bf\tilde{b}}$ – as it must be for Alfvén waves. Although one might be tempted to neglect such a term, therefore greatly simplifying the temporal integration of the AEB equations, Alfvén waves would be described only qualitatively across the Alfvén critical point. In particular, while the amplitude of outwards waves would still peak at the critical point according to the WKB approximation, it would display an exponential radial variation in the neighborhood of the critical point, rather than the algebraic one predicted by linear theory. The correction due to the large scale advection at the Alfvén speed was not taken into account in a previous version of the model (Tenerani & Velli, 2013). In the next section we provide an equivalent form of the AEB which allows to tackle the problem of the temporal integration across the Alfvén critical point of the full equations given in (34d).

5. Model equations and numerical implementation

In order to integrate eqs. (34d) we need to decouple the induction and momentum equation by introducing generalized Elasässer variables involving the vorticity ${\bf\boldsymbol{\eta}=\boldsymbol{\nabla}\times u}$ and the current density ${\bf\tilde{j}=\boldsymbol{\nabla}\times\tilde{b}}$. To this end, we make use of the Helmholtz decomposition ${\bf u=\boldsymbol{\nabla}\chi+\boldsymbol{\nabla}\times V}$, where the scalar and vector potentials $\chi$ and ${\bf V}$ are such that

[TABLE]

The sought set of equations can be found by taking the curl and the divergence of the momentum equation and the curl of the induction equation and, in analogy with the Elsässer variables, by definining ${\boldsymbol{\eta}^{\mp}}=\boldsymbol{\eta}\mp{\bf\tilde{j}}$. Evidently, we drop the large scale gradient terms in the compressible part of the velocity as this does not contribute to shear Alfvén modes. The AEB equations can therefore be rewritten for these new generalized Elsässer variables, the divergence $\xi$, pressure and density as:

{widetext}

[TABLE]

with the following definitions,

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

[TABLE]

We defer to Appendix B the corresponding AEB equations for the two-dimensional case. Notice that $\boldsymbol{\eta}^{\mp}$ represent incompressible Alfvén modes, and that now a coefficient of the form ($1\pm V_{a,x}/U$) multiplies the time derivative of outwards (upper sign) and inwards (lower sign) modes – in analogy with what happens in the fixed frame of reference but with respect to the radial variable. Such a coefficient, while being theoretically correct, introduces a formal, apparent singularity for inwards waves at the Alfvén critical point, where $V_{a,x}=U$, that in general can be removed by imposing proper boundary conditions. Inhomogeneity indeed couples counter-propagating modes, so that an outward propagating mode $\boldsymbol{\eta}^{-}$ forces an “anomalous” $\boldsymbol{\eta}^{+}$ component (which has the same phase of $\boldsymbol{\eta}^{-}$ and hence it is not singular), plus, in general, a reflected inward propagating $\boldsymbol{\eta}^{+}$ (that has a singularity at the critical point), see also Appendix C. Since the latter is moving inwards, it will never reach the critical point and the system is not, in principle, singular. The point is that the AEB is periodic, and the singular reflected component of $\boldsymbol{\eta}^{+}$ propagates within the box while approaching the critical point. This evidently introduces numerical stability problems, if the regularity condition, i.e. that at the critical point $\boldsymbol{\eta}^{+}$ must be given only by the anomalous, forced component, is not imposed. In other words, one must have that both left- and right-hand-side of eq. (39a), lower sign, are equal to zero when $U=V_{a,x}$, and this gives a constraint on the functional form of $\boldsymbol{\eta}^{+}$ at the critical point. In order to satisfy this condition, we have used an Euler semi-implicit method for the time advancement of $\boldsymbol{\eta}^{+}$, in which the propagation ${\bf V_{a}\cdot\boldsymbol{\nabla}}\boldsymbol{\eta}^{+}$ and the linear terms proportional to $\boldsymbol{\eta}^{+}$ are implicit.

This can be seen as follows: after Fourier decomposition, the evolution equation for $\boldsymbol{\eta}_{k}^{+}$ takes the following schematic form, with both linear and nonlinear contributions:

[TABLE]

Now, the Euler semi-implicit scheme with time-step $\Delta t=t_{n+1}-t_{n}$,

{widetext}

[TABLE]

allows to remove the singularity and gives the correct regularity condition on $\boldsymbol{\eta}^{+}$ at leading order in fluctuations’ amplitude. An explicit Runge-Kutta scheme is used for the time advancement of all other fields.

6. Properties of the AEB model and numerical tests

In this section we discuss the properties and limits of validity of the AEB described by eqs. (39d), and we show that it is suitable for the description of small scale plasma dynamics from the sub-Alfvénic region out to the inner heliosphere, provided the wind is supersonic. To this end we analyze the evolution of Alfvénic perturbations and the general dispersion relation of MHD waves as a function of heliocentric distance. For the sake of simplicity, and because here we are mainly concerned with the dynamics close to the Alfvén critical point, we consider fluctuations within a radial magnetic field. In particular, we take a flux tube of fast solar wind with asymptotic speed $U_{\infty}=750$ km/s, and we assume an adiabatic closure, thus the polytropic index is $\gamma=5/3$. In Fig. 2 we show the solar wind speed $U(R)$ (solid black) and the resulting Alfvén (solid blue) and sound (dashed gray) speed as a function of the heliocentric distance $R$. In our solar wind model we have placed the Alfvén critical point at $10\,R_{s}$ and the sonic point at about $3\,R_{s}$. In the following, we discuss via direct analyses that shear Alfvén and oblique magnetosonic modes propagating within this background equilibrium are well described by the AEB model. Next we show a numerical test of an Alfvén wave in parallel propagation.

6.1. Normal mode analysis

As already mentioned in the Introduction, it is known that in the case of an inhomogeneous medium with flow the conservation of wave energy density of a given mode of fluctuation is generalized to the conservation law of wave action $\mathcal{E}/\omega$, where $\mathcal{E}=1/2\rho_{0}u^{2}$ is the wave energy density and $\omega$ its intrinsic frequency (frequency in the plasma rest frame), provided frequencies are larger than $\tau_{e}^{-1}$ (WKB approximation). Within this limit, and for a steady background with spherical symmetry, conservation of the flux $S$ of the wave action across a spherical surface reads $(V_{g}+U)r^{2}\rho_{0}u^{2}/\omega=S$, where $V_{g}$ is the radial component of the group velocity (Bretherton & Garrett, 1969; Jacques, 1977). Such a conservation law leads to a radial variation of the amplitude of $u$ which, in most cases, peaks more or less where the radial group velocity crosses the flow speed. In particular, for outwards Alfvén waves $S=u^{2}(U+V_{a})^{2}/(UV_{a})$, so that their amplitude peaks at the Alfvén critical point, and next it decreases as $V_{a}\sim r^{-1}$ for $V_{a}\ll U$:

[TABLE]

The conservation law of wave action can be extended to the so-called non-WKB limit in which the coupling between counter-propagating modes induced by the background inhomogeneity cannot be neglected (Heinemann & Olbert, 1980). This is necessary especially for the lower frequencies with respect to $\tau_{e}^{-1}$. Although our derivation is rigorously valid only for $\omega\tau_{e}\gg 1$, in practice interesting results may also be derived when larger spatial and temporal scales are included in the simulations (see Grappin & Velli (1996) for the extreme case of co-rotating interaction regions). Therefore demonstrating a greater fidelity also at such lower frequencies is important and, in particular, a correct description of inwards Alfvénic modes is central to any study of turbulence in the solar wind and in particular of wave reflection driven turbulence.

In this regard, consider plane waves in parallel propagation and the two-dimensional AEB equations given in Appendix B. In this case $u_{y}\mp\tilde{b}_{y}=-(i/k)\eta^{\mp}$, where $k=k_{*}(L_{0x}/L_{x})$, and from eqs. (B3b) one recovers the well known expressions for small amplitude shear Alfvén waves ${\bf z}^{\mp}={\bf u\mp\tilde{b}}$ in the non-WKB limit:

{widetext}

[TABLE]

[TABLE]

From eqs. (50)–(51) follows our counterpart of the generalized equation of wave action conservation in the solar wind, where the radial dependence is now transformed into a time dependence,

{widetext}

[TABLE]

After some algebra,  eq. (52) leads to the known conservation law of the total wave action flux $S^{*}$,

[TABLE]

As can be seen, the large scale gradient and the Alfvén speed advection terms, in addition to the standard expansion ones, are necessary for the description of the evolution of energy densities close to the Alfvén critical point. In particular, the term $(1\pm V_{a}/U)$ which multiplies the time derivative ensures the correct radial evolution of the intrinsic frequencies.

Indeed, the evolution of energy densities with heliocentric distance crucially depends on the radial dependence of the intrinsic frequency $\omega$. The dispersion relation of Alfvén waves in the AEB can be easily obtained from eqs. (50)–(51) to leading order in $(k_{*}R_{0})^{-1}$, assuming a solution of the form $A(t)\exp[-i\int^{t}{\omega(t^{\prime})dt^{\prime}}+ik_{*}x]$. This is a standard procedure but a detailed calculation can be found also in Appendix C. We obtain in this way

[TABLE]

where $\Omega=kU$ is the (constant) frequency in the fixed frame of reference. In this regard, recall that, because of the stretching of the box, the derivative in the radial ($x$) direction is rescaled, and as a consequence the radial wave vector is $k=k_{*}(L_{0x}/L_{x})$. It can be easily verified that $kU=const.$, and that it corresponds to the absolute frequency $\Omega$ when $V_{a}/U\ll 1$. Also, $k_{abs}$ can be recognized as the wave vector of Alfvén waves in the fixed frame, and it must be considered as the effective wave vector within the AEB, being the one consistent with the intrinsic dispersion relation. If the term $\pm V_{a}/U$ that multiplies the time derivative was neglected, we would not have the correct radial evolution of $\omega(R)$ close to the critical point. This is what happens for outwards and inwards radial sound waves (represented respectively by $y^{\pm}=u_{x}\pm\tilde{p}$) at the sonic point, for which the same considerations as those for Alfvén waves hold.

To conclude this first part, we show that, as desired, the AEB reproduces the intrinsic frequencies of MHD waves to an excellent approximation, provided $U>c_{s}$. In Fig. 3 we plot the radial evolution of $\omega(R)$ for the shear Alfvén and parallel sound waves (upper panel) and of fast and slow mode waves propagating in the $\{x,y\}$ plane at two different angles (middle and bottom panels). In these plots we compare $\omega(R)$, obtained from the AEB linearized equations, with the theoretical solution, that we label $\omega_{abs}(R)$. The latter is the Doppler shifted frequency in the fixed frame of reference, $\omega_{abs}(R)=\Omega-k_{r,abs}U$, and $k_{r,abs}$ can be found as a function of $R$ from the MHD dispersion relation, using our solar wind model. We compare frequencies at the same values of $k_{x}U=\Omega$ and of perpendicular wave number $k_{y}=n/R$: the black line corresponds to $\omega_{abs}$ (the exact solution), the red-dotted line corresponds to $\omega(R)$ in the AEB, and the blue-dotted line corresponds to the intrinsic frequency in the AEB if the large scale advection correction, $\pm(V_{a}/U)\partial/\partial t$, is neglected. As can be seen, even for relatively large angles the AEB provides an excellent description of frequencies, provided the wind is supersonic, and it even captures the reflection point of fast oblique inwards modes (when the two branches of out- and inwards modes connect). On the contrary, Alfvén and fast modes would be poorly approximated in our regions of interest without the large scale advection term.

6.2. Numerical test: Alfvén wave propagation across the critical point

We consider the same fast solar wind flux tube with radial magnetic field shown in Fig. 2, from an initial heliocentric distance $R_{0}=6\,R_{s}$. The initial values of the background speeds are $U_{*}=390$ km/s, $V_{*a}=750$ km/s and $c_{*s}=100$ km/s. The inverse expansion time is $\tau_{e}^{-1}=U_{\infty}/R_{0}\simeq 2\times 10^{-4}$ s*-1*. We consider a small amplitude wave packet polarized in the $z$ direction with initial conditions

[TABLE]

The amplitude is set to $A/V_{*a}=0.001$, the initial wavelength is $\lambda_{*}=0.75\,R_{s}$, that corresponds to a frequency $\Omega=0.005$ Rad/s, and the width of the envelop is $\Delta=1\,R_{s}$. We choose a localized wave packet instead of a monochromatic wave for illustrative purposes, as this allows to better show the generation of the two components of the inward mode, the forced and the reflected one, and to check their evolution.

In Fig. 4 we show the evolution of the root-mean-square amplitude of ${z}^{-}$ and $z^{+}$ as a function of the average heliocentric distance $R$. We compare our numerical solution with that obtained by numerically solving the boundary value problem of a stationary, monochromatic Alfvén wave in the Eulerian frame at the same $\Omega$, shown as a dotted line. As can be seen in the upper panel, our results for $z^{-}$ are in agreement with theory, and the amplitude is given to an excellent approximation by $|z^{-}|^{2}\propto UV_{a}/(U+V_{a})^{2}$. The envelop of $|z^{+}|^{2}$ also evolves following the theoretical solution (bottom panel). The presence of rapid oscillations at the beginning is due to beats between the forced and the reflected components of $z^{+}$. This stems from the fact that the AEB solves an initial value problem. Therefore, if at $t=0$ only $z^{-}$ is given, then a forced and a reflected $z^{+}$ are generated accordingly, causing beats when they superpose. In Appendix C we derive the analytical solution of the wave equations (50)–(51) via a two-scale analysis that illustrates this point. The numerical scheme however behaves like an effective filter and it damps the singular reflected component while approaching the critical point, at $10\,R_{s}$. In this way the reflected mode does not propagate beyond the critical point, as expected, preserving the causal disconnection between the sub-Alfvénic and the super-Alfvénic corona.

This can be seen also by inspecting Fig. 5, where the contours of $z^{+}$ in $\{x,R\}$ space are displayed. One can see that there are two components propagating at two different phase speeds that decouple with time (heliocentric distance). These two components are tracked by the yellow-dotted lines corresponding to the characteristics of the system: the forced component of $z^{+}$, by definition, follows the characteristic of $z^{-}$ given by

[TABLE]

while the characteristic of the reflected component is given by

[TABLE]

In Fig. 6 we show for reference the waveform of $z^{+}$ at $t\Omega=0.9$ ($R=6.1\,R_{s}$, black line) and at $t\Omega=17$ ($R=8.3\,R_{s}$, blue thin line) that shows the separation of the two components of ${z^{+}}$ during the evolution of the initial wave packet, and the damping of the singular component.

7. Summary and conclusions

We have presented the accelerating expanding box that extends and improves the previous expanding box model in two directions. The original EBM is suitable for describing plasma dynamics in the co-moving frame of the solar wind when a uniform flow is assumed at large distances from the Sun, i.e., when the flow is both supersonic and super-alfvénic, $U\gg c_{s}$ and $U\gg V_{a}$. The present model is an extension and improvement in that:

$i)$ We allow for a non uniformity of the solar wind speed at first order inside the box (i.e. we retain up to first derivative of $U(R)$). This causes an additional stretching of the box in the $x$ (radial) direction and additional terms $\propto U^{\prime}$ coming from the ${\boldsymbol{\nabla}\cdot{\bf U}}$ and from the material derivative of ${\bf U}$. These terms cause the average pressure, density and azimuthal magnetic field to decrease with heliocentric distance faster than in the EBM.

$ii)$ In order to be consistent with the description of the dynamics close and across the Alfvén critical point, we retain terms of order $V_{a}/U$ (neglecting effects of the Parker spiral close to the Sun), which are not taken into account in the standard EBM. These new terms modify the decay rates of fluctuations only, and introduce additional reflection terms that allow for the conservation of the adiabatic invariants of Alfvén waves, i.e. the generalized wave action flux, given in eq. (53).

To summarize, the AEB model is described by eqs. (34d), that have been re-written in terms of generalized Elsässer variables in eqs. (39d). The latter formulation has the advantage of separating the compressible and incompressible parts of the velocity, as well as decoupling forward and backward Alfvénic modes, necessary for the numerical integration of the equations across the Alfvén critical point. As we have shown, the model equations properly describe the Alfvén critical point crossing without unphysical accumulation of reflected modes at the boundary between sub- and super-Aflvénic corona. Finally, linear analysis of the normal modes of the AEB shows that the theoretical evolution of the intrinsic frequencies of MHD waves is recovered in the supersonic wind.

A.T. thanks D. Del Sarto for his useful comments and discussions. This work was supported by the NASA program LWS, grant NNX13AF81G as well as the NASA Solar Probe Plus Observatory Scientist grant. This work used the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by National Science Foundation grant number ACI-1053575.

Appendix A Large scale gradients in the induction equation

In this Appendix we discuss the approximations done to derive the large scale gradients in the induction equation, and we show how to obtain eq. (34b) from eq. (31b) in a form consistent with $\boldsymbol{\nabla}\cdot{\bf b}=0$. We therefore consider only ${\bf\boldsymbol{\nabla}\times u\times B_{0}}$, since both the nonlinear (${\bf\boldsymbol{\nabla}\times u\times b}$) and the expansion (${\bf\boldsymbol{\nabla}\times U\times b}$) terms do not affect the conservation of the magnetic field divergence.

We develop to first order in $\epsilon$ and $\alpha$ vectors, derivatives and average magnetic field, so as to obtain

[TABLE]

[TABLE]

[TABLE]

where

[TABLE]

In writing eq. (A3) we have expanded to linear order ${\bf B_{0}}(r)={\bf B_{0}}(R)+{\bf B_{0}}^{\prime}(R)\delta x$, and ${\bf B_{0}}^{\prime}(R)$ is expressed by making use of the following evolution equations of magnetic field in the solar wind,

[TABLE]

The resulting components of ${\bf\boldsymbol{\nabla}\times u\times B_{0}}$ are:

[TABLE]

[TABLE]

[TABLE]

After collecting eqs. (A6)–(A8) and neglecting the remaining higher order terms, we obtain the following expressions for ${\bf\boldsymbol{\nabla}\times u\times B_{0}}$ in cartesian coordinates:

[TABLE]

[TABLE]

[TABLE]

The reduced components of ${\bf\boldsymbol{\nabla}\times u\times B_{0}}$ given in (A9)-(A11) have zero divergence, also in the simpler case with $\bar{B}_{0,y}=0$ that is considered specifically in this paper. It is also easy to recognize the origin of each term: the first one represents ${\bf B_{0}}({\bf\boldsymbol{\nabla}\cdot u})$; the second and third terms correspond to the advection at the Alfvén velocity ${\bf B_{0}\cdot\boldsymbol{\nabla}u}$; the last ones, of the form $\sim\bar{B}_{0}u/R$, are the large scale gradients of the magnetic field originating from ${\bf u\cdot\boldsymbol{\nabla}B_{0}}$. Although the decomposition given in eqs. (A9)-(A11) is rigorous to our order of approximation, there is now an explicit dependence on local coordinates, whereas we want to recover local homogeneity and periodicity. We therefore neglect first order geometrical terms $\propto y/R,\,z/R$, the local inhomogeneity of ${\bf B_{0}}$, which is $\propto\delta x/R$, and $2u_{x}/R$ (in any case this term does not contribute to linear Alfvén waves), while we retain the large scale gradients in the perpendicular directions and advection:

[TABLE]

The expression for the normalized magnetic field consistent with eq. (A12) is

[TABLE]

Change to the co-moving frame is now straightforward with the transformation

[TABLE]

which allows us to write

[TABLE]

The modified induction equation as given by (A15) does not conserve the divergence of the magnetic field because the last two terms depend on the velocity ${\bf\bar{u}}$. In order to conserve the divergence-free condition, we have to impose that such couplings involve the incompressible part of the velocity ${\bf\bar{u}}^{\text{inc}}$, as indeed is the case for Alfvén waves. Notice that our expanding box coordinates are time-dependent, so that the divergence of eq. (A15) leaves a residual term $U^{\prime}/U\bar{V}_{a,x}\partial\bar{u}_{x}/\partial\bar{x}$ that comes from the time derivative of coordinates. Although for small amplitude Alfvén waves this equals zero (in our approximation of radial magnetic field close to the Sun), we formally balance it by adding $U^{\prime}/U\bar{V}_{a,x}u_{x}$ along the ${\bf\hat{x}}$ direction.

Appendix B Two-dimensional model equations

For the two-dimensional case, velocity and magnetic field can be decomposed via two scalar potentials $\chi$ and $\phi$ plus a component along ${\bf\hat{z}}$:

[TABLE]

[TABLE]

In this case we integrate the variables $\eta^{\mp}=\eta\mp\tilde{j}$ and $z^{\mp}=u_{z}\mp\tilde{b}_{z}$, with $\tilde{j}=j/\sqrt{4\pi\rho_{0}}$ and $\tilde{b}_{z}=b_{z}/\sqrt{4\pi\rho_{0}}$:

[TABLE]

Appendix C Analytical solution of forced Alfvén waves in the AEB

We consider here a shear Alfvén wave in a radial magnetic field ${\bf B}_{0}$, with solar wind speed $U$ and Alfvén speed $V_{a}$. The AEB counterpart of the equations for outwards and inwards waves are:

[TABLE]

[TABLE]

We derive below the analytical solution of eqs. (C1)–(C2) with given initial condition $z^{-}(x,0)=A^{-}(0)\exp(ik_{*}x)$, $z^{+}(x,0)=0$. We employ a two-scale approach in which the small parameter is $\epsilon=1/(k_{*}R_{0})$, where $k_{*}$ is the initial wave number, $k=k_{*}(L_{0x}/L_{x})$, and we recall that $\tau_{e}=R_{0}/U_{\infty}$. For the sake of simplicity we consider monochromatic waves, and according to the two-scale analysis we will look for first order solutions in $\epsilon$ of the form

[TABLE]

[TABLE]

where the subscript “f” and “h” stand for the forced and the homogeneous part of $z^{+}$.

C.1. Outward wave

By substitution of (C3) into (C1) and by introducing normalized variables

[TABLE]

we obtain

[TABLE]

At order $\epsilon^{0}$ one gets an equation for the phase of the outwards wave

[TABLE]

where $k_{abs}$ is the wave vector in the absolute frame of reference,

[TABLE]

At the order $\epsilon$ one obtains an equation for the amplitude

[TABLE]

that can be integrated numerically by changing variable to $dR=Udt$, and has solution, as expected,

[TABLE]

where $C_{0}$ is determined from the initial condition.

C.2. Inward forced wave

Assuming that there are no inwards waves at $t=0$, we look for a solution to eq. (C2) that in general is given by a forced wave $z^{+}_{f}$, corresponding to the particular solution of eq. (C2), plus a reflected wave $z^{+}_{h}$, the latter being the solution of the homogeneous part of eq. (C2). The forced wave is of the form

[TABLE]

By inserting (C11) in eq. (C2) and with the same ordering procedure as above, one gets to first order in $\epsilon$ an algebraic equation for $A^{+}$ (remember that $k=k_{*}L_{0x}/L_{x}$):

[TABLE]

The amplitude of $z_{f}^{+}$ is therefore

[TABLE]

The homogeneous solution $z_{h}^{+}$ is necessary in order to impose the initial condition, and must satisfy the following equation,

[TABLE]

and hence, to leading order in $\epsilon$,

[TABLE]

By imposing the initial condition $z^{+}_{h}(x,0)+z^{+}_{f}(x,0)=0$ we obtain a condition on $\phi$ and $A^{+}_{h}$,

[TABLE]

The complete solution of the $z^{+}$ mode is,

[TABLE]

from which it can be seen that the total amplitude of $z^{+}$ displays beating oscillations at a frequency corresponding to the phase difference between the forced and the reflected waves:

[TABLE]