© ESO 2019

1. Introduction

Solar magnetically active regions exhibit many oscillatory phenomena. Their study is a relevant topic in astrophysics, since they can be used to infer the properties of the media where they propagate using seismology techniques (De Moortel & Nakariakov 2012), and they have also been proposed as one of the candidates to explain chromospheric and coronal heating (Alfvén 1947). Since the first detection of waves in sunspots (Beckers & Tallant 1969), many oscillatory phenomena have been reported in various magnetic structures, from the photosphere to the outer layers. In the case of sunspots, photospheric umbral oscillations show a period of five minutes (Bhatnagar et al. 1972), while at the umbral chromosphere their period is reduced to three minutes (Beckers & Schultz 1972) and they can form shocks (Lites 1986; Centeno et al. 2006). All these manifestations are associated to a global process of wave propagation in sunspots (see Khomenko & Collados 2015, for a review). In fact, several authors have provided evidence of wave propagation from the photosphere to the corona (Reznikova et al. 2012; Krishna Prasad et al. 2015; Zhao et al. 2016).

One of the basic observational findings in umbral atmospheres is the reduction of the period of the oscillations with increasing height, from five minutes in the photosphere to three minutes in the chromosphere (e.g., Lites 1992). Three mechanisms have been proposed to explain the origin of the three-minute oscillations: propagation of waves with frequency above the cutoff directly from the photosphere, a chromospheric resonant cavity, and their formation in the wake of propagating shocks. The latter is found in nonlinear one-dimensional simulations (Hansteen 1997; Bard & Carlsson 2010), although Ulmschneider et al. (2005) concluded that the shock mergers leading to this effect are unrealistic. In this work we focus on the two other models. Centeno et al. (2006) fitted the phase spectrum between the Doppler velocity from photospheric and chromospheric lines with a model of linear vertical propagation of slow magnetoacoustic waves in a stratified atmosphere with radiative losses. They conclude that waves with frequency above the cutoff value can propagate upwards, and due to the rapid drop of the density, their amplification is stronger than that of evanescent low-frequency waves (with frequency below the cutoff). This physical effect was previously modeled by Fleck & Schmitz (1991). In their numerical simulations, waves were driven at the bottom boundary by a long-period sinusoidal function in the vertical velocity. This driver also generates high-frequency waves, which propagate upward and dominate the power spectra at the chromosphere. Later, Felipe et al. (2010a) modeled this mechanism by explicitly using a driver that introduced waves in the three-minute band below the photosphere, and Felipe et al. (2011) performed observationally driven simulations that reproduce the observed chromospheric three-minute oscillations by means of high-frequency waves propagating from deeper layers.

The second mechanism explored in this work is that of a chromospheric acoustic resonator (Zhugzhda & Locans 1981; Zhugzhda et al. 1983; Gurman & Leibacher 1984). Slow magnetoacoustic waves are trapped between the height where the cutoff frequency is below the wave frequency (around the temperature minimum) and the steep temperature gradients at the transition region. These boundaries are permeable and form a resonant cavity. Zhugzhda (2008) analyzed the effect of this resonator taking into account the changes in the temperature with height. This mechanism has also been studied through numerical simulations using the ideal magnetohydrodynamic (MHD) approximation (Botha et al. 2011; Snow et al. 2015) and using a two-fluids model (Wójcik et al. 2018). These latter studies support the existence of an acoustic resonator that generates the observed chromospheric oscillations in the three-minute band. It must be noted that the resonator does not prevent the effect of the cutoff frequency.

Observational works have reported the presence of three-minute oscillations in the umbral photosphere (e.g., Lites & Thomas 1985; Abdelatif et al. 1986; Centeno et al. 2006; Felipe et al. 2010b). Chae et al. (2017) suggested that the three-minute oscillations are generated by magnetoconvection occurring in the umbra (in light bridges and umbral dots), and reported photospheric three-minute oscillations with energy flux large enough to explain the chromospheric oscillations as upwardly propagating slow magnetoacoustic waves. Stangalini et al. (2012) found a three-minute wave enhancement in the photosphere of a pore, whose amplitude increases with the magnetic field strength.

In this paper, we numerically investigate wave propagation in sunspot umbrae. We performed a set of numerical experiments that explore two of the mechanisms discussed above: the chromospheric acoustic resonator and the propagation of waves in the three-minute band directly from the photosphere. We evaluate the effect of various physical ingredients on the results. In Sect. 2 we describe the numerical simulations and the results are presented in Sect. 3. Sections 4 and 5 contain the discussion and conclusions, respectively.

2. Numerical simulations

Wave propagation from the solar interior to the corona has been computed using the code MANCHA (Khomenko & Collados 2006; Felipe et al. 2010a). In this study, we restricted the use of the code to solve the ideal nonlinear MHD equations in a one-dimensional atmosphere, but keeping the three-dimensional components of the vectors. The equilibrium state is explicitly removed from the equations. The code solves the following system of equations for perturbations:

$$\begin{array}{c}\hfill \frac{\partial {\rho }_1}{\partial t}+\mathrm{\nabla }[({\rho }_0+{\rho }_1){\mathit{v}}_{\mathbf{1}}]=0,\end{array}$$(1)

$$\begin{array}{cc}& \frac{\partial [({\rho }_0+{\rho }_1){\mathit{v}}_1]}{\partial t}+\mathrm{\nabla }·[({\rho }_0+{\rho }_1){\mathit{v}}_1{\mathit{v}}_1\hfill \\ \hfill & +(p_1+\frac{{\mathit{B}}_1^2}{2{\mu }_0}+\frac{{\mathit{B}}_1{\mathit{B}}_0}{{\mu }_0})\mathit{I}-\frac{1}{{\mu }_0}({\mathit{B}}_0{\mathit{B}}_1-{\mathit{B}}_1{\mathit{B}}_0-{\mathit{B}}_1{\mathit{B}}_1)]={\rho }_1\mathit{g},\hfill \end{array}$$(2)

$$\begin{array}{cc}& \frac{\partial e_1}{\partial t}+\mathrm{\nabla }·[{\mathit{v}}_{\mathbf{1}}((e_0+e_1)+(p_0+p_1)\hfill \\ \hfill & +\frac{|{\mathit{B}}_0+{\mathit{B}}_1{|}^2}{2{\mu }_0})-\frac{1}{{\mu }_0}({\mathit{B}}_0+{\mathit{B}}_1)({\mathit{v}}_1·({\mathit{B}}_0+{\mathit{B}}_1))]\hfill \\ \hfill & =({\rho }_0+{\rho }_1)(\mathit{g}·{\mathit{v}}_1)+Q_{\mathrm{rad}},\hfill \end{array}$$(3)

$$\begin{array}{c}\hfill \frac{\partial {\mathit{B}}_1}{\partial t}=\mathrm{\nabla }\times [{\mathit{v}}_1\times ({\mathit{B}}_0+{\mathit{B}}_1)],\end{array}$$(4)

where ρ is the density, p is the gas pressure, v is the velocity, B is the magnetic field, I is the identity tensor, and e is the total energy per unit volume, given by

$$\begin{array}{c}\hfill e=\frac{1}{2}\rho v^2+\frac{p}{\gamma -1}+\frac{B^2}{2{\mu }_0}·\end{array}$$(5)

