Can the Solar p-modes Contribute to the High-frequency Transverse Oscillations of Spicules?

Lateral motions of spicules serve as vital indicators of transverse waves in the solar atmosphere, and their study is crucial for understanding the wave-heating process of the corona. Recent observations have focused on high-frequency transverse waves (periods < 100 s), which have the potential to transport sufficient energy for coronal heating. These high-frequency spicule oscillations are distinct from granular motions, which have much longer timescales of 5–10 minutes. Instead, it is proposed that they are generated through the mode conversion from high-frequency longitudinal waves that arise from a shock-steepening process. Therefore, these oscillations may not solely be produced by the horizontal buffeting motions of granulation but also by the leakage of p-mode oscillations. To investigate the contribution of p-modes, our study employs a two-dimensional magneto-convection simulation spanning from the upper convection zone to the corona. During the course of the simulation, we introduce a p-mode-like driver at the bottom boundary. We reveal a notable increase in the mean velocity amplitude of the transverse oscillations in spicules, ranging from 10%–30%, and attribute this to the energy transfer from longitudinal to transverse waves. This effect results in an enhancement of the estimated energy flux by 30%–80%.


INTRODUCTION
The temperature of the solar corona surpasses 10 6 K, exceeding the surface temperature by several hundredfold (Grotrian 1939;Edlén 1943).The underlying cause for the high-temperature corona is believed to reside within the solar magnetic field.However, the precise mechanisms responsible for this magnetically-driven heating continue to be a topic of active scientific debate, referred to as the "coronal heating problem" (Klimchuk 2006(Klimchuk , 2015;;Reale 2014;De Moortel & Browning 2015;Cranmer & Winebarger 2019;Van Doorsselaere et al. 2020;De Pontieu et al. 2022).
Solving the coronal heating problem involves three key steps (Klimchuk 2006): understanding how mechanical energy is transferred to the corona, how energy dissipates within the corona, and how the corona thermally responds to heating events.This research aims to improve our understanding of the first step, energy transfer.Although energy transport by flux emergence possibly plays a significant role (Cranmer & van Ballegooijen 2010;Wang 2020Wang , 2022)), transverse fluctuations (or Alfvénic waves, Alfvén 1942;McIntosh et al. 2011), have been considered a promising energy carrier due to their highly efficient nature of field-aligned energy transport.Various observations and simulations investigate the heating scenario by the transverse waves (see the recent reviews by Van Doorsselaere et al. 2020;Morton et al. 2023) but its quantitative contribution is still under investigation.
The observations by the Coronal Multi-channel Polarimeter (CoMP, Tomczyk et al. 2008) have revealed that the power spectra of coronal transverse waves exhibit a distinct peak around 4 mHz.This feature is globally present in terms of both spatial and temporal characteristics (Tomczyk et al. 2007;Tomczyk & McIntosh 2009;Morton et al. 2015Morton et al. , 2019)).The peak frequency corresponds to the typical time scale of the solar p-modes, which are acoustic eigen oscillations in the convection zone (Leighton et al. 1962;Ulrich 1970;Deubner 1975).These results suggest that the horizontal motion of magnetic flux tubes due to granular buffeting (Spruit 1981;Steiner et al. 1998;Choudhuri et al. 1993;Musielak & Ulmschneider 2002;Fujimura & Tsuneta 2009) is not only a source of coronal transverse waves but also that the p-modes may contribute to their generation by transferring a fraction of longitudinal wave energy into transverse wave energy.The transfer mechanism is still under debate, with several proposed theories.These include the direct excitation process of transverse waves by the interaction between flux tubes and longitudinal waves below the equipartition layer (so called mode absorption, Bogdan et al. 1996;Hindman & Jain 2008;Riedl et al. 2019;Skirvin et al. 2023b), as well as the mode conversion process from longitudinal to transverse waves at the equipartition layer where the local sound speed and Alfvén speed become equal (Schunker & Cally 2006;Cally 2007;Jess et al. 2012;Wang et al. 2021;Shimizu et al. 2022).Additionally, another mode conversion from fast to Alfvén waves around the refraction height of the fast waves, which is inherently a three-dimensional process unlike the mode conversion at the equipartition layer and the mode absorption (Cally & Goossens 2008;Cally & Hansen 2011;Hansen & Cally 2012;Khomenko & Cally 2012).
Transverse oscillations of plasmas in coronal loops without obvious damping, which are called decayless oscillations, are widely used for proxies of transverse waves in the corona (Wang et al. 2012;Tian et al. 2012;Nisticò et al. 2013).Recent high-cadence observations conducted by the Extreme Ultraviolet Imager (EUI, Rochus et al. 2020) aboard the Solar Orbiter (SolO, Müller et al. 2020) have identified the existence of high-frequency decayless oscillations (periods < 100 s) in both active regions and the quiet Sun regions (Zhong et al. 2022;Petrova et al. 2023;Mandal et al. 2022;Shrivastav et al. 2023).Furthermore, Lim et al. (2023) have analyzed the observed decayless oscillations and revealed that highfrequency transverse waves play a more dominant role in coronal heating than low-frequency (periods > 100 s) transverse waves.
Spicules, which are characterized as chromospheric jets embedded in the corona (see the reviews by, e.g., Beckers 1968Beckers , 1972;;Sterling 2000;Tsiropoula et al. 2012;Skirvin et al. 2023a), frequently exhibit transverse oscillations, used as proxies for the transverse waves (see the recent review by Jess et al. 2023, and references therein).Several observations have provided evidence that high-frequency spicule oscillations, of which typical period ranging from 20 to 50 s, carry substantial energy flux, contributing to both chromospheric and coronal heating (Okamoto & De Pontieu 2011;Srivastava et al. 2017;Bate et al. 2022).These results support the significance of high-frequency transverse waves for coronal heating.
High-frequency spicule oscillations cannot be generated by horizontal granular motions due to their significantly longer time scale (300-1000 s, Schrijver et al. 1997).Longitudinal waves are considered one of the potential mechanisms for generating these oscillations, as demonstrated in previous numerical studies through processes such as mode conversion (Shoda & Yokoyama 2018) or mode absorption (Gao et al. 2023;Skirvin et al. 2023b).Therefore, p-modes may play a role in generating high-frequency spicule oscillations.However, the contribution of p-modes relative to granulations remains uncertain, because both p-modes and vertical granular motions generate longitudinal waves and observations cannot differentiate their contributions in terms of oscillation period.Hence, our study aims to explore the potential role of p-modes in generating high-frequency spicule oscillations using a two-dimensional magnetoconvection simulation.

