Investigating `dark' energy in the solar corona using forward modeling of MHD waves

TL;DR

This study uses 3D MHD simulations to demonstrate that the apparent energy discrepancy in solar coronal Alfvénic waves is due to line-of-sight superposition effects, revealing that the true wave energy is hidden in non-thermal line widths.

Contribution

The paper provides the first comprehensive simulation-based evidence that line-of-sight superposition accounts for the observed energy discrepancy, eliminating the need for additional 'dark' energy sources.

Findings

LOS superposition reduces Doppler velocities and increases non-thermal line widths.

Simulated correlation matches observed wedge-shaped pattern.

True wave energy is 0.2-1% of observed energy, explaining the discrepancy.

Abstract

It is now well established that the Alfv\'enic waves are ubiquitous in the solar corona. However, the Alfv\'enic wave energy estimated from the Doppler velocity measurements in the corona was found to be four orders of magnitude less than that estimated from non-thermal line widths. McIntosh & De Pontieu (2012) suggested that this discrepancy in energy might be due to the line-of-sight (LOS) superposition of the several oscillating structures, which can lead to an underestimation of the Alfv\'enic wave amplitudes and energies. McIntosh & De Pontieu (2012) termed this coronal `dark' or `hidden' energy. However, their simulations required the use of an additional, unknown source of Alfv\'enic wave energy to provide agreement with measurements of the coronal non-thermal line widths. In this study, we investigate the requirement of this unknown source of additional `dark' energy in the solar corona using gravitationally stratified 3D magnetohydrodynamic (MHD) simulations of propagating waves. We excite the transverse MHD waves and generate synthetic observations for the Fe XIII emission line. We establish that the LOS superposition greatly reduces the Doppler velocity amplitudes and increases the non-thermal line widths. Importantly, our model generates the observed wedge-shaped correlation between Doppler velocities and non-thermal line widths. We find that the observed wave energy is only 0.2-1\% of the true wave energy which explains 2-3 orders of magnitude of the energy discrepancy. We conclusively establish that the true wave energies are hidden in the non-thermal line widths. Hence, our results rule out the requirement for an additional `dark' energy in the solar corona.

Equations

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

$$\frac{\partial (\rho \mathbf{v})}{\partial t}+\nabla .(\rho \mathbf{vv}-\mathbf{BB}/{\mu }_0)+\nabla (p+{\mathbf{B}}^2/2{\mu }_0)-\rho \mathbf{g}=0,$$

$$\frac{\partial E}{\partial t}+\nabla .(\mathbf{v}E-BB.v/{\mu }_0+\mathbf{v}.(p+{\mathbf{B}}^2/2{\mu }_0))-\rho v.g=0,$$

$$\frac{\partial B}{\partial t}-\nabla \times (\mathbf{v}\times \mathbf{B})=0,$$

$$\nabla .\mathbf{B}=0,$$

$$\rho (x,y,z)=\left({\rho }_0+\sum _{i=0}^{50}{A}_i{\exp }^{-[(y-y_i{)}^2+(z-z_i{)}^2]/2{\sigma }_i^2}\right){\exp }^{-x/H(y,z)}.$$

$$v_y(x=0,t)=\sum _{i=1}^{10}{U}_i\sin ({\omega }_{i}t),$$

$$v_z(x=0,t)=\sum _{i=1}^{10}{V}_i\sin ({\omega }_{i}t),$$

$$v_{rms}={⟨\sqrt{\frac{{\sum }_{t=0}^{T-1}[v_y^2(x=0,y,z,t)+v_z^2(x=0,y,z,t)]}{T}}⟩}_{y,z};T=50.$$

$$ϵ(x,y,z,t)=\frac{A_b}{4\pi }{n}_e^2(x,y,z,t)G(n_e,T),$$

$$I(x,z,t)={\int }_{\lambda }{I}_{\lambda }(x,z,t)d\lambda ,$$

$${\lambda }_D(x,z,t)=\frac{{\int }_{\lambda }\lambda I_{\lambda }(x,z,t)d\lambda }{{\int }_{\lambda }{I}_{\lambda }(x,z,t)d\lambda }-{\lambda }_0,$$

$$\sigma (x,z,t)=\sqrt{\frac{{\int }_{\lambda }(\lambda -{\lambda }_D{)}^2{I}_{\lambda }(x,z,t)d\lambda }{{\int }_{\lambda }{I}_{\lambda }(x,z,t)d\lambda }}.$$

$$I_{\lambda }(x,z^{\prime },t_n)=\sum _{j=1}^{100}{I}_{\lambda }(x,z_j^{\prime },t_{0,j}+n\delta t);1⩽n⩽20,$$

$$rms{v}_D(x,z^{\prime })=\sqrt{\frac{{\sum }_{n=1}^{20}{v}_D(x,z^{\prime },t_n{)}^2}{20}},$$

$$mean{\sigma }_{nt}(x,z^{\prime })=\frac{{\sum }_{n=1}^{20}{\sigma }_{nt}(x,z^{\prime },t_n)}{20}.$$

$$E(x,y,z,t)=\frac{\rho (x,y,z,t)v^2(x,y,z,t)}{2}+\frac{B(x,y,z,t{)}^2}{2\mu }+\frac{p(x,y,z,t)}{\gamma -1}+\rho (x,y,z,t)\varphi (x,y,z,t),$$

$$E=⟨E(t){⟩}_t=⟨\frac{1}{V}{\int }_V(E(x,y,z,t)-E(x,y,z,0))d{V}^{\prime }{⟩}_t.$$

$$F=⟨F{⟩}_t=⟨\frac{1}{A}{\int }_{S}E(0,y,z,t)v_{g}d{A}^{\prime }{⟩}_t,$$

Full text

Sun: Corona, Sun: waves, Sun: magnetohydrodynamics

††journal: ApJ

1 Introduction

The solar corona is heated to millions of Kelvin, with the mechanism responsible for this heating having evaded researchers for decades (Walsh & Ireland, 2003; Klimchuk, 2006; Parnell & De Moortel, 2012; Cranmer & Winebarger, 2018). One of the possible mechanisms of heating is the dissipation of magnetohydrodynamic (MHD) waves in the solar atmosphere (Klimchuk, 2006; Arregui, 2015). MHD waves and their different wave modes have been ubiquitously observed in the solar atmosphere by both space and ground-based instruments (Banerjee et al., 2007). It was suggested that hydromagnetic waves in the solar atmosphere can produce non-thermal broadening of the emission lines (Hollweg, 1973; Van Doorsselaere et al., 2008). The evidence of the broadening of transition region emission lines in the quiet sun and coronal holes was first provided by the spectrograph on the Skylab (Doschek et al., 1976a, b; Feldman et al., 1976). Early observations of non-thermal broadening of coronal emission lines (formed at temperatures $\sim$ 106 K) by Hassler et al. (1990) pointed to the possible existence of the Alfvén waves in the solar atmosphere. Later, Banerjee et al. (1998) and Doyle et al. (1998) used the Solar Ultraviolet Measurements of Emitted Radiation (SUMER) spectrometer on-board Solar and Heliospheric Observatory (SOHO) to study the variation of line widths of the Si VIII emission line at different locations in northern and southern polar coronal holes. These authors observed that the line width of the Si VIII emission lines increased from 27 to 46 km s*-1* when the distance above the solar limb increased to 17 Mm to 175 Mm, respectively. They computed the energy flux in the Alfvén waves and found that it is slightly less than the energy flux (8$\times$105 ergs cm*-2* s*-1*) required to balance total coronal energy losses in coronal holes (Withbroe & Noyes, 1977). Similar studies on the nature of non-thermal broadening of coronal emission lines were carried out using the Coronal Diagnostic Spectrometer (CDS) (O’Shea et al., 2005), Ultraviolet Coronagraph Spectrometer (UVCS) on-board SOHO (Ofman & Davila, 1997a; Kohl et al., 1999; Abbo et al., 2016) and the Extreme Ultraviolet Imaging Spectrometer (EIS) on-board Hinode (Doschek et al., 2007; Banerjee et al., 2009; Hahn et al., 2012). The estimated energy flux was found to be just enough to balance the energy losses in the polar open-field regions of the solar corona. Additionally, several 1D (Lau & Siregar, 1996; Orta et al., 2003; Suzuki & Inutsuka, 2006; Cranmer et al., 2007; Oran et al., 2013; van der Holst et al., 2014; Oran et al., 2017), 2.5D (Ofman & Davila, 1997b, 1998), and 3D (using reduced MHD; van Ballegooijen et al., 2011, 2017) wave-based models driven purely by Alfvén waves have been somewhat successful in explaining the large non-thermal widths of coronal emission lines and acceleration of fast solar winds in open and closed magnetic field regions.

Resolved measurements of the propagating Alfvénic waves came from observations of the chromosphere and transition region using imaging data from the Solar Optical Telescope (SOT) on board Hinode (De Pontieu et al., 2007) and the Atmospheric Imaging Assembly (AIA) on-board the Solar Dynamics Observatory (SDO) (McIntosh et al., 2011). These authors reported waves with amplitudes of $\sim$ 20 km s*-1*, suggesting that they are capable of providing the energy flux of 100-200 W m*-2* to accelerate fast solar wind and balance radiative losses in quiet corona. The coronal counterpart to these propagating Alfvénic waves was observed using Doppler velocity data (Tomczyk et al., 2007, 2008; Tomczyk & McIntosh, 2009; Morton et al., 2015) from the Coronal Multi-Channel Polarimeter (CoMP Tomczyk et al., 2008), and through direct measurements with SDO/AIA (Thurgood et al., 2014; Weberg et al., 2018; Morton et al., 2019). Both decaying (Nakariakov et al., 1999; Aschwanden et al., 2002) and decayless (Nisticò et al., 2013; Anfinogentov et al., 2013) transverse oscillations have been observed in the solar atmosphere through direct imaging. The CoMP data can also measure coronal line widths, providing estimates of the non-thermal component that is comparable to previously reported values. Surprisingly, the measured Doppler velocity fluctuations only have amplitudes of $\sim$0.5 km s*-1*, which suggests the Alfvénic waves only have an energy flux of $\sim$0.01 W m*-2*. Furthermore, in a comparative study, Morton et al. (2015) reported an average velocity amplitude of $\sim$14 km s*-1* measured with SDO at $1.01~{}R_{\odot}$, which is at a lower height in the corona than where CoMP Doppler velocities are measured. Thus, there seems to be an apparent discrepancy in the wave energy estimated using the CoMP Doppler velocities compared to those estimated using non-thermal line widths and the imaging measurements from SOT and AIA.

This discrepancy was investigated by McIntosh & De Pontieu (2012), who used a Monte Carlo method to forward model the emission spectrum generated by several oscillating structures (termed as ‘threads’) along the line-of-sight (LOS). We note that the simulations were effectively a toy model of how Alfvénic waves would impact spectral lines and no MHD simulations were used. It was suggested that small Doppler velocities and large non-thermal line widths are the consequences of unresolved swaying motions of threads along the LOS due to optically thin solar corona. In addition to this, the authors used CoMP to demonstrate a correlation between root mean square (rms) Doppler velocities and mean non-thermal line widths that appeared as a wedge-shape. The authors could explain the wedge-shaped correlation using their model, but, only by including an additional component of non-thermal broadening whose origin was not known. Finally, these authors suggested that the ‘dark’ or ‘hidden’ energy, which is not observed by direct imaging, could be residing in the non-thermal line widths. The effects of superposition of coronal loops along the LOS on the wave amplitudes was also investigated using a 3D MHD model (De Moortel & Pascoe, 2012). In this study, the authors found that the kinetic energy measured from the LOS Doppler velocities is an underestimate of the total kinetic energy present in the model.

In spite of the developments outlined above, little work has been done to investigate the LOS effects on Doppler velocities and line widths in the solar corona in the context of the CoMP. This displays the need for an in-depth study of the wave propagation in the solar corona using MHD models, examing whether they can generate the observed values of non-thermal widths of emission lines and their variation with height through the corona. Additionally, the model should also be able to produce a wedge-shaped correlation between Doppler velocities and non-thermal line widths.

In this work, we investigate the correlation between rms Doppler velocities and non-thermal line widths in open magnetic field regions using 3D MHD simulations of propagating waves. Further, we explore the requirement of an additional non-thermal broadening, which is needed to explain the wedge-shape correlation in McIntosh & De Pontieu (2012). We forward model our MHD simulations for the Fe XIII emission line and examine the variation of the non-thermal line widths with height in the solar atmosphere. The solar plasma is inhomogeneous and gravitationally stratified, which leads to the reflection (Ferraro & Plumpton, 1958; Hollweg, 1978) and non-linear interaction of waves (Matthaeus et al., 1999; Cranmer et al., 2007) propagating through it. Recently, Magyar et al. (2017) have studied the effects of perpendicular inhomogeneities on unidirectionally propagating MHD waves and reported that such inhomogeneities lead to a generalised phase mixing that generates a turbulence-like behaviour. These authors termed this ‘uniturbulence’ since it is produced by unidirectionally propagating Alfvénic waves. Karampelas et al. (2019) have investigated the effects of gravitational stratification on the heating of coronal loops using 3D MHD simulations. These authors have reported that the inclusion of gravity increased the average temperature near the footpoint and the apex of the coronal loop compared to the simulations where gravity was excluded (see also Karampelas et al., 2017). Since the gravity and plasma inhomogeneities alter the wave propagation, it is crucial to study their effects on the observed properties of waves. Thus we extend the MHD model of unidirectional propagating waves in a perpendicularly inhomogeneous plasma (Magyar et al., 2017; Magyar & Van Doorsselaere, 2018) by including gravity and investigate the effects LOS superposition of waves on the observables such as Doppler velocities and line widths. The paper is structured as follows. The observations using CoMP are reported in section 2. The choice of parameters (in the model) and numerical set-up is discussed in Section 3. Forward modeling using Fe XIII emission line is described in Section 4. Section 5 outlines the analysis and results which is followed by a discussion and conclusions in Section 6.

2 Observations

Here we use observations taken with the CoMP instrument on the 27 March 2012. The data set has been used previously in Morton et al. (2015, 2016, 2019), where discussion of any additional post-processing to the Level 2 data files are given in detail. CoMP measures intensity at three positions across the Fe XIII 10747 Å line, from which estimates for the Doppler shifts and line widths can be made. The rms Doppler velocities and non-thermal line widths are calculated from these quantities, where we follow the prescription in McIntosh & De Pontieu (2012) to calculate the non-thermal component of the line widths (and permit direct comparison between results).

In Figure 1, we plot a two-dimensional histogram of the non-thermal line width against rms Doppler velocity for the entire field-of-view (FOV). It is seen that a similar wedge shape is obtained to that in McIntosh & De Pontieu (2012). There are some differences between the location of the wedge and its exact shape. However, this variation can likely be attributed to differences in the coronal magnetic field and plasma conditions, which leads to different wave properties (e.g, Morton et al., 2016, 2019).

It is well known that the density of the corona decreases as a function of height. CoMP data provides measurements between $1.05~{}R_{\odot}$ to $1.3~{}R_{\odot}$, which corresponds to approximately 2 pressure scale-heights for a 1.6 MK plasma. Furthermore, it is known that the amplitude of Alfvénic modes depends on the density, which from WKB theory goes as $v\propto\langle\rho\rangle^{-1/4}$. Hence, we should expect that Doppler velocities and non-thermal line widths measured with CoMP show some dependence with height. An increase in amplitude with height is visible in the full FOV images of Figure 1 in McIntosh & De Pontieu (2012). Moreover, previous observations of coronal holes, with CoMP and other spectrometers, have demonstrated an increase in amplitude for both rms velocity and non-thermal line widths (for altitudes $<1.2R_{\odot}$) that are broadly in agreement with the WKB theory (e.g, Hahn et al., 2012; Morton et al., 2016). However, we note, that the observed increase in amplitude in other regions of the corona (i.e. away from coronal holes) will not match the expected amplification for WKB Alfvénic waves due to density stratification, as the waves are known to undergo some form of damping in the corona (e.g. Verth et al., 2010, Tiwari et al., 2019). It is worth highlighting that kink waves in compressible plasmas have mixed properties that are similar to the surface Alfvén waves (Goossens et al., 2009, 2012), thus several authors used the term ‘Alfvénic’ to describe such waves. In general such waves are transverse and largely incompressible. In this work, we call these waves as transverse MHD waves.

To examine how the variation in wave amplitude with height influences the wedge shape, we show a portion of the data from a coronal hole region in the lower panel of Figure 1. In this region the magnetic field is almost radial, and hence enables us to unambiguously show the influence of change in density with height on both the rms velocity and non-thermal line widths. Figure 1 reveals that part of the contribution to the wedge shape, affecting both the range of non-thermal width values and the positive correlation with rms Doppler velocity, is due to the dependence of the transverse MHD wave amplitude on height111The larger values of non-thermal line widths at 78-94 Mm are likely due to the influence of scattered light from the CoMP occulter.. We note that this consideration was neglected from the analysis of McIntosh & De Pontieu (2012) and, as such, raises questions about the inferred values of wave properties from their Monte Carlo simulations.

3 Numerical setup

To understand the transverse MHD wave propagation in open-field regions, we perform an ideal 3D MHD simulation using MPI-AMRVAC that solves hyperbolic partial differential equations in near conservative form (Porth et al., 2014). The following equations are solved in cartesian geometry for a grid size of 128$\times$512$\times$512 that span 50 Mm$\times$5 Mm$\times$5 Mm (see Figure 2 a).

[TABLE]

where $\bf{B}$ is the magnetic field, $\bf{v}$ is the plasma velocity, $\rho$ is the density, $p=\rho\frac{k_{b}}{\mu m_{H}}T$, is the gas pressure and $E$ is total energy density defined as, $E=\frac{p}{\gamma-1}+\frac{\rho{\bf v}^{2}}{2}+\frac{{\bf B}^{2}}{2\mu_{0}}$. Furthermore, $\mu_{0}$ is the magnetic permeability in the free space, g is the acceleration due to gravity pointing along the negative $x$ axis, $\mu$ is taken to be 0.6 for coronal abundance, $m_{H}$ is the mass of the proton, $k_{b}$ is the Boltzmann constant, and $\gamma$ is chosen to be 5/3.

3.1 Initial Conditions

A uniform grid without any mesh refinement is employed for performing the simulations. The spatial resolution along the $x$, $y$, and $z$ axis is 0.39 Mm, 0.01 Mm, and 0.01 Mm, respectively, and the $x$ axis defines the vertical direction. Initially, the simulations are set up assuming a vertically isothermal atmosphere, meaning the temperature is constant with height. Thus, $\rho$ is an exponentially decaying function of the height due to the gravitational stratification. Further, in the initial setup, inhomogeneities in density are randomly placed transverse ($y-z$ plane) to the direction of the magnetic field according to the following equation,

[TABLE]

$H(y,z)$ is the scale height that depends on the temperature which is different at different locations in the $y-z$ plane. We choose $\rho_{0}=2\times$10*-13* kg m*-3*. The magnitude of the inhomogeneity is given by $A_{i}$, which is randomly drawn from a uniform distribution of [0, 5] $\rho_{0}$. The spatial extent of the inhomogeneity is controlled by $\sigma_{i}$, which is randomly chosen from [0, 250] km. The spatial location of the inhomogeneity $y_{i}$ and $z_{i}$ are chosen within the simulation domain i.e, drawn randomly from a uniform distribution of [-2.5, 2.5]. We set the gas pressure in the $y-z$ plane to constant, hence Eq. 2 also indirectly determines the initial temperature structure. The initial magnetic field is assumed to be uniform and vertical, $\textbf{B}=B\bf{\hat{x}}$, and its strength is prescribed to be 5 G. This value of the magnetic field strength is typical for coronal holes (Hollweg, 1990). We assume a low beta ($\beta$) plasma with $\beta=0.15$.

Longitudinal and transverse sections of the initial setup are shown in panels (b) and (c) of Figure 2. Figure 3 shows the variation of the density with height in the initial setup at the location ($y=0$, $z=-0.9$) marked with the dashed line in Figure 2 (b). Corresponding to the location of the dashed line in Figure 2, the scale height is 40 Mm (see, Figure 3).

3.2 Initial evolution

Since the initial setup described above is not in the pressure equilibrium, we evolve the simulations initial setup for $\sim$100 s before implementing any driving. We use a TVD second-order solver and a Woodward slope limiter to solve for $p$, ${\bf v}$, $\rho$, and ${\bf B}$. Moreover, Powell’s scheme is employed to ensure a divergence-free magnetic field. For this stage we employ open boundary conditions across all boundaries, such that any MHD waves generated can leave the simulation domain. This evolutionary step allows the initial state of the system to relax to a state of pressure equilibrium. As the plasma relaxes, the density inhomogeneities and the magnetic field expand in response (see, Figure 3). The magnetic field becomes vertically stratified and concentrates at the density inhomogeneities due to the lower gas pressure in these regions. From here on, we use $t=0$ to denote the time when the simulations have reached at the pressure equilibrium. An example of the temperature and density structure in the $y-z$ plane is shown in the top panel of Figure 4.

3.3 Wave excitation

Once the simulation reach pressure equilibrium, we excite transverse MHD waves by driving the entire bottom boundary ($x$=0 Mm) perpendicular to the direction of the background magnetic field (i.e, in $y-z$ plane). In contrast to the initial evolution, boundaries in the $y$ and $z$ directions are now set to be periodic. However, the top boundary of the simulation domain is kept open, so that transverse MHD waves can leave the domain.

The bottom boundary is driven uniformly by a velocity driver composed of a superposition of ten velocity drivers with different periodicities. The $y$ and $z$ components of the velocity driver are given by the following equations,

[TABLE]

where periods (hence, $\omega$’s) are chosen from the observed log-normal distribution (Thurgood et al., 2014; Weberg et al., 2018; Morton et al., 2015, 2019). The magnitude of $U_{i}$ and $V_{i}$ are randomly chosen from a uniform distribution of $[-U_{o},U_{o}]$. In this work, we perform simulations for three different values of $U_{o}$: $\frac{5}{\sqrt{2}}$ km s*-1*, $\frac{11}{\sqrt{2}}$ km s*-1*, and $\frac{22}{\sqrt{2}}$ km s*-1*. The root mean square (rms) value of velocity driver ($v_{rms}$) averaged over the entire bottom boundary corresponding to above three cases are $\sim$7 km s*-1*, $\sim$15 km s*-1*, and $\sim$26 km s*-1*, respectively. The $v_{rms}$ averaged over the entire bottom boundary ($x=0,y,z$) is computed using the following relation.

[TABLE]

Here $T$ is the total duration of the simulations. All wave simulations are performed for 1000 s with a cadence of 20 s, thus giving 50 snapshots for each run. Angular brackets represents the average over $y-z$ plane. Henceforth, throughout the manuscript, $v_{rms}$ represents the rms velocity of the driver at the bottom boundary. Due to the difference in periodicity and amplitudes, the resultant velocity field has a varying phase with time. An animation corresponding to the panel (c) of Figure 2 is available for the case where the lower boundary of the simulation domain is driven by a velocity driver with $U_{o}$=$\frac{11}{\sqrt{2}}$ km s*-1*. The average sound speed ($c_{s}$) and Alfvén speed (vA) at $t=0$ are $\sim$ 120 km s*-1* and 500 km s*-1* respectively. Since $v_{rms}$ $<$ $c_{s}$ $<$ vA for all cases, the excited waves are in a linear regime and weak compressible limit (or largely incompressible).

4 Forward modeling with FoMo

To compare the results of the simulations with observations taken with the CoMP (Section 2), we need to convert physical variables obtained from the simulations (e.g., density, temperature ($T$), velocity) to spectroscopic observables such as specific intensity (function of wavelength) for the Fe XIII emission line centred at 10749 Å, from which the Doppler velocities and line widths can be derived. Since the solar corona is optically thin, the specific intensity ($I$) is computed by adding the emission ($\epsilon$) of different structures along the LOS. We use the FoMo tool developed for the forward modeling of the optically thin emission from the coronal plasma (Antolin & Van Doorsselaere, 2013; Van Doorsselaere et al., 2016), specifically FoMo-C (adapted for C++), to compute intensities, $I(\lambda,x,z,t$), and generate synthetic observations comparable to CoMP. The line formation temperature and thermal width of the Fe XIII emission line are $\sim$1.6 MK and $\sim$21.78 km s*-1* respectively. For calculating the emission in Fe XIII (10749 Å), we use the CHIANTI atomic database (Dere et al., 1997), a coronal abundance of the Fe XIII relative to the hydrogen, and assume that the emitting plasma is at ionisation equilibrium. Specifically, the contribution function, $G(n_{e},T)$, for the Fe XIII is computed using the CHIANTI. Further, the emission is computed at every location using the following relation;

[TABLE]

where $A_{b}$ is the coronal abundance and $n_{e}$ is the electron density (see Antolin & Van Doorsselaere, 2013; Van Doorsselaere et al., 2016). Computation of $G(n_{e},T)$ requires the knowledge of the rates of electron excitation/de-excitation, proton excitation/de-excitation, and photoexcitation/stimulated emission (see Dere et al., 1997; Young et al., 2003). We note that photoexcitation does not influence our results (see Appendix A).

The upper right panel in Figure 4 shows the $\epsilon$ computed using Eq. 5 at $x$=20 Mm and $t$=0. It is evident that the $\epsilon$ is weak in the regions of high density and low temperature (i.e, inside the density inhomogeneties) or low density and high temperatures (i.e, outside the density inhomogeneities). However the emission is maximum where temperatures are $\sim$1.6 MK (i.e, at the boundaries of the inhomogeneities). This happens because the $G(n_{e},T)$ sharply varies with temperature and peaks at $\sim$1.6 MK. It should be noted that after performing the LOS integration of $\epsilon$ in the $y-z$ plane, inhomogenities still appear brighter than the surroundings, even though the emission is less inside the inhomogeneties. The bottom right panel of Figure 4 shows the emission at $x$=20 Mm and $t$=45 when the turbulence has developed. This lead to the formation of fine scale structures in the synthetic images obtained after LOS integration of $\epsilon$.

We choose twelve different LOS directions (0*∘, 15∘, 30∘, 45∘, 60∘, 75∘, 90∘, 105∘, 120∘, 135∘, 150∘, and 165∘) through the $y-z$ plane perpendicular to the $x$ axis. After the application of the FoMo, LOS integration reduces the data cube to two spatial dimensions and provides the monochromatic specific intensity, $I_{\lambda}$ from which we calculate the total intensity, $I$, line widths ($\sigma$), and Doppler shifts ($\lambda_{D}$). As an example, the direction of integration for LOS=0∘* is along the $y$ axis (as shown in Figure 2 (c)) and we obtain $I_{\lambda}(x,z,t$). Then, at all instances and locations we compute, $I$, $\sigma$, $\lambda_{D}$ by taking the moments of $I_{\lambda}(,x,z,t$), as given by the following relations:

[TABLE]

Here $\lambda_{0}$ is the location of the peak of the Fe XIII emission line which is 10749 Å. Further, exponential line width ($\sigma_{1/e}=\sqrt{2}\sigma$) is computed and subsequently, converted to the velocity. Similarly the Doppler shifts are also converted to velocities ($v_{D}$). A similar method is adopted for deriving specific intensities, Doppler velocities and exponential line widths for other LOSs.

The synthetic images of the total intensity, Doppler velocity, and exponential line width for LOS=0*∘* and $U_{o}$=$\frac{11}{\sqrt{2}}$ km s*-1* are shown in Figure 5. An animation of this figure is available. While performing the FoMo, we choose a spectral resolution of 4 km s*-1* and the spatial resolution is kept as that of the simulation cube. To compare the synthetic observations with the observations, we degrade the spatial resolution of the synthetic observations to that of CoMP ($\sim$3200 km). However, for simplicity, we maintain the same spectral resolution, which is higher than the spectral resolution of the CoMP ($\sim 33$ km/s).222We mention that McIntosh & De Pontieu (2012) did not find any significant differences in the $v_{rms}$ and $\sigma_{1/e}$ when synthesizing data with both low and high spectral resolutions.

5 Analysis and Results

In the left panel of Figure 6 we show an example of the Fe XIII emission line profile (blue) at the start of the simulation ($t$=0), when the velocity driver at the bottom boundary is not applied. The width of the initial line profile shown is $\sim$22.5 km s*-1*, which is slightly greater than the thermal width of the Fe XIII line ($\sim$21.78 km s*-1*). This difference is due to the temperature and density inhomogeneities along the LOS. Overplotted in green is the line profile averaged over 400 s (20 frames)333We chose 20 frames for averaging the spectra because the average time period of oscillations in our simulation is 400 s (or 20 frames)., where there is additional broadening due to the time-averaged behaviour of the transverse wave motions.

The right panel of Figure 6 presents the line profile at the same location but with a larger amplitude driver ($U_{o}$=$\frac{22}{\sqrt{2}}$ km s*-1*). The width of the initial line profiles (blue) are the same in both panels, while the time-averaged line profile (green) is significantly broader than for the simulation with the lower amplitude driver. The larger wave amplitude produces a broader emission line profile when integrated over the time period of oscillations, provided the medium is optically thin, because of the superposition of the shifted spectra due to the incoherent and spatially unresolved swaying motions of the structures.

5.1 Superposition of line profiles

Since the solar corona is optically thin, the spectrum observed at a given location in the plane-of-sky (POS) is the superposition of the spectra of different structures oscillating with random phases and different polarisations along the LOS. Given the relatively short length of our simulation box, we effectively stack multiple realisations of the simulation along the LOS to mimic the corona. The following procedure is apodted to obtain random segments (from the synthetic observations) that are oscillating in random phases and with random polarisations.

First, we randomly choose several LOSs from the twelve different LOS (as described in Section 4) with a uniform probability of choosing any LOS. This is equivalent to choosing different polarisations of the waves. Next, we randomly chose several start times ($t_{0}$) of oscillations and considered twenty consecutive frames for every start time. The random choices of the start time are equivalent to assuming that different structures are oscillating with different phases. A segment, $I_{\lambda}(x,z^{\prime},t_{0,j})$, is then defined as a 2D projection of the simulation cube, obtained by choosing a random start time ($t_{0,j}$) of oscillation and integrated over a random LOS ($LOS_{j}$) direction. $z^{\prime}$ represents the arbitrary direction that is perpendicular to the chosen random LOS and confined in the $y-z$ plane. For example, if $LOS_{j}$=0*∘*, $z^{\prime}_{j}=z$ because the LOS is along $y$ axis. Since we have several LOSs and start times, we get several such random segments. Finally, we superimpose the intensities of all such random segments to obtain a resultant intensity that will be used for further analysis. This method can be understood with the following relation:

[TABLE]

where $I_{\lambda}(,x,z^{\prime},t_{n})$ is the resultant specific intensity obtained after the superposition of the intensities of 100 random segments. $LOS_{j}$ (thus $z^{\prime}_{j}$) is chosen randomly from twelve LOSs; $t_{0,j}$ is the starting time (in frame number) chosen randomly from a uniform distribution of [1,30]. $\delta t$ is the time cadence, kept as 1 frame (or 20 s) in this study. Since $n$ varies from 1-20 frames, the maximum allowed value of the random start time cannot be greater than 30. Following this procedure, the total number of unique segments available in this study are 360 (30$\times$12). One should note, if we increase $n$, the total number of unique segments will decrease because the number of allowed t0,j will decrease. On the other hand, if $n$ is small, very few frames will be available to compute the mean and root mean square (rms) estimates of the non-thermal line widths and Doppler velocities, respectively.

5.2 Measured wave properties vs number of segments

Next, we compute the moments of $I_{\lambda}(,x,z^{\prime},t_{n})$ using Equation 6 and derive the Doppler shifts and thus the Doppler velocity, $v_{D}(x,z^{\prime},t_{n})$, and exponential line width, $\sigma_{1/e}(x,z^{\prime},t_{n})$, at each instant ($n$) for twenty consecutive frames ($n$=1 – 20). Then, the non-thermal line widths ($\sigma_{nt}$) are estimated by subtracting quadratically the line width obtained at $t$=0 (22.5 km s*-1*) when no velocity drivers are applied (see Section 4) from the exponential line width. Finally, the rms Doppler velocities, $rms~{}v_{D}(x,z^{\prime})$, and mean non-thermal line widths, $mean~{}\sigma_{nt}(x,z^{\prime})$, are computed over 20 frames using the following relations:

[TABLE]

It should be noted that the rms Doppler velocities obtained after random sampling of segments is different from the rms velocity of the driver ($v_{rms}$) at the bottom boundary, which is described by Equation 4. Additionally, we degrade the spatial solution of the simulation cube to the spatial resolution of the CoMP ($\sim$3200 km) and a similar procedure as outlined above was adopted to estimate rms Doppler velocities and mean non-thermal line widths.

Given that a segment is obtained by integrating the simulation cube along a LOS perpendicular to the direction of the magnetic field ($x$ axis), a segment is thus integrated at least over $\sim$5 Mm. Therefore, a total of 100 segments are equivalent to a distance of at least 500 Mm, which corresponds to an inclination ($\tan^{-1}(250/R_{sun})$) of 20*∘* with respect to the normal to the surface of the Sun. Since the inclination is not large, we do not choose a LOS inclined to the $x$ axis of the simulation cube. This allows us to investigate the variation of the Doppler velocities and line widths with heights without mixing the emission coming from different heights.

Furthermore, we also note that the rms Doppler velocities and mean line widths averaged over $y-z$ plane sharply decrease and increase, respectively, with an increase in the number of segments, although both reach a plateau when the number of random segments is greater than 100 (see Figure 7). Therefore, we choose 100 segments in the current study. This result is partially in agreement with the findings of McIntosh & De Pontieu (2012), where the authors reported that the rms Doppler velocities monotonically decrease while mean line widths stay constant with increasing number of ‘threads’ (see their Figure 3). This difference may arise due to different methods employed in McIntosh & De Pontieu (2012) and our study. Furthermore, it should be noted that the properties of a ‘segment’ used in this study are different from a ‘thread’ used in McIntosh & De Pontieu (2012). They define a ‘thread’ as an elementary oscillating structure which is uniformly bright and optically thin, while a ‘segment’ used in this study is obtained by integrating the transversely inhomogeneous and gravitationally stratified simulation cube over 5 Mm along a given LOS. Therefore one segment comprises several ‘thread’-like structures along the LOS but driven uniformly at the bottom boundary. It is worth noting at this point that we do not give any selective weights to the intensities of 100 segments lying along the LOS. In reality, the foreground segments may scatter the photons coming from the background segments. Therefore, the emission from foreground structures might dominate over background structures. This effect is ignored here for the simplicity.

5.3 Doppler velocity vs non-thermal line width

To examine relationships between $rms~{}v_{D}(x,z^{\prime})$ and $mean~{}\sigma_{nt}(x,z^{\prime})$ (Equation 8), we create density and scatter plots comparing the two variables for all values of $x$ and $z^{\prime}$. Figure 8 (a) displays a two dimensional density plot of the variation of the rms Doppler velocities and mean non-thermal line widths for the simulation with driver amplitude of $U_{o}$=$\frac{11}{\sqrt{2}}$ km s*-1*. A wedge shape correlation between rms Doppler velocities and mean non-thermal line widths is evident and is qualitatively similar to Figure 1. Figure 8 (b) shows the variation of rms Doppler velocities with mean non-thermal line widths segregating different height ranges (shown in different colors) for the same driver amplitude. The mean non-thermal line widths increase with increasing heights, again resembling the observations taken from the CoMP (Figure 1). Figure 8 (c) is obtained after degrading the spatial resolution of the simulation cube to the spatial resolution of the CoMP, which is about 65% of the total extent of the simulation cube.

In Figure 8 (c), several Lissajous like curves for a given $x$ can be seen. This is the result of the degradation of the simulation cube. The simulation cube extends up to 5 Mm in $y$ and $z$ directions. Since the periodic boundaries were used in $y$ and $z$ directions while performing simulations, we perform the degradation of the simulation cube by wrapping the $y$ and $z$ boundaries of the simulation cube. This lead to the closed curves (Lissajous curves) in Figure 8 (c) for a particular $x$. We confirmed this effect disappears if, instead of wrapping the cube in $y$ and $z$ directions, nearest-neighbours is used to perform degradation; the curves shown in different colors are no longer closed but open. In addition to this, the degradation of the simulation cube filtered the small spatial variations in the $rms~{}v_{D}(x,z^{\prime})$ and $mean~{}\sigma_{nt}(x,z^{\prime})$ and only the large spatial variations that are very small compared to the one sigma standard deviations, remained. Such large scale variations in $rms~{}v_{D}(x,z^{\prime})$ and $mean~{}\sigma_{nt}(x,z^{\prime})$ (due to large scale density variation) for a given $x$ lead to the observed open or closed curves. However, like the real observations, if noise is included, it may distort the observed Lissajous curves and the scatter plots may look more realistic. Such Lissajous curves do not affect the results of this study and therefore, a wedge-shape correlation can still be noted in Figure 8 (c).

The results displayed in Figure 8 (a–c) demonstrate that the mean non-thermal line widths vary from $\sim$15 km s*-1* - 24 km s*-1*. Though the simulation generates a wedge-shaped correlation (as observed in Figure 1 and McIntosh & De Pontieu, 2012), the range of the line widths obtained for $v_{rms}$ $\sim$15 km s*-1* is narrower than the observed values of non-thermal line widths (which can reach up to 40 km s*-1*). Since the non-thermal broadening of emission lines is due to the LOS superposition of different transversely oscillating segments, it is prudent to assume that line widths must depend on the velocity amplitude of these oscillations (see also right panel of Figure 6). To study the effect of the amplitude of velocity drivers on mean line widths, we perform another simulation with $U_{o}$=$\frac{22}{\sqrt{2}}$ km s*-1* that gives a rms velocity of 26 km s*-1* at the bottom boundary (obtained using Equation 4), and the results are shown in Figure 8 (d–f). It can be noted from the figure that the mean non-thermal line widths increase and so does the rms Doppler velocities. This constitutes additional evidence of the correlation between rms Doppler velocities and non-thermal line widths. Moreover, a wedge-shaped correlation similar to Figure 8 (a–c) is also clearly seen. Similarly, we perform another simulation with $U_{o}$=$\frac{5}{\sqrt{2}}$ km s*-1* ($v_{rms}$ $\sim$7 km s*-1*) and found that both the mean non-thermal line widths and rms Doppler velocities decrease (see bottom panels of Figure 8).

Figure 9 presents the variation of the rms Doppler velocities and mean non-thermal line widths for different wave amplitudes of velocity drivers. We also choose random segments from the simulations with different velocity amplitudes as shown in red and grey colors in Figure 9. The wedge shape correlation can be conspicuously noted. It can be understood from this figure that the non-thermal line widths (and rms Doppler velocities) depend on the input wave amplitudes. Our results are in agreement with McIntosh & De Pontieu (2012), where the authors have shown that non-thermal broadening of emission line depends on the input wave amplitudes. Further, in our study, the wave amplitudes depend on the density, which decays exponentially with height. Thus transverse wave amplitude increases with height in the corona, and so do the non-thermal line widths. Therefore, the wedge-shape is at least partially due to height-dependent wave amplitude of the transverse MHD waves. Overall, our simulations could reproduce the height-dependent wedge-shape correlation between rms Doppler velocities and mean non-thermal line widths as seen in the observations without artificially adding any additional non-thermal broadening as done in McIntosh & De Pontieu (2012) to match the observed non-thermal line widths. These results allow us to conclude that large non-thermal line widths (due to the unresolved wave amplitudes) in solar corona conceal large wave amplitudes (and hence energies). This relaxes the requirement for artificially adding an extra unknown source of ‘dark’ or ‘hidden’ energy in the solar corona to match the non-thermal line widths.

5.4 Non-thermal widths vs height

Next, we study the variation of the non-thermal line widths with height above the solar photosphere. Figure 10 demonstrates the increase in the mean non-thermal line widths ($mean~{}\sigma_{nt}(x,z)$, averaged over $z$ axis for 100 random segments) with height due to the density stratification. We also show this is true for different LOSs and strengths of the velocity drivers ($U_{0}$). Curves in different colors represent the three different velocity drivers. To understand the effects of the gravitational stratification on the line widths, we perform a similar analysis as described above on a simulation (with $U_{o}$=$\frac{11}{\sqrt{2}}$ km s*-1*) without including gravity. We find that the non-thermal line widths do not vary with height in that case (see the curve in black in Figure 10). Further we note that the non-thermal line widths level-off with increasing heights (Figure 10). The nature of variation is consistent with the observations of line widths (Hassler et al., 1990; Banerjee et al., 1998; Doyle et al., 1998; Banerjee et al., 2009; Hahn et al., 2012) in coronal holes and a few other MHD simulations (AWSoM; Oran et al., 2017). The levelling-off of the mean line widths could be a signature of wave damping or reflections or both (Doyle et al., 1998; Hahn et al., 2012). However, the physics behind this plateau is beyond the scope of this paper and will be the subject of a future study.

5.5 Energy estimate

Finally we estimate the total energy density, $E$, in the simulation cube averaged over the entire duration and volume of the simulation cube. Total energy, $E$ is computed using the following relations.

[TABLE]

[TABLE]

Here, $V$ represents the volume of the simulation cube and $\phi$ is the gravitational potential. We compare $E$ derived using equations 9 and 10 with the time averaged Alfvén wave energy density $\langle\rho~{}rms~{}v_{D}^{2}\rangle_{x,z^{\prime}}$, estimated using the observed rms Doppler velocities of the synthetic images obtained by random superposition of 100 segments. This expression is frequently used to estimate the energy carried by the transverse and largely incompressible (Alfvénic) waves using the observed values of rms Doppler velocity fluctuations (Tomczyk et al., 2007; McIntosh et al., 2011). This comparison allow us to estimate the amount of the underestimation of the true wave energies due to the LOS superposition.

The average total energy ($E$) of the simulations with $U_{0}$ = $\frac{5}{\sqrt{2}}$ km s*-1*, $\frac{11}{\sqrt{2}}$ km s*-1*, and $\frac{22}{\sqrt{2}}$ km s*-1* is found to be $\sim$ 1.23$\times$10${}^{-5}~{}J~{}m^{-3}$, $\sim$ 4.2$\times$10${}^{-5}~{}J~{}m^{-3}$, and $\sim$ 2.2$\times$10${}^{-4}~{}J~{}m^{-3}$, respectively. The corresponding ‘observed’ average Alfvén wave energy is $\sim$ 9$\times$10${}^{-8}~{}J~{}m^{-3}$, $\sim$ 4.1$\times$10${}^{-7}~{}J~{}m^{-3}$, and $\sim$ 2.7$\times$10${}^{-6}~{}J~{}m^{-3}$. Therefore, the observed Alfvén wave energy was found to be $\sim$0.7%, $\sim$0.9%, and $\sim$1%, respectively, of the actual wave energy. This is much less than 10%-40% reported in De Moortel & Pascoe (2012) where the authors considered only ten loops along the LOS integration.

We also estimate the time average energy flux, $F$, injected into the simulation domain through the bottom boundary ($x$=0). It can be estimated by the following relation.

[TABLE]

where $v_{g}$ is the group speed of the transverse MHD wave which is $\sim$800 km s*-1* and $A$ is the area of the bottom boundary. $F$ at the bottom boundary for the simulations with $U_{0}$ = $\frac{5}{\sqrt{2}}$ km s*-1*, $\frac{11}{\sqrt{2}}$ km s*-1*, and $\frac{22}{\sqrt{2}}$ km s*-1* is found to be $\sim$22 Wm*-2*, $\sim$61 Wm*-2*, and $\sim$262 Wm*-2* respectively.

In addition to this, we also estimate the Poynting flux ($S_{x}=\frac{-1}{\mu_{0}}|({\bf v\times B)\times B}|_{x}$) passing through the bottom boundary of the simulations. Since the simulations are ideal, we have ignored the magnetic diffusivity. However, a small numerical diffusivity will be present in the simulations. $S_{x}$ for the simulations with $U_{0}$ = $\frac{5}{\sqrt{2}}$ km s*-1*, $\frac{11}{\sqrt{2}}$ km s*-1*, and $\frac{22}{\sqrt{2}}$ km s*-1* is found to be $\sim$14 Wm*-2*, $\sim$51 Wm*-2*, and $\sim$213 Wm*-2* respectively. Its worth noting that the values of $S_{x}$ are in good agreement with the flux computed using the equation 11.

We compare the energy flux computed using equation 11 and $S_{x}$ with those estimated using synthetic images obtained by random superposition of 100 segments. It can be calculated by multiplying the time averaged Alfvén wave energy density at the bottom boundary ($\langle\rho~{}rms~{}v_{D}^{2}\rangle_{x,z^{\prime}}$) with the Alfvén speed. The time averaged observed Alfvén wave energy flux is estimated to be $\sim$0.04 Wm*-2*, $\sim$0.21 Wm*-2*, and $\sim$1.08 Wm*-2* for simulations with $U_{0}$ = $\frac{5}{\sqrt{2}}$ km s*-1*, $\frac{11}{\sqrt{2}}$ km s*-1*, and $\frac{22}{\sqrt{2}}$ km s*-1* respectively. Again, we note that the observed Alfvén wave energy is $\sim$0.2% - 0.4%, respectively, of the true wave energy flux.

Thus our simulations are able to generate large differences of 2-3 orders of magnitude in the synthetically observed and the actual wave energy flux, which had probably tempted McIntosh & De Pontieu (2012) to coin the term ‘hidden’ or ‘dark’ energy.

6 Discussion and conclusions

A recent debate on the ‘hidden’ or ‘dark’ energy in the solar corona has arisen due to the discrepancy in the wave energies measured in corona using CoMP in comparison to those measured with SDO/AIA. To understand this discrepancy, we have performed a 3D MHD simulation of transverse MHD waves propagating in a gravitationally stratified plasma with a typical setting of coronal holes (open-field regions).

Observational studies using data from the SDO and CoMP have revealed that the polar coronal holes comprise several fine-scale magnetic structures that exhibit periodic transverse motions (Thurgood et al., 2014; Morton et al., 2015). Hence, our simulations utilize a series of density inhomogeneities which they are randomly positioned perpendicular to the direction of the magnetic field. Observations of Alfvénic waves in these regions revealed a log-normal distribution of periods with a mean value $\sim$ 470 km s*-1* and a velocity distribution with a mean value $\sim$ 14 km s*-1*, that can range from $\sim$5 km s*-1* - 40 km s*-1*. In order to model the observed properties of Alfvénic waves in coronal holes, a log-normal distribution of period with a mean value of 400 s and velocity drivers with rms wave amplitude of 7 km s*-1*, 15 km s*-1*, and 26 km s*-1* (with varying phase) are used.

We find that the forward modelling of Fe XIII 10749 Å emission line leads to a sharp decrease (and increase) in the rms Doppler velocities (and mean line widths) as we integrate over an increasing number of density inhomgenetities. These results are in close agreement with the findings of McIntosh & De Pontieu (2012). For this study, we perform MHD simulations of propagating transverse MHD wave with three different values of $v_{rms}$ only. More realistic simulations representing the true distribution of wave velocity amplitudes described in, e.g. Morton et al. (2015, 2019), can be performed in the future. Our simulations could generate a wedge-shaped correlation between the rms Doppler velocities and mean non-thermal line widths. We degrade the spatial resolution in our simulations to match the spatial resolution of CoMP but found no significant differences in the rms Doppler velocities and mean non-thermal line widths. In our study, transverse MHD waves and induced uniturbulence are able to produce the observed values of large non-thermal line widths due to LOS integration without including any significant contribution from flows or torsional motions as suggested by McIntosh & De Pontieu (2012); De Moortel & Pascoe (2012). Further, unlike the method described in McIntosh & De Pontieu (2012), we have not artificially added any extra non-thermal broadening to match the observed values of non-thermal line widths (15 km s*-1* - 45 km s*-1*). The MHD model used in this study is more detailed because it employs the mechanism of 3D MHD wave propagation, compared to the Monte Carlo model of oscillating ‘threads’ in McIntosh & De Pontieu (2012).

De Moortel & Pascoe (2012) have also studied the effects of LOS integrations by employing a 3D MHD simulation. However, the authors have not performed the forward modeling for emission lines and a single pulse of velocity driver was used. In this study, we use a, more realistic, continuous multiple periodic driver with random direction of the polarisation at the bottom boundary to mimic the photospheric driving of the plumes.

There are a number of caveats that come with our results. As seen by comparing Figure 1 (and previous studies, Tomczyk et al., 2007; McIntosh & De Pontieu, 2012; Morton et al., 2015) to Figure 9, the rms Doppler velocities calculated from our synthetic emission lines are larger than those observed in the CoMP data. One possible reason could be that the CoMP instrument has a finite exposure time. In our simulations, a snapshot is recorded instantaneously (exposure time $\sim$ 0 s). The finite exposure time leads to the averaging of spectra in time, thus reducing the magnitude of measured Doppler velocity fluctuations (as occurs from spatial averaging). Unfortunately we are unable to test this with the current simulations as the available unique segments are limited by the total time of the simulation (which is in turn limited by the computational resources and time).

Moreover, although we can reproduce the wedge-shape of the joint distribution of rms Doppler Velocity and non-thermal lines widths, it should be borne in mind that here we give a proof of concept that focuses on transverse waves propagating in open field regions. The exact shape of the correlation will depend on the conditions of the solar plasma at the time of observations which is of significantly greater complexity to model and well beyond the scope of this study.

The propagating waves in this study suffers a weak damping by resonant absorption because the simulation domain (50 Mm) is smaller than the damping length ($\sim$ 200 Mm) which is comparable to the wavelength of the transverse (kink) wave (Pascoe et al., 2010). This is further confirmed by performing a simulation without any gravitational stratification, where line widths did not decrease with height. Magyar et al. (2017) have also reported that the rms wave amplitude and Alfvénic wave energy flux in their MHD simulation without the inclusion of gravity were weakly damped over 50 Mm. Furthermore, we find that the non-thermal line widths first increase and then level off with the height in our simulations. Hahn et al. (2012) investigated several different mechanisms such as the effect of scattered light, photoexcitation of the emission lines, and inhomogeneities in the temperatures to explain the levelling-off and decrease in the non-thermal line widths. These authors finally proposed that the damping of waves can explain the observed decrease in the line widths. In our study, we rule out significant wave damping due to resonant absorption, due to numerical viscosity, and resistivity because if these were the dominant wave damping mechanisms then a decrease in the non-thermal line widths must have been seen in the simulation without the gravitational stratification. For a similar reason, we rule out that a substantial energy of the waves was spent in ohmic heating of the plasma (also Magyar et al., 2017). Moreover the non-thermal line widths in the non-stratified MHD simulation increase by 1-2 km s*-1* with increasing height. It allows us to conclude that the induced uniturbulence in these simulations does not play a significant role in the non-thermal broadening. Unresolved wave amplitudes are the main reason for the non-thermal broadening of the optically thin emission lines. Apart from the damping, reflections of wave due to density stratification is also expected in our simulation. The reflection of Alfvén wave energy due to gravitational stratification leading to the levelling-off non-thermal line widths is reported in Oran et al. (2017). The detailed mechanism of wave reflection and its effect on the variation of the non-thermal line widths in our simulations is beyond the focus of the current study and will be explored in the future.

Our study indicates that a spectrograph with good spectral resolution will not be able to resolve the Doppler velocities due to LOS superposition of the emission spectra of different structures, which is the result of the optically thin nature of the solar corona. However, a good spatial resolution may be useful to resolve the POS motions (similar to SDO/AIA). CoMP suffers from both coarse spatial resolution and the LOS superposition of structures in the optically thin corona. Therefore CoMP could not resolve the true wave amplitudes in the LOS Doppler velocity fluctuations and POS velocity fluctuations. This lead to the gross underestimation of the true wave energy flux. Our study reveals that only 0.2-1% of true wave energy flux or energy density is estimated by the resolved Doppler velocity fluctuations. The unresolved wave amplitudes results in the non-thermal broadening. Therefore, we conclude that Doppler velocities are not the true representatives of the wave amplitudes and wave energies. True wave energies are hidden (in form of unresolved wave amplitudes) in the non-thermal line widths. Our study has revealed the sites of ‘hidden’ or ‘dark’ wave energies without adding any unknown source of artificial energy. Finally, to estimate the energy budget and to explain the heating of the solar corona, a relation between true wave amplitudes and the non-thermal line width is needed and will be explored in the future.

Appendix A Photoexcitation of Fe XIII 10749 Å

The electron transitions for the Fe XIII line become dominated by photoionisation at a certain height above the solar photosphere. The specific intensity of a spectral line depends on the rate of the collisional excitations and radiative excitations (or photoionisation). The collisional excitation rate varies as $n_{e}^{2}$, while radiative excitation rate varies as $n_{e}$ (Young et al., 2003; Landi et al., 2016). Therefore, in the low corona, emission of the Fe XIII is dominated by the collisional excitations. While at the greater heights, the transitions are dominated by radiative excitations. Usually the contribution of the radiative excitations are low for typically use coronal emission lines, which are in the EUV, but for the infrared Fe XIII their contribution becomes significant at distances as low as 1.1 $R_{\odot}$ (Landi et al., 2016; Del Zanna & DeLuca, 2018).

In order to examine whether this influences our results, the forward modeling was also performed including the contribution of photoionisation where the source is set at a height of $1R_{\odot}$. The photoionisation is incorporated by including the rates of radiative excitation in addition to collisional excitations while computing the ionic fraction of the Fe XIII emission line and $G(n_{e},T)$. We assume that photoionisation happens due to a background radiation field emitted by a blackbody at $T$=5700 K. The detailed mechanism of computing $G$ including photoionisation can be found in Young et al. (2003). The results presented in the proceeding sections are identical whether photoionisation is included or not. This is likely due to our simulations extended up to 1.07 $R_{\odot}$, which is less than one pressure scale height. Hence, the electron density does not become so small that photoionisation dominates over collisional excitations. Thus, throughout the manuscript, we have not included the effects of photoionisation in our calculations of synthetic radiation from the Fe XIII emission line.

We thank the anonymous referee for useful suggestions that have improved the manuscript. TVD and VP were supported by the GOA-2015-014 (KU Leuven) and the European Research Council (ERC) under the European Union’s Horizon 2020 research and innovation programme (grant agreement No 724326). We also acknowledge Dr. Hui Tian and Prof. Valery Nakariakov for insightful discussions. The authors acknowledge the work of the National Center for Atmospheric Research/High Altitude Observatory CoMP instrument team.