Subindex 0 indicates the background values, corresponding to an atmosphere in equilibrium, while variables with subindex 1 are the perturbed quantities. The gravitational acceleration and the magnetic permeability are given by g and μ₀, respectively. The term Q_rad represents the radiative losses (see Sect. 3.5). Artificial diffusion terms were added following Vögler et al. (2005).

Umbral models were constructed by combining the convectively stabilized CSM_B interior model (Schunker et al. 2011) with the photospheric and chromospheric atmosphere obtained from the Avrett (1981) umbral model. A one million Kelvin isothermal corona was added following Santamaria et al. (2015), with a sharp temperature increase from the chromospheric value to the coronal temperature in a transition region of 200 km in thickness. Various umbral models were constructed. They differ in the height of the base of the isothermal corona, which ranges from z = 1.8 Mm to z = 2.6 Mm. The photosphere is located at the height z = 0, where the optical depth at 500 nm is unity. Figure 1 shows the temperature, cutoff frequency, and density profiles of those models. The sunspot model M from Maltby et al. (1986) is also plotted as a reference. In the top panel, the color lines indicate the acoustic cutoff of various umbral models (with differences in the height of the transition region) obtained from the expression:

$$\begin{array}{c}\hfill {\omega }_{\mathrm{c}}=\frac{c_{\mathrm{s}}}{2H}(1-2\frac{\mathrm{d}H}{\mathrm{d}z}{)}^{1/2},\end{array}$$(6)