Simulation Setup
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 RA-MENS1 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 ρ is the mass density, v is the gas velocity, B is the magnetic field, e = e int + ρv 2 /2 + B 2 /8π is the total energy density, e int is the internal energy density, p is the gas pressure, g is the gravitational acceleration, and I is unit tensor.Q cnd and Q rad denote the heating by thermal conduction and radiation, respectively.The radiation Q 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 CHI-ANTI is defined in T ≥ 10 4 K, we employ the loss function from Goodman & Judge (2012) in the lowertemperature range (T ≤ 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 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 Mm in the x direction and encompasses a vertical range from 2 Mm below the surface to 14 Mm above it, resulting in a total range of L z = 16 Mm in the z direction.z = 0 Mm is defined where the horizontally averaged optical depth is unity.The grid size is uniformly set to 25 km in x direction and 25 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: where g = 2.74 × 10 4 cm s −2 , R sun = 6.96 × 10 10 cm, θ = z/r, h = r sin θ, r = 2L z /π, and ẑ is the unit vector in the z-direction.
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 × 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 pmode 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 G is initially imposed.After 3 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 s < t < 14, 400 s).
We use the snapshot of the without p-mode stage at t = 14, 400 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 (δv z ), density (δρ), and pressure perturbations (δp) to the local values of v z , ρ, and p in addition to the convection motions.They are defined as follows: where C s = γp/ρ is the sound speed and γ = 5/3 corresponds to the specific heat ratio of adiabatic gas.P p = 300 s and λ p = 4.0 Mm, which are the typical values for the period and horizontal wavelength of p-modes obtained by previous observations (Leighton et al. 1962; Following the observation of p-modes by Oba et al. (2017), v p in Equation ( 6) is a constant value (= 4.7 × 10 4 cm s −1 ) designed so that the power of v z at z = 0 Mm in the without and the with p-mode stage correlate as where the operator ⟨⟩ denotes the temporal and horizontal averaging.While Oba et al. (2017) did not observe the horizontal velocity field, for reference, we present the ratio of the power of v x at z = 0 Mm as follows, It is worth noting that v x is amplified due to the pmode-like driver despite the driver not being applied to v x because the longitudinal waves generated by the driver propagate isotropically through the convection zone to the photosphere.The simulation was run for 4 hours (14, 400 s < t < 28, 800 s), with the last hour (25, 200 s < t < 28, 800 s) specifically dedicated to analyzing the numerical data.It is worth noting that the temporal averages of Equation ( 9) and ( 10) are calculated over the duration 25, 200 s < t < 28, 800 s. Figure 1 illustrates the velocity power spectra, denoted as E fq j and E wn j , as a function of frequency and wavenumber at z = 0 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 (ν = ω/2π) and horizontal wavenumbers (k x ) in increments of ∆ω = 2π/(N t ∆t) = 1.7 × 10 −3 s −1 and ∆k x = 2π/(N x ∆x) = 3.1×10 −9 cm −1 , where ∆x represents the grid size in the x direction (which is 25 km), and ∆t is equal to 2 s.The power spectra exhibit sensitivity to the bottom boundary condition (i.e., with or without p-mode stages).The

Spicule Oscillation Analysis Method
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 spiculelike structures consisting of plasmas with chromospheric temperatures (T ∼ 10 4 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. Identifying the spicules within the simulated data.
2. Tracking the transverse oscillations of the identified spicules at a specific height.
3. Deriving periods and velocity amplitudes of the spicule oscillations.
In constructing our three-step method, we draw upon the NUWT2 code as a reference (Morton et al. 2013;Thurgood et al. 2014;Weberg et al. 2018).
Step 1 : We utilize an automated approach to identify spicules, considering them as local minima of temperature below 40, 000 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 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 ⊥ ) 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 fq sp , with respect to frequency ν defined as, where k = x, ⊥, v k,sp denotes the recorded velocity of the spicule during its lifetime, N k,sp represents the sample size of v k,sp for each spicule.We select the frequency to be that of the significant wave component, fulfilling the following condition: Note that the frequency components which show at least 3/4 of an oscillatory cycle are chosen.This approach leads to the selective detection of the high-frequency oscillations because it cannot identify oscillations with periods longer than the typical lifetime of spicules, which is on the order of a few minutes.The individual oscillation often exhibits multiple significant frequency components.We fit the time evolution of the horizontal velocity oscillations with a sinusoidal and linear function v fit (t) as v fit (t) = where the subscript i expresses the individual wave component of the superposition, N mode represents the number of oscillation modes chosen according to Equation (15), A v,i is the velocity amplitude, P i is the period (inverse of the selected frequency).We obtain A v,i , ϕ i , c 1 , and c 2 through the fitting procedure.

Energy Transfer between Waves
Figure 4 illustrates the distributions of the ⟨S z ⟩ for the stages without and with p-modes, where ⟨S z ⟩ in the with p-mode stage is enhanced above z = 1 Mm compared to that in the without p-mode stage.
In contrast, the increase in energy flux transported by longitudinal waves in the with p-mode stage initiates below z = 0 Mm as a consequence of the p-mode-like driver.(see Section 2.1).Therefore, the increase in ⟨S z ⟩ signifies the transfer of energy from the longitudinal waves triggered by the driver to the transverse waves.It's worth noting that ⟨S z ⟩ in the with p-mode stage below z=0.5 Mm is negative because of the greater amplification of ⟨S emerge . Due to the temporal and horizontal averaging applied to the energy fluxes, the amplification of ⟨S z ⟩ likely involves multiple mechanisms, such as the mode conversion or mode absorption, which are local phenomena.

Statistical Properties of Spicule Oscillations
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 os ) is at least 100.This analysis encompasses heights up to 5, 400 km, where N os ≥ 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 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 Mm obtained from v ⊥ within the spicules.These distributions are estimated using kernel density estimation (KDE) with a Gaussian kernel.The bandwidth parameter is selected through cross-validation3 , following the same method as in Morton et al. (2021).The distributions of P v reveal a dominance of high-frequency (periods < 100 s) components in both stages with a peak around 20 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 ⊥ 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 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 Pv and velocity amplitude Āv .Table 1 summarizes the oscillation parameters (N os , Pv , and Āv ) derived from perpendicular velocities v ⊥ within the detected spicules.Across different altitudes, there is no significant vari-ation 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 ⊥ (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,sp = Āv,with / Āv,without .The presence of the driver leads to an intensification of the mean velocity amplitudes by 10-30% (i.e., R v,sp = 1.1-1.3)when using Āv from v ⊥ and 10-40% (i.e., R v,sp = 1.1-1.4)when using Āv from v x .When we assess the increase in the power of velocity amplitude R 2 v,sp , we obtain R 2 v,sp = 1.3-1.8when using v ⊥ and R 2 v,sp = 1.2-1.9when using v x .This intensification exceeds the power of horizontal velocity at z = 0 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.

DISCUSSION
The velocity amplitude of high-frequency spicule oscillations is enhanced by the presence of the p-mode-like driver.Therefore, despite accounting for the refraction of transverse (fast magnetic) waves, p-modes still have the potential to contribute to the production of highfrequency spicule oscillations, although the exact mechanism of the energy transfer cannot be definitively determined.
The distributions of oscillation periods for both stages exhibit a peak around 20 s.This value is consistent with the mean crossing time of acoustic waves across the equipartition layer, denoted as ⟨τ eq ⟩, which is expressed as follows (Shoda & Yokoyama 2018): ⟨τ eq ⟩ = 21 s in the without and with p-mode stage, as indicated by the dashed lines in the top panels of Figure 5.This result can be explained as follows: pulse-like amplifications in the velocity amplitude, lasting for approximately 20 s, are generated through the mode conversion at the equipartition layer (Schunker & Cally 2006;Cally 2007), as shown in Shoda & Yokoyama (2018).This finding is a consistent quantitative feature across all the heights we examined.However, it is important to note that this factor does not conclusively rule out the presence of other energy transfer systems, such as the mode absorption (Bogdan et al. 1996;Hindman & Jain 2008;Riedl et al. 2019;Skirvin et al. 2023b).
Many observations employ the velocity amplitude of spicule oscillations to estimate the carried energy flux F (Van Doorsselaere et al. 2014;Jess et al. 2023), using the following formula: v ⊥ is obtained from the mean velocity amplitude of spicule transverse oscillations observed, while ρ and v A are derived from various models (e.g., a bright network chromospheric model by Vernazza et al. 1981).Therefore, the squared value of the velocity amplitude serves as a primary indicator of the energy flux propagating within spicules.When we assume that all the obtained values of Āv,without and Āv,with represent the velocity amplitudes of propagating waves, the energy flux carried by them through spicules is amplified by 30-80% (i.e., R 2 v,sp = 1.3-1.8).To analyze the nature of the detected oscillations, i.e., whether they are propagating or standing waves, we calculate their phase velocity v ph as follows: where ∆ϕ represents the phase lag of the individual waves ϕ at two different heights (5, 400 km and 4, 800 km), as obtained from Equation ( 15), and d represents the height difference, i.e., d = 600 km.ϕ and P v are derived from v ⊥ within the spicules.To detect features of the same wave at different heights, we derive the properties of individual waves at a specific height and search for similar ones at an adjacent height, following the criteria of Bate et al. (2022).The considered properties (and criteria) include: (i) the x-position of the spicule averaged over its lifetime (with a difference not exceeding 250 km, i.e., 10 grids), (ii) the lifetime of the spicule (ensuring the duration of its existence between the considered heights overlaps), (iii) the duration of the oscillation (with a difference not exceeding 50%), and (iv) the oscillation period (with a difference not exceeding 10%).Based on these criteria, we have identified 13 upwardly propagating, 7 downwardly propagating, and 6 standing waves in the without pmode stage.In the with p-mode stage, we have identified 6 upwardly propagating, 2 downwardly propagating, and 10 standing waves.Waves with a phase speed v ph > 500 km s −1 are classified as standing waves, following the same criteria applied by Okamoto & De Pontieu (2011).It should be noted that the number of the upward Āv downward Āv standing Āv driver (km s −1 ) (km s −1 ) (km s −1 ) without 3.7 ± 0.5 3.7 ± 0.4 4.0 ± 0.9 with 5.9 ± 0.8 6.9 ± 1.2 3.8 ± 0.8 identified waves is much smaller than the total number of the detected oscillations N os .This is because the transverse waves in our two-dimensional simulation are fast mode waves, which propagate not only along but also across the magnetic field lines inside the spicules.Considering the obtained values of Āv,without and Āv,with in our simulation (averaged between those at z = 5, 400 km and z = 4, 800 km) presented in Table 2, the energy flux propagating through spicules is intensified by around 150% (R 2 v,sp = 2.5±1.4) for the upward waves and 250% for the downward waves (R 2 v,sp = 3.5±1.9).On the other hand, there is no clear increase for standing waves (i.e., R 2 v,sp = 0.9 ± 0.8).Thus, propagating waves within spicules are significantly enhanced by the p-mode-like driver, while standing waves are not affected.However, care must be taken in the classification of the standing waves.Upwardly and downwardly propagating waves are superposed in the temporal evolutions of the spicule oscillations obtained by our method, potentially leading to misidentification of standing waves.The more sophisticated classification, achieved by filtering upwardly and downwardly propagating waves using spatiotemporal Fourier analysis (Tomczyk & McIntosh 2009), may change the result.
The oscillation periods obtained from horizontal velocities in simulated spicules (Figure 5) fall within the observed range of high-frequency spicule oscillations, as indicated by previous studies (Okamoto & De Pontieu 2011;Bate et al. 2022).The velocity amplitudes of the simulated spicule oscillations (Table 1) are smaller than many observed values (median ∼ 10 km s −1 , Okamoto & De Pontieu 2011; Pereira et al. 2012;Bate et al. 2022), while a few studies show good correspondence (Jafarzadeh et al. 2017;Yoshida et al. 2019).This discrepancy is likely to come from the top boundary condition in our simulation.In our simulation, the top boundary condition is implemented such that the magnetic field becomes vertical, i.e., B x = 0 G, which consequently suppresses the transverse velocities of spicules.
The top boundary condition (B x = 0 G) also results in the reflection of transverse waves.To compare the en-ergy flux carried by the reflected wave at the top boundary and the incident wave at the transition region, we have estimated them by using Equation ( 18).We verify that ⟨F ⟩ at the upper boundary is 5% of that at the transition region in the absence of the p-mode-like driver, and 3% in the presence of the driver.Hence, we can disregard the influence of the top boundary on the transverse spicule oscillations (i.e., below the transition region).However, the wave modulation in the corona arising from this reflection is not to be disregarded.Consequently, our numerical model does not have the capacity to investigate the generation of the 4 mHz peak observed in the coronal power spectra, as demonstrated in CoMP (Tomczyk et al. 2007;Tomczyk & McIntosh 2009;Morton et al. 2015Morton et al. , 2019)).
It is important to note that our two-dimensional simulation does not account for the three-dimensional effects such as the other mode conversion, responsible for generating Alfvén waves from fast waves (Cally & Goossens 2008;Cally & Hansen 2011;Hansen & Cally 2012;Khomenko & Cally 2012).In addition, spicules display three-dimensional motions, encompassing not only bulk transverse (kink-mode-like) oscillations but also rotational and cross-sectional oscillations (Sharma et al. 2017(Sharma et al. , 2018)).These additional factors may indeed have an impact on the properties of high-frequency spicule oscillations produced by the p-modes.Further investigations using three-dimensional simulations are necessary to explore these effects in detail.

