Finding the mechanism of wave energy flux damping in solar pores using numerical simulations
Julia M. Riedl, Caitlin A. Gilchrist-Millar, Tom Van Doorsselaere, David B. Jess, Samuel D. T. Grant
AAstronomy & Astrophysics manuscript no. main © ESO 2021February 25, 2021
Finding the mechanism of wave energy flux damping in solar poresusing numerical simulations
J.M. Riedl , C.A. Gilchrist-Millar , T. Van Doorsselaere , D.B. Jess , and S.D.T. Grant Centre for mathematical Plasma Astrophysics (CmPA), KU LeuvenCelestijnenlaan 200B bus 2400, 3001 Leuven, Belgiume-mail: [email protected] Astrophysics Research Centre, School of Mathematics and Physics, Queen’s University BelfastBelfast, BT7 1NN, U.K.Received; accepted
ABSTRACT
Context.
Solar magnetic pores are, due to their concentrated magnetic fields, suitable guides for magnetoacoustic waves. Recentobservations have shown that propagating energy flux in pores is subject to strong damping with height; however, the reason is stillunclear.
Aims.
We investigate possible damping mechanisms numerically to explain the observations.
Methods.
We performed 2D numerical magnetohydrodynamic (MHD) simulations, starting from an equilibrium model of a singlepore inspired by the observed properties. Energy was inserted into the bottom of the domain via di ff erent vertical drivers with a periodof 30s. Simulations were performed with both ideal MHD and non-ideal e ff ects. Results.
While the analysis of the energy flux for ideal and non-ideal MHD simulations with a plane driver cannot reproduce theobserved damping, the numerically predicted damping for a localized driver closely corresponds with the observations. The strongdamping in simulations with localized driver was caused by two geometric e ff ects, geometric spreading due to diverging field linesand lateral wave leakage. Key words.
Waves – Methods: numerical – Sun: photosphere – Sun: oscillations – sunspots – Magnetohydrodynamics (MHD)
1. Introduction
Solar pores are macroscopic features resembling small sunspotslacking a penumbra, but can also occur as a precursor or remnantof sunspots (Garcia de La Rosa 1987; Sobotka 2003; Thomas& Weiss 2004). Given their nearly circular symmetry and highmagnetic field concentrations, pores act as e ffi cient wave guidesfor magnetoacoustic waves, allowing wave flux to enter higherregions of the solar atmosphere (Jess et al. 2015) where the en-ergy can then be dissipated (Grant et al. 2018).The observational evidence of waves in solar pores is vast.Photospheric sausage modes in pores were identified by, e.g.,Fujimura & Tsuneta (2009) (sausage and / or kink waves), Mor-ton et al. (2011) (being fast waves according to Moreels et al.2013), Dorotoviˇc et al. (2014) (standing slow and fast waves),Grant et al. (2015) (propagating slow surface waves), Keys et al.(2018) (surface and body waves), and Gilchrist-Millar et al.(2021) (propagating slow surface and body waves, hereafter re-ferred to as GM21). These authors all found evidence of waveperiods of around 3 and / or 5 minutes, indicating the likely roleof photospheric p-modes as a driver for the waves.The propagation to the chromosphere was studied byBalthasar et al. (2000), who confirmed, by using the VacuumTower Telescope (VTT) on Tenerife, the presence of five-minuteoscillations for the magnetic field in the deep photosphere, asseen in other observations. Using the Transition Region andCoronal Explorer (TRACE) observations, they found a peak ata period of three minutes in the chromosphere. Stangalini et al.(2011) reported longitudinal acoustic waves reaching the chro- mosphere in both three- and five-minute bands. They underlinethe strong connection between wave transmission and magneticfield geometry, which suggests that for pore models special at-tention should be paid to the definition of the magnetic field, asalso suggested by Jess et al. (2013).However, how far waves in solar pores propagate into higherlayers of the solar atmosphere is still unclear. Khomenko & Col-lados (2006) conducted numerical simulations of waves in asmall sunspot; they used a localized wave source to study wavepropagation, refraction, and mode conversion. They found thatdue to the vertical and horizontal stratification of the Alfvénspeed, (low β ) fast waves are refracted in the chromosphere backdown to the photosphere, while slow modes continue propagat-ing up. Recent simulations of a chromospheric resonance layerabove a sunspot done by Felipe et al. (2020) show that actualwave propagation only takes place between the photosphere andchromosphere. A chromospheric resonance layer was previouslyalso simulated by Botha et al. (2011) and observed by Jess et al.(2020). On the other hand, Riedl et al. (2019) showed in concen-trated flux tubes that plane waves are converted to tube (sausageand kink) waves that are able to propagate to the corona sincethe tube structure greatly a ff ects the wave propagation (Cally &Khomenko 2019; Khomenko & Cally 2019).Grant et al. (2015), and more recently GM21, measured waveenergy throughout the lower atmosphere of solar pores, and in-deed report significant energy flux damping as a function ofheight. Analytic calculations of Yu et al. (2017) show that theobserved damping could at least be partly explained by reso-nant damping of slow sausage waves. Although they find that Article number, page 1 of 12 a r X i v : . [ a s t r o - ph . S R ] F e b & A proofs: manuscript no. main this damping mechanism is stronger than previously expected,the numerical studies of Chen et al. (2018), validated by ana-lytic calculations of Geeraerts et al. (2020), show that dampingdue to electrical resisitvity is much more potent than that dueto resonant absorption. However, this alone is not enough to ac-count for the damping. Flux could also be lost due to leaky tubewaves (Cally 1986). Leaking waves had already been observedby Stangalini et al. (2011) and Morton et al. (2012). Grant et al.(2015) mentioned that part of the waves in their observationsare reflected, which fits the simulations of Khomenko & Colla-dos (2006), and that mode conversion (Cally 2001; Bogdan et al.2003) might play a role. Frequency dependent damping for slowmagnetoacousic waves in sunspot umbrae is discussed by Kr-ishna Prasad et al. (2017), who find that higher frequencies aredamped more strongly. The authors suspect this behavior occursdue to radiative and / or conductive losses.Another important factor to consider is the cuto ff frequencypresent in stratified media (Lamb 1909). Acoustic waves withlower frequencies than the cuto ff frequency cannot propagate,but are evanescent standing waves. This e ff ect can be used todetermine the cuto ff frequency of the solar atmosphere (Felipeet al. 2018), which indicates that five-minute waves like thoseobserved by GM21 should be below the cuto ff . However, thephase lag and propagation speed between di ff erent heights sug-gest that the observed waves in GM21 are indeed propagat-ing as evanescent waves should not show any phase di ff erences(Carlsson & Stein 1997). On the other hand, the picture is notcompletely black and white. Centeno et al. (2006) summarizedthe e ff ects of the cuto ff frequency; when radiative losses aretaken into account, they find that there is no clear value forthe cuto ff frequency, but a transition between mainly evanescentand mainly propagating waves. Therefore, it is possible that thewaves in GM21 are partly subject to the cuto ff , which could ac-count for at least part of the observed damping. For the sake ofthis study, however, we assume the waves to be 100% propagat-ing.In this paper we aim to expand our understanding of thedamping mechanisms in solar pores by explaining the observeddamping with simple two-dimensional (2D) numerical simula-tions, using a model inspired by the observational parametersobtained by GM21 for their pore 3. We insert propagating wavesat the bottom of the domain with a vertical velocity driver witha frequency above the cuto ff frequency, and study the result-ing wave energy flux with height in comparison with the datafrom GM21 for di ff erent setups. In Section 2 we briefly reiter-ate the most important points of GM21 before introducing themodel, the numerical setup, and the approach for calculating thewave energy flux. The results, distinguished by driver location,are presented in Section 3 and thoroughly discussed in Section4. A short discussion about the case of a driver with frequencylower than the cuto ff frequency is presented in Appendix A.
2. Methods
The model developed in this work is inspired by the observationsdetailed in GM21, who utilized data obtained by the Facility In-frared Spectropolarimeter (FIRS; Jaeggli et al. 2010) based atthe National Science Foundation’s Dunn Solar Telescope (DST),Sacramento Peak, New Mexico. The FIRS data consist of sit-and-stare slit-based spectropolarimetric observations of the de-caying active region NOAA 12564, which was captured between14:09 - 15:59 UT on 2016 July 12 in the Si I 10827 Å spectral
Fig. 1.
Energy flux across all five observed pores as a function of height.The color scale is logarithmic. Pore boundaries are shown as whitedashed lines. The green solid line shows the inclination angle of themagnetic field. From GM21. line. The observations acquired contain a set of five solar poresthat were positioned along a unique straight-line configuration.To cover all pores in a single FIRS exposure, the DST coude ta-ble was rotated so the spectrograph slit passed through the centerof each photospheric pore boundary.An examination of the photospheric Si I 10827 Å line bisec-tor velocities showed periods on the order of five minutes acrossall pore structures. Through spectropolarimetric inversions us-ing the Stokes Inversion based on Response functions (SIR; RuizCobo & del Toro Iniesta 1992) code, the local plasma densities,magnetic field strengths, and temperatures were deduced as afunction of atmospheric height spanning the range 0 – 500 km.The central pore (pore 3 in GM21) exhibited the best signal-to-noise ratio, and so was selected for comparison with the presenttheoretical work.For pore 3 documented by GM21, the magnetic fields werefound to be close to vertical toward the pore center, with fieldstrengths of 2400 G and 1000 G at atmospheric heights of 0km and 500 km, respectively. Temperatures ranged from 5000K to 3500 K and densities spanned from 8.5x10 − kg m − to9.8x10 − kg m − across the same height range. Parameters de-rived from the inversions were combined with mean square ve-locities to calculate energy flux estimates as a function of at-mospheric height (Equation 12). The energy flux across all fivepores as a function of height is displayed in Figure 1. Pore 3 wasfound to exhibit considerable energy damping with an averageenergy flux on the order of 25 kW m − at an atmospheric heightof 100 km, dropping to 1.5 kW m − at 500 km. The dampingmechanisms producing this drop in energy flux remain elusive.In addition, an increase in energy flux toward the boundaries ofpore 3 indicated the presence of surface mode waves. In order to investigate the wave damping in solar pores withnumerical simulations, we first need to create a 2D gravita-tionally stratified magnetohydrostatic (MHS) equilibrium atmo-sphere that resembles the observational data. For a MHS equlib-
Article number, page 2 of 12.M. Riedl et al.: Finding the mechanism of wave energy flux damping in solar pores using numerical simulations rium the following condition must be fulfilled ∇ p − µ ( ∇ × B ) × B − ρ g = , (1)where p is the gas pressure, ρ is the density, B is the magneticfield, µ is the magnetic permeability, and g is the gravitationalacceleration.We start by choosing the magnetic field in the z -direction B z ( x , z ) = a ( z ) (cid:34) arctan (cid:32) x + r ( z ) s ( z ) (cid:33) − arctan (cid:32) x − r ( z ) s ( z ) (cid:33)(cid:35) + b ( z ) , (2)which results in a bundle of strong vertical magnetic field of a ( z )inside the pore above the background field b ( z ) outside the pore,resembling the observations. The parameter r ( z ) describes the ra-dius of the pore, while s ( z ) is the smoothness parameter, whichdefines the thickness of the transition between pore and back-ground. For the sake of simplicity, the written dependence onthe vertical coordinate ( z ) is omitted from now on for these fourparameters.The parameters defining Equation 2 are a = . G axis [T] , r = × √ G axis [m] , b = . G side [T] , s = . r [m] , (3)with exponential functions approximating the observed magneticfield strength of pore 3 from GM21 at the axis of the pore G axis = . − z / + .
07 [T] and the side of the pore G side = . − z / + .
02 [T].Because div B =
0, we know that in 2D ∂ B x ∂ x = − ∂ B z ∂ z . (4)Therefore, B x ( x , z ) = − (cid:90) ∂ B z ∂ z dx + h ( z ) = − dadz (cid:34) s (cid:32) ( x − r ) + s ( x + r ) + s (cid:33) + ( x + r ) arctan (cid:18) x + rs (cid:19) − ( x − r ) arctan (cid:18) x − rs (cid:19)(cid:21) − a (cid:34) dsdz ln (cid:32) ( x − r ) + s ( x + r ) + s (cid:33) + drdz arctan (cid:18) x + rs (cid:19) + drdz arctan (cid:18) x − rs (cid:19)(cid:35) − dbdz x , (5)where we assume that h ( z ) = x = ∂ p ∂ x ∂ z = ∂ p ∂ z ∂ x (6)must be true. By di ff erentiating the x -component of Equation1 with respect to z and the z -component with respect to x , andcombining the resulting derivatives with Equation 6, we find aconstraint for the density, ∂ρ∂ x = µ g (cid:34) ∂ B x ∂ x (cid:32) ∂ B z ∂ x − ∂ B x ∂ z (cid:33) + B x (cid:32) ∂ B z ∂ x − ∂ B x ∂ z ∂ x (cid:33) + ∂ B z ∂ z (cid:32) ∂ B z ∂ x − ∂ B x ∂ z (cid:33) + B z (cid:32) ∂ B z ∂ x ∂ z − ∂ B x ∂ z (cid:33)(cid:35) , (7) assuming g = [0 , − g ] with g =
274 m / s. The density can then beobtained with ρ ( x , z ) = (cid:90) ∂ρ∂ x dx + f ( z ) . (8)The function f ( z ) is of great importance here as it defines thegravitational stratification of the density. We therefore set f ( z ) tobe equal to the average density obtained from the observationsof GM21 for their pore 3. Since ∂ρ/∂ x also has a dependence on z , we add a small constant to ρ to ensure its non-negativity. Dueto the complexity of Equation 8 it is solved numerically.From the second component of Equation 1, the pressure canbe calculated with p ( x , z ) = (cid:90) ∂ p ∂ z dz + j ( x ) . (9)As long as the pressure is symmetric around the pore axis at x =
0, there is no need to add a function j ( x ). However, we adda constant to ensure a positive pressure. This equation is alsosolved numerically.Theoretically, the described model is in MHS equilibrium.However, numerical calculations as used in the solution of themodel and in the simulation code itself are imperfect, often re-sulting in somewhat unstable behavior, especially when gravityis involved. Therefore, using the boundary conditions describedin Section 2.3, we simulate the model without driver for a phys-ical time of 1300 seconds to let it settle down. After this time,there are no significant changes to density, magnetic field, orpressure on a timescale compared to a few driver periods. Thisslightly relaxed atmosphere is then used as the initial conditionfor our simulations. We note, however, that even after the slightrelaxation there are still significant velocities within the domain,meaning that the resulting model atmosphere has not completelysettled to a MHS equilibrium.The top panel of Figure 2 shows the vertical magnetic fieldcomponent of the initial atmosphere, with field lines shown inorange. Due to the symmetry of the problem, only half of thepore is included in our model, with the pore axis being located at x =
0. The pore itself is located on the left side of the plot, wherethe magnetic field is strong and mainly vertical. The deviation ofthe horizontal profile from the arctan-shape of Equation 2 oc-curs because of the equilibration process. The comparison of themodel with the observations of GM21 (Figure 2 bottom) showsgreat similarity. It should be noted, however, that the horizontalextent of our model pore (FWHM radius ≈ .
44 Mm) is smallerthan the pores in the observations (radius ≈ . ff erent radii into account, thefield inclinations also coincide quite well (see Figure 3). Theplasma- β in our model is higher than unity everywhere, with val-ues ranging from 2 to 6.5 inside the pore and higher values up to40 and higher outside.Similarly, Figure 4 shows the density of the settled modeland the comparison to observations, where the density struc-tures seen appear during the equilibration process. It is imme-diately apparent, that the model density is much less stratifiedwith height than the observations, even though we added the ob-servational density as stratification in Equation 8. This is causedby the e ff ect of ∂ρ/∂ x calculated by Equation 7 already havinga dependence on z , which in total decreases the stratification. Inaddition, it is also slightly decreased when the atmosphere is al-lowed to settle down. However, it should be noted that for the Article number, page 3 of 12 & A proofs: manuscript no. main
Fig. 2.
Top:
Vertical component of the magnetic field of the settledmodel atmosphere. Orange lines depict the magnetic magnetic fieldlines.
Bottom:
Comparison of the model atmosphere with the observa-tions from GM21 of pore 3. The maximum observational value withinthe pore is shown by obs. peak , while the horizontal average across thewhole pore is shown by obs. average . The green lines show the modelvalues for the indicated lines in the top figure (line 1 at pore axis). observations in GM21, the density is not a direct output of theinversions, but is instead determined through solving equationsof state using inferred inversion outputs, under the assumptionof hydrostatic (HS) equilibrium. This simplifying assumption isproblematic in strong magnetic fields as it ignores the Lorentzforce, thus providing notable uncertainties on the densities inputinto the model, of up to an order of magnitude (Borrero et al.2019).Nonetheless, the density values from GM21 are still consis-tent with those from semi-empirical models like that of Maltbyet al. (1986), who considered a magnetized atmosphere at thecenter of a sunspot umbra. They also assumed a HS equilibrium;however, this assumption is valid for the center of an axiallysymmetric sunspot as the magnetic terms in Equation 1 vanish.Therefore, we have to assume that the observational values ofthe density are more reliable than the model values.The smaller pore radius and less stratified density in ourmodel compared to the observations are due to compromises be-ing made when solving Equation 1. Once a non-force-free mag-netic field is chosen, the density or pressure cannot be freelychosen, but only manipulated through the addition of integrationconstants, as can be seen in Equations 8 and 9. Therefore, inorder to obtain a stable model for our simulations, certain con-cessions have to be made. In addition, due to the same reasons,our model also results in a plasma- β > β found by GM21 within the pores. The impactof the di ff erences between observations and theoretical model onour results is discussed in Section 4.1. Fig. 3.
Comparison of the magnetic field inclination between model(red) and observations for pore 3 of GM21 (black) as a function of hor-izontal distance normalized to the pore radius. The pore radius for themodel was assumed to be 0.44 Mm, while the radius of the observedpore is 2.5 Mm. The observational values are taken along the line per-pendicular to the slit. Model values are mirrored around x = x = Fig. 4.
Top:
Logarithm of the density of the settled model atmosphere.Orange lines depict the magnetic magnetic field lines.
Bottom:
Compar-ison of the model atmosphere with the observations from GM21 of pore3. The maximum observational value within the pore is shown by obs.peak , while the horizontal average across the whole pore is shown by obs. average . The green lines show the model values for the indicatedlines in the top figure (line 1 at pore axis).
All our simulations are conducted using the PLUTO code(Mignone et al. 2007), which solves the magnetohydrodynamic(MHD) equations when using the respective module. Fluxes are
Article number, page 4 of 12.M. Riedl et al.: Finding the mechanism of wave energy flux damping in solar pores using numerical simulations computed using a linearized Roe Riemann solver, while the timestep is advanced using an unsplit second-order accurate charac-teristic tracing method, which is less dissipative than multi-stepalgorithms. To deal with the inevitable occurrence of div B weuse the mixed hyperbolic / parabolic divergence cleaning tech-nique of Dedner et al. (2002), which is further discussed inMignone et al. (2010). Gravity is added using a body force withconstant acceleration toward the negative z -direction.Keeping a model atmosphere stable when gravity is includedcan often prove di ffi cult and is highly contingent on the bound-ary conditions at boundaries perpendicular to the gravitationalacceleration. In our case it was not possible to set fully openboundary conditions at the top boundary. We therefore expandthe model atmosphere at the top to add a thick high-viscositylayer to absorb all outgoing waves, e ff ectively having an openboundary. We use the same boundary condition for the rightboundary. The viscosity is treated with an explicit time integra-tion. Including the viscous layers the domain ranges from 0 to3 Mm in the x -direction and from -0.095 to 0.795 Mm in the z -direction with 1000 ×
297 cells, leading to a spatial resolution of3 km in both directions. Excluding the viscous layers, a physicaldomain remains ranging from 0 to 2.68 Mm with 894 cells in the x -direction and from -0.095 to 0.475 Mm with 191 cells in the z -direction. Only this physical domain is used for the analysis andfigures. The height of 0 Mm is defined as the bottom of the pho-tosphere. After settling the calculated model from Section 2.2for 1300 s (defined as t = v z , driver = A sin (cid:32) π T t (cid:33) , (10)with the amplitude A =
160 m / s and the period T =
30 s. Sincethe driver purely perturbs the velocity, some of the driver energyimmediately flows into pressure and density perturbations. Be-tween the ghost cells (additional cells outside the computationdomain to enable numerical integration) including the driver andthe first cell of the domain, the root mean square of the velocityperturbation is therefore reduced to levels observed by GM21 atthe bottom of the pores of about 50 m / s. This short period for thedriver was chosen because a typical p-mode period of 300 s isclose to the cuto ff period in our model, leading to the formationof standing waves due to reflections. However, we want to in-vestigate propagating waves and their damping. In addition, forlonger periods the wavelength would increase accordingly, caus-ing the resulting waves to not fit into the domain. For the sakeof completeness, we also did simulations with a low-frequencydriver below the cuto ff frequency, and we show a crude analysisin Appendix A.For some of our simulations we include non-ideal e ff ects likeviscosity, resistivity, and thermal conduction. Those e ff ects wereadded using explicit time integration, and for expected values inthe photosphere ( R e and R m taken from Ossendrijver 2003). Inthe case of the simulations with viscosity, where viscosity wasalso present in the physical domain, simulations were done withexaggerated values for the viscous shear coe ffi cient. The energy flux can be calculated as (e.g., Goedbloed & Poedts2004) E = − µ ( v × B ) × B + (cid:32) ρ v + ρ Φ + γγ − p (cid:33) v , (11)where Φ = − gz + const . is the gravitational potential. The leftterm of Equation 11 is the Poynting flux, which is the magneticcomponent of the energy flux, whereas the other terms describethe hydrodynamic component. Since in our model β > / s within the wholephysical domain or up to ∼
350 m / s within the pore after settlingthe atmosphere for 1300 s. These velocities are higher than thedriver amplitude. Thus, in addition to simulations with a driver,we also conduct simulations without a driver, allowing us to ex-tract the e ff ects caused by the input waves alone. This is done bysubtracting all the variables of the simulations without a driverfrom the variables of the simulations with a driver, e ff ectivelygiving us the perturbed variables ρ (cid:48) = ρ driver − ρ nodriver , p (cid:48) = p driver − p nodriver , v (cid:48) = v driver − v nodriver , B (cid:48) = B driver − B nodriver . To obtain the wave energy flux, these perturbed variables are putinto Equation 11, in a process that is similar to linearization.In GM21, on the other hand, the wave energy flux was cal-culated as E = ρ v g (cid:104) v (cid:105) , (12)with v g being the group speed and (cid:104) v (cid:105) being the mean squarevelocity. For our simulations, Equations 11 and 12 yield simi-lar trends with absolute values in the same order of magnitude.Using Equation 11 facilitates a more detailed analysis, which ispossible due to the much more detailed knowledge of the data insimulations compared to observations.
3. Results
We conducted a range of simulations, including and removingnon-ideal e ff ects, and applying di ff ering drivers. Depending onthe driver location, the results can be divided into two distinctgroups, which are discussed in the following. We applied the velocity driver described in Equation 10 on thewhole bottom boundary, resulting in plane fast waves propagat-ing upward at approximately the sound speed. A single snapshotof the vertical velocity perturbation after two driver periods isshown in Figure 5. The amplitude of the waves increases withheight, as is expected due to the conservation of energy in astratified plasma. The wave fronts are not completely horizon-tal, but have a jagged form at the pore location. This happensdue to di ff ering wave speeds at di ff erent locations. The verticalwave ridges visible at x ≈ . Article number, page 5 of 12 & A proofs: manuscript no. main
Fig. 5.
Snapshot of the vertical velocity perturbation after two periodsfor the full driver and ideal MHD. The gray lines show magnetic fieldlines. The blue line highlights the field line considered for the analysisin Figure 6. The full time sequence is available as a movie online.
Fig. 6.
Relative wave energy flux parallel to the magnetic field averagedover time as a function of height for the full driver. The solid green lineshows the energy flux along the pore axis, whereas the dashed purpleline shows the average flux across the pore up to the field line high-lighted in Figure 5. The observational data (black line with symbols)are from pore 3 of GM21. All fluxes are normalized to the first obser-vational data point. observational data were normalized to the data point at z = . ff ects to a ff ect the waves. Therefore, we fail to re- Fig. 7.
Snapshot of the vertical velocity perturbation after two periodsfor the localized driver for ideal MHD. The gray lines show magneticfield lines. The red bar below the x -axis indicates the driver location.The blue line highlights the field line considered for the analysis in Fig-ure 8. The full time sequence is available as a movie online. produce the observed damping with a plain driver located at thewhole bottom boundary. Solar pores are magnetic structures that do not form in the pho-tosphere but are already present below the solar surface. As so-lar pores are good wave guides, it is valid to assume that onlythe pore itself may be driven. Numerical simulations (Cameronet al. 2007) supported by observations (Cho et al. 2013) sug-gest that rapid cooling within pores could lead to downflows thatcollide with the plasma of lower layers to produce reboundingupflows, which further motivates the assumption of a localizeddriver. Moreover, previous simulations (Kato et al. 2016) showthat photospheric bu ff eting by turbulent motions lead to the e ffi -cient excitation of waves. We therefore alter our driver to a step-function driver that is only present in the inner part of the pore v z , driver = A sin (cid:16) π T t (cid:17) x ≤ . x > . . Figure 7 shows a snapshot of the vertical velocity perturba-tion after two periods for the step-function driver. In contrast tothe respective figure for the full driver, the wave fronts are nothorizontal and the maximum amplitude is lower because the at-mosphere is only driven at one location (indicated by the redbar below the x -axis). The blue line highlights a field line rootedslightly outside the driver location at x = .
22 Mm. There areclearly waves present beyond this field line, suggesting that thewaves do not purely propagate along the magnetic field.If we now study the wave energy flux as a function of heightfor the simulation with localized driver, as shown in Figure 8(left), it is immediately apparent that the energy flux is nowstrongly damped, in stark contrast to the simulations with the fulldriver. This sudden drop in wave energy flux with height by justchanging the driver location can be explained by two geometricmechanisms: geometric spreading and lateral wave leakage.
The magnetic field lines in our model diverge with height. There-fore, if the waves were perfectly propagating along the field
Article number, page 6 of 12.M. Riedl et al.: Finding the mechanism of wave energy flux damping in solar pores using numerical simulations
Fig. 8.
Left:
Relative wave energy flux parallel to the magnetic field av-eraged over time as a function of height for the localized driver. Thesolid green line shows the energy flux along the pore axis, whereas thedashed purple line shows the average flux across the pore up to the fieldline highlighted in Figure 7. The observational data (black line withsymbols) are from pore 3 of GM21. All fluxes are normalized to thefirst observational data point.
Right:
Comparison of flux damping in thesimulation with localized driver with the e ff ects of geometric damping.The solid green line shows the same data as in the left plot for compar-ison. The other lines are explained in the text. lines, the flux along a single field line, as well as the average fluxat each height within the pore, would be expected to drop dueto the flux being distributed across a wider area with increasingheights. The decrease in flux with height due to this mechanismis proportional to 1 / R in 2D geometry, where R is the distancebetween the pore axis and a specific field line. Such a curve isshown in Figure 8 (right, dotted red line) for the field line high-lighted in Figure 7. Since this curve drops substantially less withheight than the wave flux, there must be another mechanism withapproximately equal significance.In addition, if only geometric spreading caused the damp-ing, the wave flux parallel to the magnetic field integrated acrossthe pore should be constant with height because the same totalamount of flux would be contained inside the pore at all heights.This is not the case, which can be seen with the dash-dotted or-ange line in Figure 8. Therefore, flux must escape from the porethrough its edges. In our simulations with a localized driver, we observe wavespropagating out of the solar pore, which decreases the flux in-side the pore. This is the case because magnetoacoustic wavescan propagate at an inclined angle with respect to the magneticfield. In a homogeneous plasma, the phase speed of fast and slowmagnetoacoustic waves is (e.g., Goedbloed et al. 2019) v fa / sl ( θ ) = (cid:113) v s + v A √ ± − v c cos θ v s + v A / / , (13)where v s is the sound speed, v A the Alfvén speed, v c = v A v s / ( v A + v s ) / the cusp speed, and θ the angle between the propaga-tion direction and the magnetic field. The positive (negative)sign is for the calculation of the phase speed of the fast (slow)wave. In a plasma where v s > v A (approximately β > v fa ( θ = , π ) = v s in the magnetic field directionand v fa ( θ = π/ , π/ = ( v A + v s ) / perpendicular to it. Onthe other hand, slow waves take the shape of double quasicir-cles with v sl ( θ = , π ) = v A in the magnetic field directionand v sl ( θ = π/ , π/ = ff ect was previously used by Nakariakov & Zi-movets (2011) to explain flare ribbon propagation.Assuming local homogeneity and utilizing Equation 13, wecan apply the Huygens-Fresnel principle to theoretically predictthe locations of fast and slow wave fronts. In order to do this weassume that the wave originates from a point source. In this pointthe phase speed in all directions is calculated, supplying us withinformation of the wave front location in the next snapshot. Forall subsequent snapshots we calculate the phase speed in eachpoint of the previous wave front. The next fast (slow) wave frontis then the outer edge of all fast-wave quasicircles (slow-wavedouble quasicircles).We let our theoretical wave fronts for both fast and slowwaves propagate from two point sources at the bottom of the do-main: one at the pore axis at ( x , z ) = (0 , − . x , z ) = (0 . , − . t = t = x ≈ . x ≈ . ff erence in densityleads to a di ff erence in phase speed.Although there are clearly waves leaking out of the pore,most of the flux is contained within the pore, following themagnetic field lines. To estimate the e ff ect of lateral leakage onthe damping of energy flux with height, we compare the time-integrated total flux present along the field line highlighted inblue in Figures 7 and 9 (which is the total flux lost laterally)with the time-integrated total flux inside the pore at the bottomof the domain (which is the total incoming flux). The time inte-gration of the flux is calculated for the first wave front over oneperiod T for all locations E t ( x , z ) = (cid:90) t ( z ) t ( z ) E ( x , z , t ) dt , (14) Article number, page 7 of 12 & A proofs: manuscript no. main
Fig. 9.
Snapshot of the wave energy flux parallel (left) and perpendicular (right) with respect to the magnetic field. The color range is saturated.The solid (dashed) green lines show the first theoretical wave fronts of the fast (slow) waves; the gray lines show the magnetic field lines. The redbars below the x -axis indicate the driver location. The blue lines highlight the field line considered for the analysis in Figure 8. Movies of the fulltime sequence are available online. where E is the wave energy flux according to Equation 11 and t ( z ) and t ( z ) is the time of the beginning and end, respectively,of the first wave front at height z , with t ( z ) = t ( z ) + T .For the calculation of the escaped flux we chose a field linerooted slightly outside the driver region in order to be sure thatall flux at that location has exited the pore. We then integratethe time-integrated flux components along this field line fromthe root of the field line at the bottom of the domain until height z , before calculating the time-integrated total flux. The total es-caped flux is then E t , esc ( z ) = (cid:34)(cid:90) l ( z )0 E t , (cid:107) ( x , z ) dl (cid:35) + (cid:34)(cid:90) l ( z )0 E t , ⊥ ( x , z ) dl (cid:35) / , (15)where l ( z ) is the length of the field line at height z and the in-tegrals of the fluxes are taken along the field line. Here E t , (cid:107) and E t , ⊥ are the parallel and perpendicular components of the time-integrated energy flux (Equation 14) with respect to the magneticfield (and therefore the field line), with E t , (cid:107) ⊥ E t , ⊥ . The integra-tion is done before calculating the absolute value to allow fluxwith opposing signs to cancel out.Similarly, the total flux contained in the pore at the bottomof the domain is calculated by E t , bot = (cid:34)(cid:90) x l E t , x ( x , z = z bot ) dx (cid:35) + (cid:34)(cid:90) x l E t , z ( x , z = z bot ) dx (cid:35) / , (16)where x l = .
22 Mm is the x -position of the field line root, z bot = − .
095 Mm is the z -location of the bottom of the domain, andthe integrals are taken horizontally across the pore at the bottomof the domain. Here E t , x and E t , z are the x - and z -components ofthe time-integrated energy flux, with E t , x ⊥ E t , z .The e ff ect of wave leakage on the damping is then estimatedby e damp ( z ) = − E t , esc ( z ) E t , bot . (17)The result of Equation 17 is shown in Figure 8 (right) as theblue dashed line. There is a significant di ff erence between this line and the line showing missing flux when only consideringgeometric spreading (orange dash-dotted line). Both methods areestimates, and we expect the actual e ff ect of lateral wave leakageto lie between these lines.
4. Conclusions and discussion
We created a MHS model close to equilibrium, which was in-spired by observational data of a solar pore (GM21) and investi-gated possible damping mechanisms by driving the model witha vertical velocity perturbation at the bottom of the domain. Wefound that, even if viscosity, resistivity, or thermal conductionare included, the strong damping from the observations couldnot be reproduced at all by using a driver that covers the wholebottom boundary. When switching to a localized driver, however,the results show strong damping in our simulations. This damp-ing occurs because of a) geometric spreading, where the flux isspread over a wider area due to diverging field lines and b) lat-eral wave leakage, where waves leave the pore. Therefore, evenif only considering classic wave e ff ects, significant damping canbe achieved. Wave leakage at the edge of a solar pore was indeedalready observed by Stangalini et al. (2011). It was mentioned in Section 2.2 there are di ff erences betweenour model and the observational data example pore (GM21, pore3). The di ff erences in density and pressure profiles mainly leadto di ff erences in characteristic wave speeds. This does not a ff ectthe damping due to geometric spreading, as this damping mech-anism is only dependent on the magnetic field structure, which issimilar to the observations, with nearly vertical inclination insidethe pore and nearly horizontal field lines outside.An important point we have to note, however, is the soundspeed profile, as shown in Figure 10. In our model, the soundspeed generally increases with height, whereas it is the oppo-site for the observations. In addition, there is a strong horizontalstructuring, with lower speeds at the center and the border of thepore. From applying Snell’s refraction law, as also discussed in Article number, page 8 of 12.M. Riedl et al.: Finding the mechanism of wave energy flux damping in solar pores using numerical simulations
Fig. 10.
Sound speed of the initial atmosphere. Gray lines show mag-netic field lines. The red bar below the x -axis indicates the driver loca-tion for simulations with localized driver. The blue line highlights thefield line considered for the analysis in Figure 8. Contours for the soundspeed are shown in thick black lines. the context of sunspots by Khomenko & Collados (2006), weknow that waves travelling into a medium with higher phasespeed refract away from the line perpendicular to the constant-phase-speed-line. If in our simulations the fast (acoustic) wavesare propagating along the diverging field lines, they are refractedaway from the pore. Therefore, should the fast lateral wavesin our simulations exclusively occur because of refraction, wewould not expect acoustic waves escaping laterally for the ob-servations of pores like in GM21. The e ff ect of lateral leakingfor magnetic waves should be the same, however, as the Alfvénspeed profile in our simulations is similar to the observations.Evidence of at least some fast wave refraction occurring inour simulations is seen in the amplitude of the wave energy flux,where the amplitude is increased at the center and the bound-ary of the pore compared to the region in between. Those re-gions coincide with regions of lower sound speeds and wavesare therefore refracted toward those regions. The increased am-plitude in the pore boundary therefore does not occur becauseof sausage surface waves. Since the sound speed is higher at thelocation just outside the pore, waves that are located outside thepore would be refracted into the pore. This could be one of thereasons why the energy flux profile increases with height for thefull driver (Figure 6), as there are ample waves present outsidethe pore to be refracted. In addition, fast wave energy flux thatescaped from the pore was eventually refracted down toward thebottom of the domain in the simulations with localized driver.This can be seen in Figure 9 (right), where the perpendicular fluxcomponent for the fast waves outside the pore is mainly positiveand therefore directed downward, considering the nearly hori-zontal field lines. This refraction of fast waves is similar to whatwas found by Khomenko & Collados (2006).The observations of pore 3 GM21 also show higher energyflux concentrations at the pore boundaries. Contrary to the eventsin our simulations, it was found that these flux concentrationsare due to surface sausage modes. This could possibly promoteadditional lateral wave leakage as flux already present at the edgeof the pore could more easily escape.A crucial di ff erence between observations and simulationsis that due to the cadence of the instruments, GM21 were onlyable to investigate slow waves, whereas in this paper we have a combination of slow and fast waves. By splitting the energy fluxinto magnetic (Poynting) and hydrostatic contributions, slow andfast waves could have been studied separately. However, mostof the slow waves in our simulations with localized driver stemfrom the sharp edge of the step-function, causing most of theslow waves being concentrated just inside and atop the consid-ered field line marking the boundary of the pore in our analysis(blue highlighted field line in e.g. Figure 9), with little slow waveflux inside the rest of the pore. We therefore only considered thetotal flux for our analysis, as our estimate for the influence ofdamping due to lateral wave leakage (Equation 17) would nothave worked for slow waves alone. On the other hand, there wasno need to exclude slow waves from the same analysis as themagnitude of the Poynting flux is about three orders of magni-tude smaller than the hydrostatic component.In our model β > β < ff erent shapethan for fast waves. While, according to our results, slow wavesalso leave the pore, it is possible that due to this di ff erent shapefewer low β slow acoustic waves (observations) would leave thepore than fast acoustic waves in our simulations. However, theslow acoustic waves in observations are still comparable to thefast acoustic waves simulated here. It is therefore reasonable toassume that within the pore the slow wave energy flux wouldbe dominant over the fast wave energy flux if our model atmo-sphere had β < ffi culties in observing fastmodes: the fast wave flux would be overshadowed by the slowwave flux. In addition, having a low plasma- β inside the poreinevitably leads to a β = v s = v A ) layer at the border ofthe pore with high β outside. In these layers waves are stronglysubjected to mode conversion (Cally 2005, 2006; Schunker &Cally 2006; Hansen et al. 2015). Whether these mode conver-sions increase the amount of energy flux escaping from the poreor have a channeling e ff ect in the pore will have to be determinedin future work. In this work, we did not account for any radiative losses. Ac-cording to Carlsson & Stein (2002), acoustic waves in the pho-tosphere are much more damped at higher frequencies, meaningthat the impact of this damping mechanism in our simulationswould be larger than for the observations of GM21, who observelonger periods.Our simulations were done on a 2D Cartesian grid. In 2D, the“area” inside the pore at each height is just a 1D line. Therefore,we estimated the damping due to geometric spreading to be pro-portional to 1 / R ( z ) with R ( z ) the distance between the pore axisand a field line. In 3D, however, we expect the wave energy fluxdue to this e ff ect to decrease with 1 / R ( z ) . Estimating the changein e ff ect from 2D to 3D for wave leakage is more di ffi cult. Weassume that it is dependent on the ratio of the area inside the poreto the area that has been available for flux to escape, which is themantle of the pore up to a specific height. This ratio is R ( z ) / l ( z ) in2D and R ( z ) π/ (2 R ( z ) π l ( z )) in 3D, with l ( z ) describing the lengthof the considered field line from the root up to a certain height z .Therefore, the dependence R ( z ) / l ( z ) can also be assumed for 3D. Article number, page 9 of 12 & A proofs: manuscript no. main
The increase in e ffi ciency of geometric spreading for 3D couldaccount for the di ff erence between the damping in our simula-tions and the observed damping. We note that by assuming a 2Dgeometry in our simulations we have excluded the possibility ofAlfvén waves. As discussed above, there are both slow and fast waves presentin our simulations. The slow waves are predominantly excitedat the edge of the step-function driver. Simulations using aGaussian-shaped driver instead show that slow waves are excitedat the flank of the Gaussian, mostly at the steepest location. Thisleads to the conclusion that any kind of localized vertical driverwould excite both slow and fast waves. Therefore, we also expectboth kinds of waves to be present in the photosphere at all times.While slow modes have been observed in the photosphere manytimes, temporal resolution has so far limited similar studies forfast waves. However, future instruments on the Daniel K. InouyeSolar Telescope (DKIST), European Solar Telescope (EST), andNational Large Solar Telescope (NLST) might provide the ca-dence needed to observe fast waves propagating at an inclinedangle with respect to the magnetic field.Observing the leaking waves as seen in our simulationsmight be challenging as the magnitude of the vertical (line-of-sight) velocity perturbations is roughly a factor of ten lower thanthe perturbations inside the pore. However, since the wave frontsof the leaking waves are inclined from the vertical (as seen inFigure 9), an observer from above would see the integrated ef-fects of waves in di ff erent phases (i.e., positive and negative ve-locities within the same pixel). This would lead to spectral linebroadening. The possibility to observe the leaking waves us-ing this e ff ect can be investigated using forward modeling tech-niques, such as the FoMo code developed by Van Doorsselaereet al. (2016). Acknowledgements.
The authors thank the referee for their constructive com-ments. JMR and TVD have received funding from the European Research Coun-cil (ERC) under the European Union’s Horizon 2020 research and innovationprogramme (grant agreement No. 724326). CAG-M, DBJ, and SDTG are grate-ful to Invest NI and Randox Laboratories Ltd. for the award of a Research & De-velopment Grant (059RDEN-1) that allowed the research framework employedto be developed. DBJ and SDTG also wish to acknowledge the UK Science andTechnology Facilities Council (STFC) for funding under the Consolidated GrantST / T00021X / References
Balthasar, H., Collados, M., & Muglach, K. 2000, Astronomische Nachrichten,321, 121Bogdan, T. J., Carlsson, M., Hansteen, V. H., et al. 2003, ApJ, 599, 626Borrero, J. M., Pastor Yabar, A., Rempel, M., & Ruiz Cobo, B. 2019, A&A, 632,A111Botha, G. J. J., Arber, T. D., Nakariakov, V. M., & Zhugzhda, Y. D. 2011, ApJ,728, 84Cally, P. S. 1986, Sol. Phys., 103, 277Cally, P. S. 2001, ApJ, 548, 473Cally, P. S. 2005, MNRAS, 358, 353Cally, P. S. 2006, Philosophical Transactions of the Royal Society of LondonSeries A, 364, 333Cally, P. S. & Khomenko, E. 2019, ApJ, 885, 58Cameron, R., Schüssler, M., Vögler, A., & Zakharov, V. 2007, A&A, 474, 261Carlsson, M. & Stein, R. F. 1997, ApJ, 481, 500Carlsson, M. & Stein, R. F. 2002, in ESA Special Publication, Vol. 505, SOL-MAG 2002. Proceedings of the Magnetic Coupling of the Solar AtmosphereEuroconference, ed. H. Sawaya-Lacoste, 293–300Centeno, R., Collados, M., & Trujillo Bueno, J. 2006, ApJ, 640, 1153Chen, S.-X., Li, B., Shi, M., & Yu, H. 2018, ApJ, 868, 5Cho, K. S., Bong, S. C., Chae, J., et al. 2013, Sol. Phys., 288, 23 Dedner, A., Kemm, F., Kröner, D., et al. 2002, Journal of Computational Physics,175, 645Dorotoviˇc, I., Erdélyi, R., Freij, N., Karlovský, V., & Márquez, I. 2014, A&A,563, A12Felipe, T., Kuckein, C., González Manrique, S. J., Milic, I., & Sangeetha, C. R.2020, ApJ, 900, L29Felipe, T., Kuckein, C., & Thaler, I. 2018, A&A, 617, A39Fujimura, D. & Tsuneta, S. 2009, ApJ, 702, 1443Garcia de La Rosa, J. I. 1987, Sol. Phys., 112, 49Geeraerts, M., Van Doorsselaere, T., Chen, S.-X., & Li, B. 2020, ApJ, 897, 120Gilchrist-Millar, C. A., Jess, D. B., Grant, S. D. T., et al. 2021, PhilosophicalTransactions of the Royal Society of London Series A, 379, 20200172Goedbloed, J., Keppens, R., & Poedts, S. 2019, Magnetohydrodynamics of Lab-oratory and Astrophysical Plasmas (Cambridge University Press; Cambridge)Goedbloed, J. P. H. & Poedts, S. 2004, Principles of Magnetohydrodynamics:With Applications to Laboratory and Astrophysical Plasmas (Cambridge Uni-versity Press)Grant, S. D. T., Jess, D. B., Moreels, M. G., et al. 2015, ApJ, 806, 132Grant, S. D. T., Jess, D. B., Zaqarashvili, T. V., et al. 2018, Nature Physics, 14,480Hansen, S. C., Cally, P. S., & Donea, A.-C. 2015, Monthly Notices of the RoyalAstronomical Society, 456, 1826Jaeggli, S. A., Lin, H., Mickey, D. L., et al. 2010, Mem. Soc. Astron. Italiana,81, 763Jess, D. B., Morton, R. J., Verth, G., et al. 2015, Space Sci. Rev., 190, 103Jess, D. B., Reznikova, V. E., Van Doorsselaere, T., Keys, P. H., & Mackay, D. H.2013, ApJ, 779, 168Jess, D. B., Snow, B., Houston, S. J., et al. 2020, Nature Astronomy, 4, 220Kato, Y., Steiner, O., Hansteen, V., et al. 2016, ApJ, 827, 7Keys, P. H., Morton, R. J., Jess, D. B., et al. 2018, ApJ, 857, 28Khomenko, E. & Cally, P. S. 2019, ApJ, 883, 179Khomenko, E. & Collados, M. 2006, ApJ, 653, 739Krishna Prasad, S., Jess, D. B., Van Doorsselaere, T., et al. 2017, ApJ, 847, 5Lamb, H. 1909, Proceedings of the London Mathematical Society, s2-7, 122Maltby, P., Avrett, E. H., Carlsson, M., et al. 1986, ApJ, 306, 284Mignone, A., Bodo, G., Massaglia, S., et al. 2007, ApJS, 170, 228Mignone, A., Tzeferacos, P., & Bodo, G. 2010, Journal of ComputationalPhysics, 229, 5896Moreels, M. G., Goossens, M., & Van Doorsselaere, T. 2013, A&A, 555, A75Morton, R. J., Erdélyi, R., Jess, D. B., & Mathioudakis, M. 2011, ApJ, 729, L18Morton, R. J., Verth, G., Jess, D. B., et al. 2012, Nature Communications, 3,1315Nakariakov, V. M. & Zimovets, I. V. 2011, ApJ, 730, L27Ossendrijver, M. 2003, A&A Rev., 11, 287Riedl, J. M., Van Doorsselaere, T., & Santamaria, I. C. 2019, A&A, 625, A144Roberts, B. 2006, Philosophical Transactions of the Royal Society of LondonSeries A, 364, 447Ruiz Cobo, B. & del Toro Iniesta, J. C. 1992, ApJ, 398, 375Schmitz, F. & Fleck, B. 1998, A&A, 337, 487Schunker, H. & Cally, P. S. 2006, MNRAS, 372, 551Sobotka, M. 2003, Astronomische Nachrichten, 324, 369Stangalini, M., Del Moro, D., Berrilli, F., & Je ff eries, S. M. 2011, A&A, 534,A65Thomas, J. H. & Weiss, N. O. 2004, ARA&A, 42, 517Van Doorsselaere, T., Antolin, P., Yuan, D., Reznikova, V., & Magyar, N. 2016,Frontiers in Astronomy and Space Sciences, 3, 4Yu, D. J., Van Doorsselaere, T., & Goossens, M. 2017, A&A, 602, A108 Article number, page 10 of 12.M. Riedl et al.: Finding the mechanism of wave energy flux damping in solar pores using numerical simulations
Appendix A: Simulations with driver below thecutoff frequency
In order to focus on propagating waves all simulations in thispaper have so far been conducted with a driver frequency wellabove the expected cuto ff frequency of the model atmosphere,and therefore also with a period much smaller than the five-minute waves observed by GM21. This means that for the cur-rent study, all e ff ects of the cuto ff frequency have been ignored.However, as mentioned in Section 1, even if the waves of GM21definitely have a propagating character, they might be altered topartly evanescent waves by the existence of the cuto ff frequency(Centeno et al. 2006). In this section we explore the possibilityof damping due to evanescent waves by conducting the same twoexperiments as before, namely simulations with full driver andlocalized driver, but with a lower driver frequency. Appendix A.1: Cutoff frequency and new driver period
It is commonly accepted that acoustic waves with frequenciesbelow the cuto ff frequency are not allowed to propagate, but arestanding and evanescent. However, it is di ffi cult to define an ex-act value for the cuto ff frequency, and numerous di ff erent def-initions exist. Centeno et al. (2006) show that when radiativelosses are involved there is no clear cuto ff frequency that dis-tinguishes between fully propagating or fully evanescent waves.Felipe et al. (2018) compared analytical definitions for the cut-o ff frequency suitable for sunspot umbrae from Lamb (1909),Schmitz & Fleck (1998), and Roberts (2006) to the observedcuto ff . The results generally agree. Using the same analyticalexpressions as discussed in Felipe et al. (2018) on the observa-tional data obtained by GM21 for pore 3 shows that waves witha period of five minutes indeed have a lower frequency than thecuto ff frequency for at least most of the observed domain.According to the same equations, a driver period of five min-utes would still result in a frequency above the cuto ff frequencyfor our model atmosphere. To mimic the conditions of the obser-vations, we choose a longer driver period of T = Appendix A.2: Results
Figure A.1 shows the height–time graph of the wave energy fluxparallel to the magnetic field at the axis of the pore for the sim-ulation with localized driver. The characteristic speeds (startingfrom steepest: fast speed v fa ( θ = π/ = ( v A + v s ) / , sound speed,Alfvén speed, cusp speed) are plotted as black lines, while thecontour at value zero is shown in red. The initial part of the firstwave (i.e., the initial disturbance where the flux is above zero forthe first time) propagates with the sound speed (black dashed lineoverplotted on first red line) as it did for the propagating waves inSection 3. Then, however, the waves get altered by the e ff ects ofthe cuto ff frequency to approximately standing waves within lessthan half a driver period, as can be seen from the nearly verticalfeatures in the figure. This is not what was observed in GM21,who found clear evidence of propagating waves. The di ff erencemight be accounted for by the neglect of radiative losses in oursimulations.We performed the same study for the wave energy flux damp-ing as in Sections 3.1 and 3.2, but for the low-frequency driver.Figure A.2 shows the results for the full driver, while Figure A.3shows the results for the localized driver. It is immediately appar-ent that the energy flux for the full driver is now heavily damped Fig. A.1.
Parallel wave energy flux as a function of height and timeat the pore axis for the simulation with localized driver with a periodof 7 minutes. The black lines show (from steepest to flattest) the fastspeed (dash-dotted line), sound speed (dashed line), Alfvén speed (dot-ted line), and cusp speed (solid line). The red lines show the contoursfor zero flux. The frequency of the energy flux is approximately dou-bled compared to the driver period because of phase di ff erence between p (cid:48) and v (cid:48) (see Equation 11). Fig. A.2.
Same as Figure 6, but for a driver period of 7 minutes. as well, about the same amount as the energy flux for the local-ized high-frequency driver (Figure 8). The energy flux for the lo-calized low-frequency driver (Figure A.3) is damped even more,probably because the damping with height is not decreased byinward refracted waves as for the full driver.
Appendix A.3: Discussion
It is obvious that the choice of driver frequency strongly a ff ectsthe damping in our simulations. However, whether this is purelydue to evanescent waves is not fully clear.On the one hand, the dash-dotted orange curve in Figure A.3,which shows the damping without e ff ect of geometric spreading,strongly follows the solid green line, which is the full dampingin our simulation with the localized low-frequency driver. Thishints that geometric spreading has little to no e ff ect in this case.At the same time the dashed blue line, which is an estimate forthe influence of lateral leakage, is nearly constant, meaning thatthis e ff ect is also not very strong. Therefore, a crucial dampingmechanism is missing, which is likely the reflection of wavesdue to the cuto ff frequency. Article number, page 11 of 12 & A proofs: manuscript no. main
Fig. A.3.
Same as Figure 8, but for a driver period of 7 minutes.
On the other hand, these new simulations and their analysisare subject to some limiting factors. First of all, due to the lowfrequency, the wavelengths of the resulting waves are signifi-cantly longer than the size of the computational domain. Thiscould lead to strange boundary e ff ects influencing the results.Since the ratio of the wavelength to the size of the pore (which issmaller in our model than in the observations) also changes dras-tically, this could account for the decreased e ff ects of geometricdamping and lateral wave leakage. In addition, due to the wavesstarting at some final time t , there are no waves present in thedomain before the first waves reach a certain height (i.e., left ofthe first red line in Figure A.1). Therefore, when integrating thewave energy flux over time, the lower integration boundary t ( z )was chosen by using a relative threshold to determine the onsetof the first wave at every height. This line basically coincideswith the sound speed line (dashed) in Figure A.1. The upper in-tegration boundary was then determined by t ( z ) = t ( z ) + T , with T being the driver period. E ff ectively, the time integration for thesimulations with high-frequency driver was done over the firstperiod of the wave, as a translation of t ( z ) by T =
30 s resultedin a t ( z ) being located right in front of the next wave train. Thisis not the case for the low-frequency waves because they changefrom propagating to standing waves within the first wave period,meaning that their steepness changes in Figure A.1. Therefore,it is not clear over which time period the integration should beperformed, and the choice might a ff ect the shape of the dampingcurves in Figures A.2 and A.3.Moreover, even if the limitations listed above have little to noe ff ect, there are still no propagating waves in our low-frequencysimulations, as opposed to the observations of GM21. Therefore,the damping in the low-frequency simulations due to evanescentwaves is expected to be much stronger than for the observations,where the waves were at least partly propagating. This validatesthe study of the other damping mechanisms presented in this pa-per.ect, there are still no propagating waves in our low-frequencysimulations, as opposed to the observations of GM21. Therefore,the damping in the low-frequency simulations due to evanescentwaves is expected to be much stronger than for the observations,where the waves were at least partly propagating. This validatesthe study of the other damping mechanisms presented in this pa-per.