Fig. 1. Temperature (top panel), acoustic cutoff frequency (middle panel), and density (bottom panel) of the sunspot models. They differ in the height of the base of the corona, which is located at 1.8 Mm (orange), 2.0 Mm (green), 2.2 Mm (light blue), 2.4 Mm (dark blue, and 2.6 Mm (violet). The black dashed line corresponds to the magnetoacoustic cutoff of the umbral model with the base of the corona at 2.0 Mm and a vertical magnetic field of 2500 G. The red line is the model without corona. The dotted black line corresponds to the sunspot model M from Maltby et al. (1986).

where $H=c_{\text{s}}^2/(\gamma g)$ is the pressure scale height, γ is the adiabatic index, c_s is the sound speed, and g is the gravity. For the model with the base of the corona at z = 2.0 Mm, we also show the cutoff frequency for magnetoacoustic waves in an atmosphere permeated by a uniform vertical magnetic field of B = 2500 G, as derived by Roberts (2006):

$$\begin{array}{c}\hfill {\omega }_{\mathrm{c}}^B=c_{\mathrm{t}}[\frac{1}{4H^2}(\frac{c_{\mathrm{t}}}{c_{\mathrm{s}}}{)}^4-\frac{1}{2}\gamma g\frac{\partial }{\partial z}(\frac{c_{\mathrm{t}}^2}{c_{\mathrm{s}}^4})+\frac{1}{v_{\mathrm{A}}^2}(N^2+\frac{g}{H}\frac{c_{\mathrm{t}}^2}{c_{\mathrm{s}}^2}){]}^{1/2},\end{array}$$(7)

where

$$\begin{array}{c}\hfill N^2=-g(\frac{1}{\rho }\frac{\partial \rho }{\partial z}-\frac{1}{\gamma p}\frac{\partial p}{\partial z})\end{array}$$(8)

is the squared Brunt–Väisälä frequency, v_A = B/(4πρ)^1/2 is the Alfvén speed, and $c_{\text{t}}=c_{\text{s}}{v}_{\text{A}}/{(c_{\text{s}}^2+v_{\text{A}}^2)}^{1/2}$ is the cusp speed.

We performed several numerical experiments. In some of them, we aim to explore the effect of a chromospheric cavity formed between the temperature minimum and the transition region. In those cases, the top of the computational domain is placed at z = 5.4 Mm. At the top boundary we imposed a perfect matched layer (PML, Berenger 1994), which damps the waves and avoids the presence of spurious reflections. In the PML domain, the MHD equations (Eqs. (1)–(4)) are split into a set of three coupled one-dimensional equations, and a damping coefficient is added to each equation only near the boundary of the corresponding direction. Numerical tests with differences in the height and strength of the PML top boundary were performed, and no significant differences were found in the results. For further details about the PML, we refer the reader to Berenger (1994), Hu (1996), and Felipe et al. (2010a).

In order to study wave propagation in a model with no chromospheric cavity, we used an umbral atmosphere with constant chromospheric temperature above z = 1.4 Mm and with the top boundary located at z = 2 Mm (red line in Fig. 1). In all simulations, the bottom boundary is placed at z = −5.2 Mm. The code employs a variable spatial resolution. The grid is stretched in the vertical coordinate according to the local sound speed. A coarser resolution is used in regions with high sound speed and a finer resolution in regions with low sound speed. The travel time of acoustic waves between adjacent cells is the same throughout the domain, except at layers above the temperature minimum (including chromosphere, transition region, and corona), where we have forced the finest resolution (Δz = 12.4 km). The lowest resolution (Δz = 49.5 km) is set at the bottom boundary.

Waves are driven by specifying a perturbation as a function of time in a few grid points at the bottom boundary. We employed two different drivers. First, we implemented a monochromatic acoustic-gravity wave of a selected period. We imposed perturbations in pressure, density, and velocity obtained from analytical calculations (Mihalas & Mihalas 1984). These perturbations do not take into account the temperature gradient and magnetic field (whose effect is negligible at this depth). According to Khomenko & Cally (2012), the analytical expressions for the perturbations in velocity, pressure, and density are given by

$$\begin{array}{c}\hfill v_{z1}=Vexp(\frac{z}{2H}+k_{\mathit{zi}}z)sin(\omega t-k_{\mathit{zr}}z-k_{x}x),\end{array}$$(9)

$$\begin{array}{c}\hfill \frac{p_1}{p_0}=V|P|exp(\frac{z}{2H}+k_{\mathit{zi}}z)sin(\omega t-k_{\mathit{zr}}z-k_{x}x+{\varphi }_P),\end{array}$$(10)

$$\begin{array}{c}\hfill \frac{{\rho }_1}{{\rho }_0}=V|R|exp(\frac{z}{2H}+k_{\mathit{zi}}z)sin(\omega t-k_{\mathit{zr}}z-k_{x}x+{\varphi }_{\mathrm{R}}),\end{array}$$(11)

$$\begin{array}{c}\hfill v_{x1}=V|U|exp(\frac{z}{2H}+k_{\mathit{zi}}z)sin(\omega t-k_{\mathit{zr}}z-k_{x}x+{\varphi }_U),\end{array}$$(12)

where the amplitude of the wave is controlled by the parameter V. The relative phase shifts and amplitudes follow

$$\begin{array}{c}\hfill |P|=\frac{\gamma \omega }{{\omega }^2-c_{\mathrm{s}}^2{k}_x^2}\sqrt{k_{\mathit{zr}}^2+(k_{\mathit{zi}}+\frac{1}{2H}\frac{(\gamma -2)}{\gamma }{)}^2},\end{array}$$(13)

$$\begin{array}{c}\hfill |R|=\frac{\omega }{{\omega }^2-c_{\mathrm{s}}^2{k}_x^2}\sqrt{k_{\mathit{zr}}^2+(k_{\mathit{zi}}+\frac{(\gamma -1)}{\gamma H}\frac{c_{\mathrm{s}}^2{k}_x^2}{{\omega }^2}-\frac{1}{2H}{)}^2},\end{array}$$(14)

$$\begin{array}{c}\hfill |U|=\frac{c_{\mathrm{s}}^2{k}_x}{\gamma \omega }|P|,\end{array}$$(15)

$$\begin{array}{c}\hfill {\varphi }_P={\varphi }_U=arctan(\frac{k_{\mathit{zi}}}{k_{\mathit{zr}}}+\frac{1}{2H{k}_{\mathit{zr}}}\frac{(\gamma -2)}{\gamma }),\end{array}$$(16)

$$\begin{array}{c}\hfill {\varphi }_{\mathrm{R}}=arctan(\frac{k_{\mathit{zi}}}{k_{\mathit{zr}}}+\frac{(\gamma -1)}{\gamma H{k}_{\mathit{zr}}}\frac{c_{\mathrm{s}}^2{k}_x^2}{{\omega }^2}-\frac{1}{2H{k}_{\mathit{zr}}}).\end{array}$$(17)

For a given frequency ω, the vertical wavenumber k_z is obtained from the dispersion relation for acoustic-gravity waves

$$\begin{array}{c}\hfill k_z=k_{\mathit{zr}}+i{k}_{\mathit{zi}}=\sqrt{\frac{({\omega }^2-{\omega }_{\mathrm{c}}^2)}{c_{\mathrm{s}}^2}-\frac{k_x^2({\omega }^2-N^2)}{{\omega }^2}},\end{array}$$(18)

since we have imposed vertical wave propagation and, thus, the horizontal wavenumber k_x has been set to zero. All variables oscillate with a sinusoidal function for the duration of the simulation, the period of which was set to P = 2π/ω = 300 s.

The second kind of driver is a Ricker wavelet introduced in the vertical velocity at the bottom boundary. Its temporal variation is given by

$$\begin{array}{c}\hfill v_z(t)=A(1-2{\tau }_0^2)e^{-{\tau }_0^2},\end{array}$$(19)

where τ₀ = 0.5ωt − π and A is the amplitude. This driver excites a broadband wave spectrum with its central peak at frequency ω, which was chosen to correspond to a period P = 300 s. Drivers with this temporal evolution have been shown to excite a wave spectrum which resembles the solar observations (e.g., Parchevsky et al. 2008; Felipe et al. 2016).

In this study we performed a set of one-dimensional numerical experiments of wave propagation in a sunspot umbra. The setup of the simulations includes differences in the background atmosphere, the location of the top boundary, the presence of magnetic field and radiative losses, and the amplitude and nature of the wave-driving mechanism. These simulations allow us to explore different aspects about the origin of three-minute oscillations in umbral chromospheres. Tables 1–6 show a summary of the simulations discussed in each section. Some of the details about the implementation of the numerical simulations are addressed in the corresponding section.

3. Results

3.1. Monochromatic driver

In this section we present the results from the simulations shown in Table 1. We introduced a sinusoidal acoustic wave to these numerical simulations with a period of 300 s in the solar interior. The computational domain extends to the corona. Between the temperature minimum and the transition region a chromospheric cavity is formed, as shown by the local values of the acoustic cutoff frequency (green line from Fig. 1). At around z = 1 Mm the acoustic cutoff exhibits a local maximum around 6 mHz. This value is in agreement with the maximum cutoff determined observationally for a sunspot umbra by Felipe et al. (2018). At the transition region (in this umbral model around z = 1.8 Mm) the acoustic cutoff shows a prominent peak. Acoustic waves with frequency in the range 4.5 − 6 mHz can propagate between z ∼ 1 Mm and the transition region, but they are trapped there. The amplitude of the wave has been chosen to be arbitrarily small in order to keep the simulation in the linear regime, and the magnetic field strength has been set to zero. The effects of magnetic field and realistic amplitudes is analyzed below.

Figure 2 illustrates the temporal evolution of the vertical velocity at the photosphere (top panel) and the chromosphere (bottom panel). The photospheric velocity shows the regular pattern of a sinusoidal wave with a five-minute period and constant amplitude. This pattern continues for the duration of the simulation (500 min; Fig. 2 only shows the first half of the simulation), except for the initial ∼50 min. The first wavefront driven at the bottom boundary z = −5.2 Mm needs approximately six minutes to reach the photosphere, and during this time there is no wave perturbation. The first photospheric wavefronts exhibit changes in their amplitude. The changes at the beginning of the simulation are more evident at the chromosphere (bottom panel of Fig. 2). The initial chromospheric wavefronts show high amplitude waves with a period of around three minutes. The amplitude of these three-minute waves progressively decreases, until they completely vanish at t ∼ 230 min. After that time, the chromospheric velocity shows a steady sinusoidal wave with a five-minute period. In the stationary stage of the simulation, an amplitude increase is found between the photosphere and chromosphere due to the drop of the density. This amplification is low because five-minute waves are evanescent (their frequency is below the acoustic cutoff).

Fig. 2. Temporal evolution of the vertical velocity at z = 0.3 Mm (top panel) and at z = 1.4 Mm (bottom panel) for Simulation 3.1.1.

The change of the dominant period with time is clearly illustrated in the wavelet analysis (Torrence & Compo 1998) plotted in Fig. 3. The top panel shows the distribution of the power at the chromosphere as a function of time. Initially, strong power is concentrated in the three-minute band (blue line in the middle panel of Fig. 3). Oscillations in the five-minute band are present throughout the simulation (after the first wavefronts reach the chromosphere), and they become dominant when the three-minute oscillations fade away (red line in the middle panel of Fig. 3).

Fig. 3. Top panel: power as a function of period (vertical axis) and time (horizontal axis) obtained from the wavelet analysis of the vertical velocity at z = 1.4 Mm from Simulation 3.1.1. The power is normalized and shown in a logarithmic scale. The results are reliable for the region within the two white lines. The vertical dashed lines indicate the time steps plotted in the middle panel. Middle panel: normalized power at t = 67 min (blue line) and at t = 430 min (red line) obtained from the wavelet analysis. Bottom panel: normalized power for the full temporal series (500 min) obtained from the FFT.

Figure 4 shows the results of the spectral analysis of Simulation 3.1.2. This simulation has the same setup as the previous case, except that this model does not include a transition region or corona. Similarly to the previous simulation, power in the three-minute band appears at the beginning of the computation, but it rapidly vanishes. The stationary state of the simulation is reached after just 30 min. During most of the simulation the chromospheric oscillations are dominated by a sharp five-minute peak. Even though we are introducing a monochromatic sinusoidal oscillation in the bottom boundary, the driven wave is not purely monochromatic during the first time steps. The atmosphere is initially at rest, and the wave is only present at a few grid points. Due to the finite extension of the wave pattern, other frequencies are inevitably excited, including high-frequency waves. These high-frequency waves propagate to higher layers and their amplitude increases due to the drop of the density. Above z = 0, this increment is higher than that of the evanescent five-minute waves (see Fig. 1), and they become dominant at the chromosphere. It is not only waves with frequency above the local maximum in the acoustic cutoff frequency (∼6 mHz) that can reach the chromosphere. Waves with frequency between 4.5 and 6 mHz are evanescent in a small region, but they can leak to the chromosphere where they can freely propagate again.

Fig. 4. Same as Fig. 3 but for Simulation 3.1.2. In this simulation there is no transition region.

In Simulation 3.1.2 (without corona), the waves escape the computational domain through the top boundary, and a stationary simulation is soon obtained as the wave driven at the bottom boundary becomes purely monochromatic. However, in the simulation with transition region and corona (Simulation 3.1.1), high-frequency waves are partially reflected at the transition region due to the strong temperature gradient. These reflected waves travel downward and reach the local maximum in the acoustic cutoff frequency (Fig. 1), where those waves with frequencies below 6 mHz are partially reflected again. This way, a resonance is formed in the chromospheric cavity. The power of waves with a period of 184.3 s is greatly enhanced, and they produce a sharp power peak, as seen in the power obtained from the fast Fourier transform (FFT) of the whole simulation (bottom panel of Fig. 3). Since the temperature gradients are not perfect reflecting surfaces, the three-minute waves progressively leak from the cavity, until their power vanishes.

In summary, the driver employed in these simulations provides a constant source of five-minute waves, whereas some small amount of energy goes to other frequencies during the initial time steps. Those frequencies above the cutoff value can reach the chromosphere and excite the resonant mode. As they leak the cavity, their amplitude decreases (see bottom panel of Fig. 2), and the simulation reaches a stationary state where only five-minute power is found.

3.2. Monochromatic driver with realistic wave amplitudes

Observations show that vertical velocity oscillations in a sunspot umbra have an amplitude of a few hundred meters per second (e.g., Bhatnagar et al. 1972; Soltau et al. 1976; Lites et al. 1982; Centeno et al. 2006; Felipe et al. 2010b). In the previous section we show that monochromatic five-minute waves with low amplitude (in the linear regime) cannot excite chromospheric three-minute oscillations through resonances in a chromospheric cavity. In this section we aim to explore this effect using realistic amplitudes.

The amplitude of the waves in the atmosphere is controlled through the amplitude imposed on the driver at the bottom boundary. We have chosen the appropriate value which produces a wave with 200 m s⁻¹ amplitude at z = 0 Mm (the amplitude of the vertical velocity at the bottom boundary is 0.3 m s⁻¹). As discussed previously, the sinusoidal acoustic wave introduced as a driver initially excites a broad wave spectrum, and the chromosphere is dominated by waves in the three-minute band during the first ∼100 min as a result of the resonant cavity. The use of a high-amplitude driver in combination with the resonance generates highly nonlinear waves with period in the three-minute band. These shocks are produced by the propagation of waves in the three-minute band (above the acoustic cutoff) from deeper layers generated by the driver during the initial stages of the simulation (see previous section). Their amplitude increases according to the drop of the density (this increase is higher than that of nonpropagating evanescent waves) and develops into shocks at the chromosphere. These shocks are not the result of the unrealistic merging of high-frequency waves in one-dimensional simulations (Ulmschneider et al. 2005). Due to these nonlinearities, we have not been able to compute a simulation initialized with high amplitude and long enough to reach the stationary regime. Instead, this simulation was started with a low amplitude (similar to Simulation 3.1.1) and then amplitude was gradually increased. When the appropriate amplitude was reached for the five-minute wave, the driver no longer excited other frequencies, and the five-minute wave was purely monochromatic. This way, we were able to explore the effect of a monochromatic wave with a five-minute period and realistic amplitude on a chromospheric cavity.

The duration of the simulation is 830 min. During the initial 400 min the amplitude of the bottom boundary driver is increasing, until it generates photospheric five-minute oscillations with 200 m s⁻¹ amplitude. The amplitude of the driver then remains constant, and the simulation reaches a stationary regime. Figure 5 shows the power spectra of the stationary part of the simulation at the photosphere (blue line) and the chromosphere (red line). At both heights the power is concentrated at five minutes. These waves are evanescent, and their amplitude increases with height at a rate slower than that of propagating high-frequency waves. No trace of waves in the three-minute band is found. We can conclude that five-minute waves cannot excite oscillations in the three-minute band, even with the presence of a chromospheric cavity and realistic wave amplitudes. The excitation of standing waves in the chromospheric cavity requires the use of a broadband driver (Botha et al. 2011). For simplicity, in the following we focus on waves with low amplitude.

Fig. 5. Power spectra of the vertical velocity from the stationary part of Simulation 3.2.1 at z = 0.3 Mm (blue line) and z = 1.4 Mm (red line). Both spectra are normalized to the maximum value at z = 1.4 Mm.

3.3. Dependence of the resonance on the size of the cavity

In this section we aim to evaluate the frequency of the resonance for various atmospheric models. The steep temperature gradient of the transition region (which produces the reflection at the top of the chromospheric cavity) is shifted to different heights. We have introduced the broadband driver defined in Eq. (19), since in the previous sections we showed that the excitation of the resonant mode requires the presence of that frequency in the driving mechanism.

The top panel of Fig. 6 illustrates the power spectra of the chromospheric velocity for five different models (Simulations 3.3.1–3.3.5). All of them show a sharp peak, but they differ in the period of the peak, the power, and the gradient of the spectra towards higher frequencies. The period of the resonance is approximately given by (e.g., Freij et al. 2016)

$$\begin{array}{c}\hfill P_{\mathrm{r}}\approx \frac{2L}{{\overline{v}}_{c}n},\end{array}$$(20)

Fig. 6. Top panel: power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with waves excited by a broadband driver and with the base of the corona at different heights. The color code is the same as that in Fig. 1, and the base of the corona is located at 1.8 Mm (orange, Simulation 3.3.1), 2.0 Mm (green, Simulation 3.3.2), 2.2 Mm (light blue, Simulation 3.3.3), 2.4 Mm (dark blue, Simulation 3.3.4), and 2.6 Mm (violet, Simulation 3.3.5). Bottom panel: power spectra of the vertical velocity at z = 1.4 Mm from Simulation 3.3.6 using a broadband driver and without transition region and corona. All the power spectra are normalized to the maximum value of the power from Simulation 3.3.2 (height of the bottom of the corona at 2.0 Mm, nonmagnetic, and without radiative losses).

where L in the length of the cavity, ${\overline{v}}_{\mathrm{c}}$ is the average speed of the wave (for acoustic waves, like this simulation, this is the average of the sound speed ${\overline{c}}_{\mathrm{s}}$), and n is the harmonic number. As expected, the longer the cavity the higher the period of the resonance. The measured power peak corresponds to the fundamental mode (n = 1). For each simulation, assuming that the top of the cavity is given by the maximum of the acoustic cutoff (z₁, Fig. 1), we can iteratively determine L = z₁ − z₀, where z₀ is the height of the bottom of the cavity. We start with an initial guess for z₀ and compute ${\overline{c}}_{\mathrm{s}}$ between that height and z₁. This value of the average sound speed and the P_r obtained from the simulation are introduced in Eq. (20) in order to derive a new z₀. This process is repeated until convergence is found. The bottom of the cavity ranges from z = 813 to z = 1291 km above the photosphere, being lower for the atmospheres with the transition region at a lower height. This is in agreement with the acoustic cutoff profiles. Models with the transition region at a lower height exhibit a resonance with higher frequency, and these waves can propagate to lower heights inside the cavity (Fig. 1).

The relative power of the peaks obtained from each simulation depends on the duration of the simulation. Figure 6 shows the power spectra for the whole temporal series of 500 min. During the initial temporal steps the power of the cases with a lower transition region is stronger, but near the end of those simulations their waves have leaked from the chromospheric cavity. For the models with a higher transition region, the power inside the cavity lasts longer. For example, according to a wavelet analysis, the chromospheric power of the model with the highest transition region (Simulation 3.3.5) is reduced to a 10% of its maximum value after 67 min, whereas for the model with the lowest transition region (Simulation 3.3.1) the same power reduction is found after 40 min. We speculate that escaping through the transition region is easier for the higher frequencies of the resonant mode of the atmospheres with a lower transition region. In the case of the atmospheres with a higher transition region, lower-frequency waves are reflected more efficiently, and they remain for longer inside the chromospheric cavity. The power spectrum exhibits a tail with high power starting from the resonance power peak towards higher frequencies. This power enhancement is more prominent for the models with the transition region located at lower heights. This result is in agreement with the simulations from Snow et al. (2015).

As a reference, the bottom panel of Fig. 6 exhibits the power spectrum obtained from a simulation with the same driver (with its maximum power at a period of five minutes) but without a chromospheric cavity formed by the presence of the transition region (Simulation 3.3.6). In this case, the power is also shifted towards the three-minute band. However, there are some significant differences. First, the power peak is located at around 6 mHz, just above the local maximum of the acoustic cutoff frequency. As discussed previously, this change in the dominant frequency is due to the linear propagation of waves, since the amplitude increase of waves above the acoustic cutoff is higher than that of evanescent waves. Second, the spectrum is broader and there are no sharp peaks as a result of the resonance. Third, the power is almost two orders of magnitude lower than in some simulations with a chromospheric cavity. In this simulation the high-frequency waves travel through the domain and escape at the top boundary of the computational domain, so the wave power after the initial time steps is negligible. In addition, the resonance produces a remarkable enhancement of the power at a specific period, which makes the contrast even stronger.

3.4. Magnetized atmospheres

In this set of simulations (see Table 4), the atmospheric model with the base of the corona located at z = 2.0 Mm is permeated with a homogeneous vertical magnetic field. A different value of the magnetic field is chosen for each simulation. In order to avoid the small temporal step imposed by the extremely high coronal Alfvén speed, we constrained the strength of the Lorentz force following Rempel et al. (2009). The Alfvén speed is artificially limited to 10 000 km s⁻¹. At chromospheric layers and up to the transition region, the reduction of the Lorentz force is negligible. The whole chromospheric cavity is not affected by this Alfvén speed limiter. For one of the cases with strong magnetic field, we performed another simulation with the Alfvén speed limited to 20 000 km s⁻¹ (instead of 10 000 km s⁻¹), and we found no significant differences in the results.

The presence of the magnetic field changes the nature of the waves. Whereas in previous experiments, with the magnetic field set to zero, we analyzed the behavior of acoustic waves, in this section we study the effect of the chromospheric cavity on new wave modes. The theoretical picture for a homogeneous atmosphere with constant magnetic field indicates that three wave modes are present: fast and slow magnetoacoustic waves and the Alfvén wave. In a inhomogenous atmosphere (like those analyzed in this study), the three wave modes are not completely decoupled, but still it is common to refer to these terms to discuss the waves. We are going to follow this approach. It is a good approximation in the regions where the Alfvén speed is much higher than the sound speed (or vice versa), but not so good around $\beta =c_{\text{s}}^2/v_{\text{A}}^2=1$.

In our simulations, the driver is located in a high-β region (where the sound speed is much higher than the Alfvén speed). The vertical velocity introduced by the driver generates a fast magnetoacoustic wave whose behavior is similar to an acoustic wave. At the layer where the Alfvén speed approaches the sound speed (around the photosphere, depending on the magnetic field strength), the incident fast wave can be converted into a fast or a slow magnetoacoustic mode. In all the simulations discussed here we have imposed a vertical magnetic field. Since the direction of propagation of the fast wave (in the high-β region) is aligned with the direction of the magnetic field, conversion from fast to slow magnetoacoustic mode is favored at the height where β ∼ 1. The slow magnetoacoustic wave in the low-β region behaves similarly to an acoustic wave, but it propagates along field lines.

Figure 7 shows the power spectra of the chromospheric velocity for various numerical simulations with different vertical magnetic field strengths, from 0 to 2500 G. The magnetic field produces small differences in the period of the power peak. For small values of the magnetic field strength (up to 1000 G), an increase in the field strength is associated to an increase in the period. For stronger magnetic field strengths, the period is reduced as higher values of the magnetic field strength are considered. According to Eq. (20), the period of the fundamental mode of the resonance is determined by the size of the cavity and the wave speed. The wave speed of slow magnetoacoustic waves propagating along field lines in the low-β region is the sound speed. However, part of the chromospheric cavity is located in a region where β ∼ 1, especially in the cases with low magnetic field strength. The wave speed in that region will take some value between the sound speed and the Alfvén speed. This way, in the cases with magnetic field, ${\overline{v}}_{\mathrm{c}}$ will be generally higher than in the field-free cases. According to Eq. (20), the increase of P_r with the magnetic field for low field strengths (up to 1000 G) must be due to an increase in the size of the cavity.

Fig. 7. Power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with waves excited by a broadband driver and various magnetic field strengths. The vertical magnetic field strength is 0 (black line, Simulation 3.4.1), 10 G (violet line, Simulation 3.4.2), 250 G (dark blue, Simulation 3.4.3), 500 G (light blue, Simulation 3.4.4), 1000 G (dark green, Simulation 3.4.5), 1500 G (light green, Simulation 3.4.6), 2000 G (orange, Simulation 3.4.7), and 2500 G (red line, Simulation 3.4.8). Solid lines indicate simulations with the base of the corona located at 2.0 Mm. Dashed lines indicate simulations without transition region or corona, with a vertical magnetic field strength of 0 (black line, bottom panel, 3.4.10) and 2500 G (red lines, both panels, Simulation 3.4.9). All the power spectra are normalized to the same value used for the normalization of the power spectra illustrated in Fig. 6.

The resonance in atmospheres with stronger magnetic field produces broader peaks with lower power. This indicates that the magnetic field enhances the leaking of the waves out of the chromospheric cavity. In order to test the capability of the waves to pass through the bottom boundary of the cavity (the local maximum in the magnetoacoustic cutoff frequency, Fig. 1), we compared two simulations without transition region or corona, one of them with a 2500 G magnetic field strength and the other field free. In these cases there is no transition region resonance, since waves escape through the top boundary. The chromospheric power is plotted in the bottom panel of Fig. 7. In the simulation with B = 2500 G the power is significantly higher, and the maximum power is located at a lower frequency. These results confirm that the magnetic field reduces the effective cutoff frequency, and low-frequency waves can travel more easily through the temperature minimum when the magnetic field is present. As the magnetic field strength increases, the filtering of the lower frequencies in the spectrum is reduced. This is in qualitative agreement with the lower cutoff frequency derived for magnetoacoustic waves in an atmosphere with a vertical magnetic field (Fig. 1). Thus, in strongly magnetized atmospheres, the chromospheric cavity behaves like a resonator with losses.

3.5. Effect of radiative losses

We explore the impact of the radiative losses on the power spectra obtained from the resonant cavity. We computed various simulations with differences in the characteristic radiative cooling time τ_R (see Table 5). The radiative energy losses were implemented following Newton’s cooling law:

$$\begin{array}{c}\hfill Q_{\mathrm{rad}}=-c_{\mathrm{v}}\frac{T_1}{{\tau }_{\mathrm{R}}},\end{array}$$(21)

where T₁ is the temperature perturbation and c_v is the specific heat at constant volume. For each simulation we set a constant value for τ_R above z = 0 (in the simulations discussed in previous sections τ_R → ∞). Radiative losses are neglected below z = 0. In this section we made no effort to reproduce realistic values of the radiative losses in the actual umbral atmosphere. We aim to perform a parametric study on the effect of constant radiative losses on the chromospheric power spectra.

Figure 8 illustrates the power spectra of the vertical velocity at z = 1.4 Mm for all the cases with constant τ_R above the photosphere. When the radiative losses are neglected (black solid line in top panel of Fig. 8), the power spectrum exhibits a sharp peak with a period of 184.3 s, as discussed in previous sections. The introduction of low radiative losses (high values of τ_R) produces a striking reduction in the strength of the power peak and a broadening of the power distribution. With the presence of radiative losses there is no strict acoustic cutoff, and all frequencies are partially reflected and partially propagating (Roberts 1983; Centeno et al. 2006; Khomenko et al. 2008). The ratio of propagation to reflection increases with frequency. In a similar way to the effects of the magnetic field on the chromospheric cavity, the radiative losses also increase the losses of the resonator, and waves can escape from the cavity.

Fig. 8. Power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with various values of the radiative cooling time. In all cases waves are excited by a broadband driver and the magnetic field is set to zero. Each panel illustrates a different set of simulations (we note the differences in the range of the vertical axis). From top to bottom, the radiative cooling time is infinity (black line, Simulation 3.5.1), 3000 s (dark red line, Simulation 3.5.11), 1200 s (red line, Simulation 3.5.10), 600 s (orange line, Simulation 3.5.9), 300 s (yellow line, Simulation 3.5.8), 180 s (light green line, Simulation 3.5.7), 100 s (green line, Simulation 3.5.6), 70 s (dark green line, Simulation 3.5.5), 40 s (light blue line, Simulation 3.5.4), 20 s (dark blue line, Simulation 3.5.3), and 10 s (violet line, Simulation 3.5.2). Solid lines indicate simulations with the base of the corona located at 2.0 Mm. Dashed lines indicate simulations without a transition region or corona, with a radiative cooling time of 10 s (violet dashed line, Simulation 3.5.12), 20 s (dark blue dashed line, Simulation 3.5.13), or infinity (black dashed line, Simulation 3.5.14). All the power spectra are normalized to the same value used for the normalization of the power spectra illustrated in Fig. 6.

The minimum power is found for τ_R = 40 s. There are several aspects that determine the power of the waves at chromospheric heights. First, lower radiative losses produce a better resonant cavity, and the power at the resonant frequency (around 184 s) is enhanced more efficiently. Second, higher radiative losses increase the damping of the waves, since they lose energy as they propagate from the photosphere to the chromosphere. Third, higher radiative losses reduce the effective acoustic cutoff frequency, allowing more wave energy to reach the chromosphere. This explains the increase of the maximum power in the cases with τ_R <  40 s (dark blue and violet lines in bottom panel of Fig. 8). Thanks to the reduction of the acoustic cutoff, waves with a period around 200 s can propagate from the photosphere to the chromosphere. Since the broadband driver introduces more wave energy at this period than at shorter periods, the chromospheric power of the cases with very low τ_R is stronger despite the enhanced radiative losses.

The power spectra as a function of the radiative cooling time exhibit several regimes. For the cases with longer radiative cooling time (above ∼180 s), the spectra are dominated by the resonant mode, with a narrow band at 184 s and a tail with high power towards shorter periods. The strength of the resonant peak decreases with shorter radiative cooling times. For those cases with very short radiative cooling times, the chromospheric spectra are not dominated by the resonance, but by the power of waves propagating from deeper layers up to the chromosphere instead. This produces a shift of the power to longer periods, since the acoustic cutoff frequency is reduced due to radiative losses.

In order to confirm this interpretation, we analyzed some cases without a transition region or corona. The dashed lines in the bottom panel of Fig. 8 show the power spectra for two simulations without a transition region and short radiative cooling times (dark blue and violet dashed lines, Simulations 3.5.12 and 3.5.13). The power distribution is similar to that obtained for the simulations with a transition region (dark blue and violet solid lines). However, when the transition region is present the power is stronger, since the sharp temperature gradient at the transition region acts as a barrier and partially reflects the waves. This way, more wave energy remains at the chromosphere. The dashed black line (bottom panel of Fig. 8) shows the power spectra from a simulation without radiative losses and without a transition region (Simulation 3.5.14). In this case the chromospheric power is concentrated at shorter periods than in the cases with radiative losses (dark blue and violet dashed lines), because in the latter the effective acoustic cutoff frequencies are reduced.

3.6. Chromospheric cavity with radiative losses and magnetic field

Following Spiegel (1957), the radiative cooling time can be obtained as

$$\begin{array}{c}\hfill {\tau }_{\mathrm{R}}=\frac{\rho c_{\mathrm{v}}}{16\chi {\sigma }_{\mathrm{R}}{T}^3},\end{array}$$(22)

where T is the temperature, χ is the gray absorption coefficient, and σ_R is the Stefan-Boltzmann constant. This expression is defined in the optically thin limit and assuming local thermodynamic equilibrium and, thus, it is only applicable at photospheric heights. Figure 9 shows the vertical stratification of τ_R in our model atmosphere. It takes values of the order of seconds at the photosphere and rapidly increases with height, showing values around τ_R ∼ 3500 s at z = 0.9 Mm. At higher layers τ_R decreases with height, with τ_R = 85 s at z = 1.15 Mm (where the temperature of the atmospheric model is 6360 K). This radiative cooling time is in agreement with the chromospheric value determined by Giovanelli (1978). Above that layer, Giovanelli (1978) estimated an increase of τ_R up to 390 s for T = 8400 K, whereas according to the formula of Spiegel (1957) the radiative cooling time continues to decrease. For simplicity, the radiative loss function is only applied in the region between z = 0.2 and z = 1.15 Mm. Out of this region the applicability of Eq. (22) is uncertain. Assuming adiabatic propagation is justified for chromospheric layers, that is, those with low optical depth (Schmieder 1977). Below z = 0.2 Mm, Spiegel’s formula exhibits unrealistically low values for the radiative cooling time due to the optically thin treatment. We followed the same approach from Ulmschneider (1971), assuming a completely optically thick atmosphere below a certain arbitrary height. Ulmschneider (1971) chose a height between z = 0.11 and z = 0.14 Mm as the lower limit of applicability of Eq. (22). We set it a z = 0.2 Mm. At that height, the minimum value of τ_R used in these simulations is found (τ_R = 24 s).

Fig. 9. Radiative cooling time in the model atmosphere according to Spiegel (1957).

Figure 10 illustrates the chromospheric power spectra measured from simulations with the radiative losses implemented following Eq. (22) (see Table 6). We compare nonmagnetized (red lines) and magnetized atmospheres (2500 G, blue lines) for the simulations with a chromospheric cavity (solid line) and without a transition region (dashed lines). A comparison of the two simulations with the magnetic field set to zero shows remarkable differences. When the waves are reflected at the top of the chromosphere by the steep temperature gradient (solid line), a resonance is formed with the power peak at a period of 183.8 s. On the contrary, if waves are able to escape through the top boundary, the chromospheric power is significantly lower, since it is only due to the linear propagation of high-frequency waves from deeper layers. In this case, the chromospheric power spectrum is dominated by waves with frequency above the effective magnetoacoustic cutoff frequency (taking into account the effect of the radiative losses) and exhibits a maximum value at a period of 171.5 s.

Fig. 10. Power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with the photospheric radiative losses estimated from Spiegel (1957). Red lines correspond to cases with the magnetic field set to zero, and blue lines to a magnetic field strength of 2500 G. Solid lines represent the cases with the base of the corona located at z = 2.0 Mm, and dashed lines to cases without a corona or transition region. All the power spectra are normalized to the same value used for the normalization of the power spectra illustrated in Fig. 6.

In the case of the simulations with the atmosphere permeated by a vertical magnetic field, the maximum of the power is obtained at 186.0 and 182.6 s for the simulations with and without a transition region, respectively. Significant differences are found in the power between both cases, although these differences are not as striking as in the unmagnetized cases. The combined action of radiative losses and magnetic field strength facilitates the leaking of the waves, and the resonance produces a broader power spectrum with lower power. With realistic photospheric radiative losses and realistic values of the umbral magnetic field, the resonance is less efficient.

4. Discussion

In this paper, we study the chromospheric power spectra from numerical simulations of wave propagation in a sunspot umbra. We explore the response of the umbral atmosphere to wave excitation in cases with different parameters. Our numerical experiments include changes in the driving mechanism, the height of the transition region and corona, the vertical magnetic field strength, and radiative losses.

One of the basic approaches of this paper is the comparison of the power spectra in simulations with a chromospheric cavity (when the transition region and corona are present) and without a resonant cavity. In the latter, waves escape through the top boundary of the computational domain, which is located at the upper chromosphere. These two numerical setups allow us to evaluate two mechanisms proposed to explain the three-minute oscillations observed at chromospheric heights. First, the action of a chromospheric resonant cavity (Zhugzhda & Locans 1981; Zhugzhda et al. 1983; Gurman & Leibacher 1984). Between the temperature gradient of the photosphere and that at the transition region, the acoustic (and magnetoacoustic) cutoff frequency exhibits a local minimum. Slow magnetoacoustic waves in a certain frequency range are partially trapped inside this region. The resonant modes are enhanced, while the frequencies outside the resonance range slowly leak out. Second, the linear propagation of waves with frequencies above the cutoff frequency. Their amplitude increases due to the drop of the density, and at the chromosphere they dominate over the evanescent low-frequency waves (Fleck & Schmitz 1991; Centeno et al. 2006).

Our simulations show that a monochromatic wave with a five-minute period traveling from the interior to the corona cannot generate the three-minute oscillations observed in umbral chromospheres. This finding contrasts with the recent conclusions from Kraśkiewicz et al. (2019; in ideal MHD) and Wójcik et al. (2018; using a two-fluid model), who report that the energy from a driving period of five minutes can be transferred to shorter periods as a result of wave reflection. Our simulations reproduce the same effect reported by Kraśkiewicz et al. (2019) and Wójcik et al. (2018). However, chromospheric three-minute waves are only found if the driver excites these frequencies. In our simulations, introducing a sinusoidal acoustic wave (including the perturbations in velocity, pressure, and density with the appropriate amplitude and phase relations) as a driver at the bottom boundary, we have proven that the oscillations in the three-minute band are only excited during the initial stages of the simulation. We interpret that they are the result of the small power introduced in the three-minute band by the driver with a five-minute period from the head of the wave during the transient at the initial times of the simulations. These three-minute waves are free to propagate upward, increasing their amplitude, and they initially dominate the chromospheric power spectra due to the spatial attenuation of the evanescent low-frequency waves (five-minute band), until they escape from the computational domain. If the transition region is present, a chromospheric cavity is formed, and the power of the three-minute waves is greatly enhanced by the resonance. This power is progressively reduced as the waves leak from the chromospheric cavity. When the stationary stage of the simulation is reached (after 230 min in the simulation using an acoustic wave as driver), we obtain a pure monochromatic driver and five-minute oscillations throughout the domain.

One fundamental difference of our numerical set up with those from Kraśkiewicz et al. (2019) and Wójcik et al. (2018) is the moving piston introduced at the bottom boundary. They used a sinusoidal wave in the vertical velocity, without including the perturbations in the rest of the variables with the amplitude and phase relations for an acoustic wave. As discussed by Fleck & Schmitz (1991), oscillations at the three-minute band are already present in this driver due to the relations between the thermodynamic quantities. As a sanity check, we repeated the simulations using as a driver a sinusoidal wave only in the vertical velocity. Those simulations require a much longer time to reach a stationary regime, and the chromospheric three-minute oscillations dominate for a long time. We can conclude that both mechanisms discussed in this paper for generating the three-minute power observed at the chromosphere (linear propagation of high-frequency waves and chromospheric resonant cavity) require the excitation of three-minute waves at lower heights.

We present a novel exploration of the impact on the resonant cavity of two key ingredients for the modeling of umbral atmospheres, namely the magnetic field and radiative losses. Until now, radiative losses have been neglected by analytical and numerical studies. With regards to the magnetic field, it was assumed that its only role was to serve as a waveguide for slow magnetoacoustic waves, and low magnetic field strengths were chosen (∼10 G, Botha et al. 2011; Snow et al. 2015). However, magnetic field is known to modify the cutoff frequency of magnetoacoustic waves (e.g., Thomas 1982, 1983; Campos 1986; Stark & Musielak 1993; Roberts 2006). Umbral atmospheres with realistic chromospheric magnetic field strength values (∼2000 G, e.g., de la Cruz Rodríguez et al. 2013) exhibit notable differences in the permeability of the bottom boundary of the resonant cavity. The power spectra of the resonant mode produced by magnetized atmospheres shows slight variations in the frequency of the power peak, and striking differences in the strength and broadening. The magnetic field enhances the wave leakage through the photospheric boundary of the resonant cavity, and it behaves like a resonator with significant losses. In this work we have focused on sunspot umbrae and have only analyzed vertical magnetic fields, but the magnetoacoustic cutoff frequency can also be lowered by inclined magnetic fields (Bel & Leroy 1977; Jefferies et al. 2006).

Radiative energy losses also produce a reduction of the cutoff frequency (Roberts 1983; Centeno et al. 2006; Khomenko et al. 2008). Their effect on the resonance is similar to that of the magnetic field, since they contribute to the leaking of the waves out of the resonant cavity. The changes in the cutoff frequencies are also evident in the simulations without a resonant cavity produced by the transition region. When the radiative losses and/or the magnetic field are present, the wave power that reaches the chromosphere is shifted to lower frequencies (longer periods).

5. Conclusions

The numerical investigation presented in this paper shows that two of the mechanisms proposed to explain the chromospheric three-minute oscillations (linear propagation of high-frequency waves and chromospheric resonant cavity) can actually reproduce the observed shift in the period of the wave power with height. We conclude that both mechanisms require the excitation of three-minute oscillations at deeper layers. Since waves are driven below the photosphere by turbulent plasma motions produced by magnetoconvection, it is expected that broadband excitation will always be present and, thus, both mechanisms can take place at the same time. The question remains of what is the relative importance of each of them. We have proven that the radiative losses and the magnetic field favor the leaking of the waves through the boundaries of the resonant cavity, and that they reduce the enhancement of the wave power due to the resonance. The simulations using realistic values for the magnetic field strength and the photospheric radiative relaxation time demonstrate that chromospheric three-minute oscillations are obtained without a chromospheric cavity, but when it is present their power is enhanced by almost a factor of three.

Acknowledgments

Financial support from the State Research Agency (AEI) of the Spanish Ministry of Science, Innovation and Universities (MCIU) and the European Regional Development Fund (FEDER) under grant with reference PGC2018-097611-A-I00 is gratefully acknowledged. Wavelet software was provided by C. Torrence and G. Compo and is available at http://paos.colorado.edu/research/wavelets/. The author wishes to acknowledge the contribution of Teide High-Performance Computing facilities to the results of this research. TeideHPC facilities are provided by the Instituto Tecnológico y de Energías Renovables (ITER, SA). URL: http://teidehpc.iter.es

References

All Tables

All Figures

Fig. 1. Temperature (top panel), acoustic cutoff frequency (middle panel), and density (bottom panel) of the sunspot models. They differ in the height of the base of the corona, which is located at 1.8 Mm (orange), 2.0 Mm (green), 2.2 Mm (light blue), 2.4 Mm (dark blue, and 2.6 Mm (violet). The black dashed line corresponds to the magnetoacoustic cutoff of the umbral model with the base of the corona at 2.0 Mm and a vertical magnetic field of 2500 G. The red line is the model without corona. The dotted black line corresponds to the sunspot model M from Maltby et al. (1986).

In the text

	Fig. 2. Temporal evolution of the vertical velocity at z = 0.3 Mm (top panel) and at z = 1.4 Mm (bottom panel) for Simulation 3.1.1.
In the text

Fig. 3. Top panel: power as a function of period (vertical axis) and time (horizontal axis) obtained from the wavelet analysis of the vertical velocity at z = 1.4 Mm from Simulation 3.1.1. The power is normalized and shown in a logarithmic scale. The results are reliable for the region within the two white lines. The vertical dashed lines indicate the time steps plotted in the middle panel. Middle panel: normalized power at t = 67 min (blue line) and at t = 430 min (red line) obtained from the wavelet analysis. Bottom panel: normalized power for the full temporal series (500 min) obtained from the FFT.

In the text

	Fig. 4. Same as Fig. 3 but for Simulation 3.1.2. In this simulation there is no transition region.
In the text

	Fig. 5. Power spectra of the vertical velocity from the stationary part of Simulation 3.2.1 at z = 0.3 Mm (blue line) and z = 1.4 Mm (red line). Both spectra are normalized to the maximum value at z = 1.4 Mm.
In the text

Fig. 6. Top panel: power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with waves excited by a broadband driver and with the base of the corona at different heights. The color code is the same as that in Fig. 1, and the base of the corona is located at 1.8 Mm (orange, Simulation 3.3.1), 2.0 Mm (green, Simulation 3.3.2), 2.2 Mm (light blue, Simulation 3.3.3), 2.4 Mm (dark blue, Simulation 3.3.4), and 2.6 Mm (violet, Simulation 3.3.5). Bottom panel: power spectra of the vertical velocity at z = 1.4 Mm from Simulation 3.3.6 using a broadband driver and without transition region and corona. All the power spectra are normalized to the maximum value of the power from Simulation 3.3.2 (height of the bottom of the corona at 2.0 Mm, nonmagnetic, and without radiative losses).

In the text

Fig. 7. Power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with waves excited by a broadband driver and various magnetic field strengths. The vertical magnetic field strength is 0 (black line, Simulation 3.4.1), 10 G (violet line, Simulation 3.4.2), 250 G (dark blue, Simulation 3.4.3), 500 G (light blue, Simulation 3.4.4), 1000 G (dark green, Simulation 3.4.5), 1500 G (light green, Simulation 3.4.6), 2000 G (orange, Simulation 3.4.7), and 2500 G (red line, Simulation 3.4.8). Solid lines indicate simulations with the base of the corona located at 2.0 Mm. Dashed lines indicate simulations without transition region or corona, with a vertical magnetic field strength of 0 (black line, bottom panel, 3.4.10) and 2500 G (red lines, both panels, Simulation 3.4.9). All the power spectra are normalized to the same value used for the normalization of the power spectra illustrated in Fig. 6.

In the text

Fig. 8. Power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with various values of the radiative cooling time. In all cases waves are excited by a broadband driver and the magnetic field is set to zero. Each panel illustrates a different set of simulations (we note the differences in the range of the vertical axis). From top to bottom, the radiative cooling time is infinity (black line, Simulation 3.5.1), 3000 s (dark red line, Simulation 3.5.11), 1200 s (red line, Simulation 3.5.10), 600 s (orange line, Simulation 3.5.9), 300 s (yellow line, Simulation 3.5.8), 180 s (light green line, Simulation 3.5.7), 100 s (green line, Simulation 3.5.6), 70 s (dark green line, Simulation 3.5.5), 40 s (light blue line, Simulation 3.5.4), 20 s (dark blue line, Simulation 3.5.3), and 10 s (violet line, Simulation 3.5.2). Solid lines indicate simulations with the base of the corona located at 2.0 Mm. Dashed lines indicate simulations without a transition region or corona, with a radiative cooling time of 10 s (violet dashed line, Simulation 3.5.12), 20 s (dark blue dashed line, Simulation 3.5.13), or infinity (black dashed line, Simulation 3.5.14). All the power spectra are normalized to the same value used for the normalization of the power spectra illustrated in Fig. 6.

In the text

	Fig. 9. Radiative cooling time in the model atmosphere according to Spiegel (1957).
In the text

Fig. 10. Power spectra of the vertical velocity at z = 1.4 Mm from numerical simulations with the photospheric radiative losses estimated from Spiegel (1957). Red lines correspond to cases with the magnetic field set to zero, and blue lines to a magnetic field strength of 2500 G. Solid lines represent the cases with the base of the corona located at z = 2.0 Mm, and dashed lines to cases without a corona or transition region. All the power spectra are normalized to the same value used for the normalization of the power spectra illustrated in Fig. 6.

In the text