CONCLUSION
Our study focuses on examining the potential role of solar p-modes in generating high-frequency transverse spicule oscillations.To investigate this, we utilize a twodimensional magneto-convection simulation that encompasses the upper convection zone to the corona.We introduce a driver in the middle of our simulation duration to generate longitudinal waves with periods and wavelengths resembling the typical values of p-modes.
We then proceed to analyze and compare the spicule oscillation characteristics between the stages without and with the p-mode-like driver.We find that the mean velocity amplitude increases by 10-30% by the driver, leading to the estimated enhancement of the energy flux in the spicules by 30-80%.This result implies that both p-modes and granulations play a substantial role in generating high-frequency spicule oscillations.Further investigations utilizing three-dimensional simulations are necessary to explore the contribution of p-modes to producing three-dimensional spicule transverse oscillations.
We would like to convey our sincere appreciation to the anonymous referee for providing valuable feedback.Numerical computations were carried out on the Cray XC50 at the Center for Computational Astrophysics (CfCA), National Astronomical Observatory of Japan.M.S. is supported by JSPS KAKENHI Grant Number JP22K14077.R.J.M. is supported by a UKRI Future Leader Fellowship (RiPSAW MR/T019891/1).H. K. is also grateful for travel support provided by the UKRI Future Leader Fellowship (RiPSAW MR/T019891/1).T.Y. is supported by the JSPS KAKENHI Grant Number JP21H01124, JP20KK0072, and JP21H04492.This work was supported by NAOJ Research Coordination Committee, NINS, Grant Number NAOJ-RCC-2301-0301.

