Can the solar p-modes contribute to the high-frequency transverse oscillations of spicules?

We perform a two-dimensional numerical simulation that seamlessly covers the upper part of the solar convection zone and the corona. To this end, we use RAMENS¹¹1RAdiation Magnetohydrodynamics Extensive Numerical Solver code (Iijima & Yokoyama, 2015; Iijima, 2016; Iijima & Yokoyama, 2017; Wang et al., 2021; Kuniyoshi et al., 2023), in which we solve the compressible magnetohydrodynamic equations with gravity, radiation, and thermal conduction. The basic equations are given in the conservation form as follows.

where $\rho$ is the mass density, $\boldsymbol{v}$ is the gas velocity, $\boldsymbol{B}$ is the magnetic field, $e=e_{\rm int}+\rho\boldsymbol{v}^{2}/2+\boldsymbol{B}^{2}/8\pi$ is the total energy density, $e_{\mathrm{int}}$ is the internal energy density, $p$ is the gas pressure, $\boldsymbol{g}$ is the gravitational acceleration, and $\boldsymbol{\underline{I}}$ is unit tensor. $Q_{\mathrm{cnd}}$ and $Q_{\mathrm{rad}}$ denote the heating by thermal conduction and radiation, respectively.

The radiation $Q_{\rm rad}$ is determined through a combination of optically thick and thin components using a bridging law (Iijima, 2016). In calculating the optically-thick radiation, the frequency-averaged (i.e., grey-approximated) radiative transfer is directly solved under the local thermodynamic equilibrium (LTE) approximation, using the Rosseland mean opacity obtained from the OPAL opacity (Iglesias & Rogers, 1996). The optically-thin radiation is calculated from the loss function retrieved from the CHIANTI atomic database ver. 7.1, assuming the coronal abundance (Dere et al., 1997; Landi et al., 2012). Since the loss function from CHIANTI is defined in $T\geq 10^{4}$ K, we employ the loss function from Goodman & Judge (2012) in the lower-temperature range ($T\leq 10^{4}$ K) and smoothly connect the two functions using a bridging law (see Iijima, 2016, for detail). The equation of state is calculated under the LTE assumption, taking into account the six most abundant elements in the solar atmosphere (H, He, C, N, O, Ne). For specific details regarding the abundance and states of each element, refer to Iijima (2016). A significant portion of the internal energy near the solar surface is derived from the latent heat generated by changes in the internal states of atoms and molecules. The balance between this latent heat and radiative transfer plays a crucial role in accurately reproducing the granulation (Stein & Nordlund, 1998). The field-aligned thermal conduction of a fully-ionized plasma (Spitzer & Härm, 1953) is employed to calculate $Q_{\rm cnd}$. Although the assumption of full ionization is invalid in the chromosphere, it does not significantly influence the simulation because the thermal conduction in the chromosphere is minor. The detailed numerical procedure is found in Iijima (2016).

Letting $x$-axis be horizontal and $z$-axis be vertical, the simulation domain covers a spatial extent of $L_{x}=20\ \rm Mm$ in the $x$ direction and encompasses a vertical range from $2\ \rm Mm$ below the surface to $14\ \rm Mm$ above it, resulting in a total range of $L_{z}=16\ \rm Mm$ in the $z$ direction. $z=0\ \rm Mm$ is defined where the horizontally averaged optical depth is unity. The grid size is uniformly set to $25\ \rm km$ in $x$ direction and $25\ \rm km$ in $z$ direction. The periodic boundary conditions are applied in $x$ direction. We consider the loop-aligned simulation domain that extends from the upper convection zone to the top of the coronal loop, corresponding to one half of a symmetric closed loop. It consists of multiple flux tubes, associated with the network magnetic fields (Gabriel, 1976). Following Matsumoto (2016), we assume a half circle loop model for the gravitational acceleration as follows:

The bottom boundary condition is open for flow, mimicking the convective energy transport from the deep convection zone (see Iijima (2016) for detail). To ensure complete reflection of the Poynting flux at the top boundary (top of the coronal loop), we apply a reflective boundary that sets $v_{z}$ and $B_{x}$ to zero, while the other variables take on the same values as those one grid below. The reflected Poynting flux corresponds to the one injected from the other side of the loop. To sustain the coronal temperature, we impose artificial heating at the top boundary so that the temperature $T$ is fixed to $1.0\times 10^{6}$ K at the top. This artificial heating does not violate the scope of this work because our interest is in the transverse dynamics of spicules, not the amount of heating in the corona.

We conduct the simulation in two stages: the first stage is performed without considering p-modes (without p-mode stage), and the second stage incorporates the p-modes (with p-mode stage). In the without p-mode stage, the initial ($t=0$) condition in the convection zone is given by Model S (Christensen-Dalsgaard et al., 1996). Above the surface, the initial condition is calculated by the isothermal stratification. A uniform vertical magnetic field with a strength of $10\ \rm G$ is initially imposed. After $3\ \mathrm{hours}$ of integration, the convection is relaxed to a quasi-steady state, in which the enthalpy flux injected from the bottom boundary nearly equals the radiative flux. Following this, we perform an additional one hour of integration and utilize the numerical data obtained during this period $(10,800\ {\rm s}<t<14,400\ {\rm s})$.

We use the snapshot of the without p-mode stage at $t=14,400\ \rm s$ as the initial condition for the with p-mode stage. In this stage, we modify the bottom boundary condition to account for the influence of the p-modes because the depth of the convection zone in our simulation is not sufficient for the p-modes to develop (Finley et al., 2022). We introduce longitudinal waves into the simulation by adding vertical velocity ($\delta v_{z}$), density ($\delta\rho$), and pressure perturbations ($\delta p$) to the local values of $v_{z}$, $\rho$, and $p$ in addition to the convection motions. They are defined as follows:

where $C_{s}=\sqrt{\gamma p/\rho}$ is the sound speed and $\gamma=5/3$ corresponds to the specific heat ratio of adiabatic gas. $P_{\rm p}=300\ \rm s$ and $\lambda_{\rm p}=4.0\ \rm Mm$, which are the typical values for the period and horizontal wavelength of p-modes obtained by previous observations (Leighton et al., 1962; Christensen-Dalsgaard, 2002; Oba et al., 2017; McClure et al., 2019). The longitudinal waves produced by the p-mode-like driver propagate upward and reach the photosphere. Following the observation of p-modes by Oba et al. (2017), $v_{\rm p}$ in Equation (6) is a constant value ($=4.7\times 10^{4}\ \rm cm\ s^{-1}$) designed so that the power of $v_{z}$ at $z=0\ \rm Mm$ in the without and the with p-mode stage correlate as

Figure 1 illustrates the velocity power spectra, denoted as $E_{j}^{\rm fq}$ and $E_{j}^{\rm wn}$, as a function of frequency and wavenumber at $z=0\ \rm Mm$. Here, $j=x,z$, and these quantities are defined as follows:

where the indices for the time ($t$) and horizontal ($x$) directions are denoted as $n_{t}$ and $n_{x}$, respectively, with lengths $N_{t}$ and $N_{x}$. We sampled frequencies ($\nu=\omega/2\pi$) and horizontal wavenumbers ($k_{x}$) in increments of $\Delta\omega=2\pi/(N_{t}\Delta t)=1.7\times 10^{-3}\ \rm s^{-1}$ and $\Delta k_{x}=2\pi/(N_{x}\Delta x)=3.1\times 10^{-9}\ \rm cm^{-1}$, where $\Delta x$ represents the grid size in the $x$ direction (which is $25\ \rm km$), and $\Delta t$ is equal to $2\ \rm s$. The power spectra exhibit sensitivity to the bottom boundary condition (i.e., with or without p-mode stages). The power spectra for the with p-mode stage exhibit distinct peaks at locations that correspond to the typical values of the p-modes, $\nu=P_{\rm p}^{-1}$ and $k_{x}=2\pi\lambda_{\rm p}^{-1}$.

Figure 2 displays the snapshots of the magnetic field and temperature distribution in both the without and with p-mode stages. After the relaxation of the convection, expanding magnetic flux tubes are formed of which footpoints are concentrated into kilogauss magnetic fields. In addition, Figure 2 reveals several spicule-like structures consisting of plasmas with chromospheric temperatures ($T\sim 10^{4}\ \rm K$). We utilize a fully automatic method with three steps to detect spicules and measure the period and velocity amplitude of the transverse oscillations in them. The method consists of the following steps:

1. 1.

   Identifying the spicules within the simulated data.

2. 2.

   Tracking the transverse oscillations of the identified spicules at a specific height.

3. 3.

   Deriving periods and velocity amplitudes of the spicule oscillations.

Step 1: We utilize an automated approach to identify spicules, considering them as local minima of temperature below $40,000\ \rm K$ across both temporal ($t$) and horizontal ($x$) dimensions at a specific height. Note that only local minima with widths in the $x$-direction less than $1,100\ \rm km$ are categorized as spicules. This criterion aligns with the maximum spicule width observed in the quiet Sun atmosphere (Pasachoff et al., 2009). The lifetime of each spicule is determined by identifying the start and end times of the local minimum. In the left panel of Figure 3, a temperature map is displayed, emphasizing the edges of the spicule by white lines.