APPENDIX
A. ADDITIONAL FIGURES AND TABLES Figure 6 shows the density distributions of velocity amplitudes and oscillation period of spicules similar to those shown in Figure 5, but obtained from horizontal velocities v x within the detected spicules.Table 3 presents the spicule oscillation parameters (N os , Pv , and Āv ) derived from horizontal velocities.

Figure 1 .
Figure 1.Power spectra of vx and vz versus the frequency and horizontal wavenumber at z = 0 Mm.The black dashed lines correspond to the typical frequency ν of P −1 p = 3.3 × 10 −3 Hz (top row) and wavenumber kx/(2π) of λ −1 p = 2.6 × 10 −9 cm −1 (bottom row) of the p-modes.Christensen-Dalsgaard 2002; Oba et al. 2017; McClure et al. 2019).The longitudinal waves produced by the pmode-like driver propagate upward and reach the photosphere.Following the observation of p-modes byOba et al. (2017), v p in Equation (6) is a constant value (= 4.7 × 10 4 cm s −1 ) designed so that the power of v z at z = 0 Mm in the without and the with p-mode stage correlate as

Figure 2 .
Figure2.The snapshots consist of the (left column) magnetic field strength maps in color with magnetic field lines in grey lines and the transition region defined by T = 10 4 K highlighted by black lines, and (right column) temperature maps in color with the transition region highlighted by white lines.Both columns include cases for both the without and with p-mode stages.An animation of this figure is available that shows the temporal evolutions of the magnetic field and temperature distributions and the positions of the transition region over a period of 1, 800 s.power spectra for the with p-mode stage exhibit distinct peaks at locations that correspond to the typical values of the p-modes, ν = P −1 p and k x = 2πλ −1 p .