Step 2: At the given height and time, we calculate the mean horizontal ($v_{x}$) and perpendicular velocities (to the local magnetic field, $v_{\perp}$) between the right and left edges of the spicules. We then proceed to record the temporal evolution of them. In the right panel of Figure 3, we present an example of the recorded $v_{x}$ along the displacement, demonstrating an oscillatory pattern.

Step 3: To get the oscillation period, we utilize the fast Fourier transform (FFT) on the recorded temporal variation and generate the power spectra, denoted as $E^{\rm fq}_{\rm sp}$, with respect to frequency $\nu$ defined as,

Figure 4 illustrates the distributions of the $\langle S_{z}\rangle$ for the stages without and with p-modes, where

In this section, we investigate whether high-frequency spicule oscillations exhibit any signatures of energy transfer resulting from the longitudinal waves initiated by the p-mode-like driver, as shown in Section 3.1. We have applied our spicule oscillation analysis method at various heights, focusing on regions where the number of oscillations ($N_{\rm os}$) is at least 100. This analysis encompasses heights up to $5,400\ \rm km$, where $N_{\rm os}\geq 100$ in both the without and with stages, and the lowest height considered corresponds to the lower limit of spicule height observed in the quiet Sun ($4,200\ \rm km$, Pasachoff et al., 2009). As an illustrative example, Figure 5 displays the density distributions of velocity amplitudes $A_{v}$ and oscillation periods $P_{v}$ at $z=4.5\ \rm Mm$ obtained from $v_{\perp}$ within the spicules. These distributions are estimated using kernel density estimation (KDE) with a Gaussian kernel. The bandwidth parameter is selected through cross-validation³³3Conducted using Scikit-learn (Pedregosa et al., 2011), following the same method as in Morton et al. (2021). The distributions of $P_{v}$ reveal a dominance of high-frequency (periods $<100\ \rm s$) components in both stages with a peak around $20\ \rm s$. Moreover, the distributions of $P_{v}$ demonstrate minimal difference between the two stages, while that of $A_{v}$ indicates an enhancement in the with p-mode stage compared to the without p-mode stage. These are common characteristics across all considered heights. The distributions of $A_{v}$ and $P_{v}$ obtained from $v_{x}$ and $v_{\perp}$ within the spicules exhibit a striking degree of similarity (see the distributions derived from $v_{x}$ as depicted in Figure 6 in the Appendix). This consistency arises primarily because a lot of spicules in our simulation are vertical due to our top boundary condition, which enforces a vertical magnetic field, i.e., $B_{x}=0\ \rm G$, although some of them do exhibit clear inclinations.

At each of the heights investigated, we have computed the mean values of the oscillation period $\bar{P}_{v}$ and velocity amplitude $\bar{A}_{v}$. Table 1 summarizes the oscillation parameters ($N_{\rm os}$, $\bar{P}_{v}$, and $\bar{A}_{v}$) derived from perpendicular velocities $v_{\perp}$ within the detected spicules. Across different altitudes, there is no significant variation in mean parameters, and this remains within the range of standard errors. It is worth noting that the mean parameters derived from $v_{x}$ closely resemble those obtained from $v_{\perp}$ (refer to Table 3 in the Appendix). Furthermore, to assess the impact of the p-mode-like driver on the mean velocity amplitude, we have calculated the ratio of the mean velocity amplitude between the with and without p-mode stage, denoted as $R_{v,\rm sp}=\bar{A}_{v,\rm with}/\bar{A}_{v,\rm without}$. The presence of the driver leads to an intensification of the mean velocity amplitudes by $10$–$30\%$ (i.e., $R_{v,\rm sp}=1.1$–$1.3$) when using $\bar{A}_{v}$ from $v_{\perp}$ and $10$–$40\%$ (i.e., $R_{v,\rm sp}=1.1$–$1.4$) when using $\bar{A}_{v}$ from $v_{x}$. When we assess the increase in the power of velocity amplitude $R_{v,\rm sp}^{2}$, we obtain $R_{v,\rm sp}^{2}=1.3$–$1.8$ when using $v_{\perp}$ and $R_{v,\rm sp}^{2}=1.2$–$1.9$ when using $v_{x}$. This intensification exceeds the power of horizontal velocity at $z=0\ \rm Mm$ ($30\%$, see Equation (10)). Hence, this result is unlikely to be solely attributed to the increase in horizontal velocity caused by the p-mode-like driver. Instead, it is more likely due to the energy transfer from the longitudinal wave energy generated by the driver to transverse wave energy, as explained in Section 3.1.