Figure 3 .
Figure 3. (Left) A zoomed-in view of the temperature map, specifically focusing on a spicule.The white line indicates the transition region.The black dashed line shows the height corresponding to the right panel (z = 4.5 Mm).(Right) The mean horizontal velocity vx within a spicule is represented by the color of the dots, while the position of each dot corresponds to the mean horizontal coordinate (x) between the right and left edges during the lifetime of the spicule.The black dashed line shows the time corresponding to the left panel (t = 11, 010 s).

Figure 4 .
Figure 4.The height dependence of the averaged vertical Poynting flux (denoted as ⟨Sz⟩) during the without and with p-mode stage.

Figure 5 .
Figure 5.The kernel density estimation for the distributions of (top) the oscillation periods Pv and (bottom) the velocity amplitude Av derived from perpendicular velocities v ⊥ .Each panel displays parameters in the without (blue line) and with p-mode stage (orange line) with a 95% confidence interval (CI) for the estimates calculated using bootstrapping.The black dotted line in the top panel shows ⟨τeq⟩.

Figure 6 .
Figure 6.Same as Figure 5, but derived from horizontal velocities vx within the detected spicules.

Table 1 .
Spicule oscillation properties in the without and with p-mode stage derived from perpendicular velocities v ⊥ .without (mean ± standard error) with (mean ± standard error) height N os,without Pv,without Āv,without N os,with Pv,with

Table 2 .
Velocity amplitudes and phase velocities of the propagating or standing waves in the without and with pmode stage derived from perpendicular velocities v ⊥ within the spicules.

Table 3 .
Same as Table 1, but derived from horizontal velocities vx within the detected spicules.without (mean ± standard error) with (mean ± standard error) height N os,without Pv,without Āv,without N os,with Pv,with