Linear Acoustic Waves in a Nonisothermal Atmosphere. II. Photospheric Resonator Model of Three-minute Umbral Oscillations

We consider an atmosphere that is in hydrostatic equilibrium,

$$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}p}{{dz}}=-\displaystyle \frac{\gamma g}{{c}^{2}},\end{eqnarray} \tag{ 3 }$$

and consists of several layers. The gravitational acceleration g is assumed to be constant. As is well known, the sound speed c is related to pressure p, density ρ, and temperature T,

$$\begin{eqnarray}&&c=\sqrt{\displaystyle \frac{\gamma p}{\rho }}=\sqrt{\gamma {RT}},\end{eqnarray} \tag{ 4 }$$

with $R\equiv {k}_{B}/\mu {m}_{{\rm{H}}}$, where k_B is the Boltzmann constant, μ is the mean molecular weight, and m_H is the hydrogen mass.

One can freely choose the reference pressure p_r that specifies the reference atmospheric level, and the corresponding sound speed c_r at the reference level. By introducing the scale length

$$\begin{eqnarray}&&{H}_{r}\equiv \displaystyle \frac{2{c}_{r}^{2}}{\gamma g},\end{eqnarray} \tag{ 5 }$$

which is equal to twice the local pressure scale height at the reference level, we define the reduced height,

$$\begin{eqnarray}&&h\equiv -\displaystyle \frac{1}{2}{H}_{r}\mathrm{ln}\displaystyle \frac{p}{{p}_{r}},\end{eqnarray} \tag{ 6 }$$

which is measured from the reference level. Note that we measure the geometrical height z from this reference level as well. In the special case of an isothermal atmosphere, h is equal to z.

In Paper I we have considered the wave solution in a simple layer where dc²/dh is constant. We generalize this condition to the multilayer atmosphere. We require that it is constant in each layer indexed by l and is specified by the dimensionless parameter κ_l, which is defined as

$$\begin{eqnarray}&&{\kappa }_{l}\equiv \displaystyle \frac{{H}_{r}}{{c}_{l}^{2}}\displaystyle \frac{{{dc}}^{2}}{{dh}}\end{eqnarray} \tag{ 7 }$$

for h_l ≤ h < h_l + 1, where h_l refers to the bottom of the lth layer, and c_l is the sound speed there. This equation is integrated to give

$$\begin{eqnarray}&&\displaystyle \frac{{c}^{2}}{{c}_{l}^{2}}=1+{\kappa }_{l}\displaystyle \frac{h-{h}_{l}}{{H}_{r}}.\end{eqnarray} \tag{ 8 }$$

By combining Equations (3), (5), (6), and (8), we obtain the relationship between h and z:

$$\begin{eqnarray}&&\displaystyle \frac{{dz}}{{dh}}=\left(1+{\kappa }_{l}\displaystyle \frac{h-{h}_{l}}{{H}_{r}}\right)\displaystyle \frac{{c}_{l}^{2}}{{c}_{r}^{2}},\end{eqnarray} \tag{ 9 }$$

which is integrated to give

$$\begin{eqnarray}&&z-{z}_{l}=(h-{h}_{l})\left(1+\displaystyle \frac{1}{2}{\kappa }_{l}\displaystyle \frac{h-{h}_{l}}{{H}_{r}}\right)\displaystyle \frac{{c}_{l}^{2}}{{c}_{r}^{2}},\end{eqnarray} \tag{ 10 }$$

and, after inversion

$$\begin{eqnarray}&&h-{h}_{l}=\displaystyle \frac{{H}_{r}}{{\kappa }_{l}}\left(-1+\sqrt{1+2\displaystyle \frac{{\kappa }_{l}}{{H}_{r}}\displaystyle \frac{{c}_{r}^{2}}{{c}_{l}^{2}}(z-{z}_{l})}\right).\end{eqnarray} \tag{ 11 }$$

The general solution for the linear acoustic waves can be easily constructed using the principle of superposition. We express pressure fluctuation and velocity fluctuation for each layer indexed by l as

$$\begin{eqnarray}&&{\rm{\Delta }}p(h,t;l)={p}^{\tfrac{1}{2}}(h){\int }_{-\infty }^{\infty }{\psi }_{\omega }(h;l)\exp (-i\omega t)d\omega \end{eqnarray} \tag{ 12 }$$

$$\begin{eqnarray}&&v(h,t;l)={p}^{-\tfrac{1}{2}}(h){\int }_{-\infty }^{\infty }{\varphi }_{\omega }(h;l)\exp (-i\omega t)d\omega \end{eqnarray} \tag{ 13 }$$

using the complex wave functions ψ_ω (h; l) and φ_ω (h; l) that are written as

$$\begin{eqnarray}&&{\psi }_{\omega }(h;l)={C}_{l}\left({f}_{\omega }(h;l)+{r}_{\omega ,l}{f}_{\omega }^{* }(h;l)\right)\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&{\varphi }_{\omega }(h;l)={D}_{l}\left({g}_{\omega }(h;l)+{r}_{\omega ,l}{g}_{\omega }^{* }(h;l)\right),\end{eqnarray} \tag{ 15 }$$

where f_ω (h; l) and g_ω(h; l) in general represent the progressive components of pressure fluctuation and velocity fluctuation, respectively, and their conjugates ${f}_{\omega }^{* }(h;l)$ and ${g}_{\omega }^{* }(h;l)$, the regressive components. The reflection coefficient r_ω,l is the coefficient ratio of the regressive component to the progressive component.

Note that Equations (12) and (13) are similar to Equations (18) and (19) of Paper I, but are different from these in that Equations (12) and (13) are the full Fourier transform representation, while Equations (18) and (19) of Paper I are a simple summation representation. This difference mostly accounts for the difference in the formulae and expressions between the current work and our previous one in Paper I.

According to Paper I, in the lth layer, the coefficients C_l and D_l are related to each other as

$$\begin{eqnarray}&&\displaystyle \frac{{D}_{l}}{{C}_{l}}=\displaystyle \frac{{ig}}{2\omega }.\end{eqnarray} \tag{ 16 }$$

Following Paper I, we change the independent variable h in the lth layer into a new dimensionless variable ζ that is related to ω, h and l,

$$\begin{eqnarray}&&\zeta \equiv \displaystyle \frac{1}{{a}_{l}^{2}}\left(\displaystyle \frac{{\omega }^{2}}{{\omega }_{l}^{2}}-1\right)+\displaystyle \frac{{a}_{l}}{{H}_{r}}(h-{h}_{l}),\end{eqnarray} \tag{ 17 }$$

with the Lamb frequency at the bottom of the lth layer

$$\begin{eqnarray}&&{\omega }_{l}=\displaystyle \frac{\gamma g}{2{c}_{l}}\end{eqnarray} \tag{ 18 }$$

and a dimensionless parameter a_l defined as

$$\begin{eqnarray}&&{a}_{l}={s}_{l}{\left(| {\kappa }_{l}| \displaystyle \frac{{\omega }^{2}}{{\omega }_{l}^{2}}\right)}^{1/3},\end{eqnarray} \tag{ 19 }$$

where s_l is the sign of κ_l. Note that a_l depends on ω as well as κ_l and ω_l. It is obvious that the value of ζ at the bottom of the lth layer (h = h_l) is given by

$$\begin{eqnarray}&&\zeta (\omega ,{h}_{l};l)=\displaystyle \frac{1}{{a}_{l}^{2}}\left(\displaystyle \frac{{\omega }^{2}}{{\omega }_{l}^{2}}-1\right),\end{eqnarray} \tag{ 20 }$$

and after some arrangements, one can show that ζ at the top of the same layer is given by

$$\begin{eqnarray}&&\zeta (\omega ,{h}_{l+1};l)=\displaystyle \frac{1}{{a}_{l}^{2}}\left(\displaystyle \frac{{\omega }^{2}}{{\omega }_{l+1}^{2}}-1\right).\end{eqnarray} \tag{ 21 }$$

Note that h is continuous across the interface between the lth layer and l + 1th layer, but ζ is not in general: $\zeta (\omega ,{h}_{l+1};l)\ne \zeta (\omega ,{h}_{l+1};l+1)$ unless a_l = a_l + 1.

It has been shown in Paper I that the equation describing the function f(ζ) = f_ω(h; l) is given by

$$\begin{eqnarray}&&\displaystyle \frac{{d}^{2}f}{d{\zeta }^{2}}+\zeta f=0\,\end{eqnarray} \tag{ 22 }$$

and the progressive solution representing the waves propagating upwards is given by

$$\begin{eqnarray}&&f(\zeta )={f}_{J}(\zeta )+{{is}}_{l}{f}_{Y}(\zeta ),\end{eqnarray} \tag{ 23 }$$

with two real function f_J and f_Y that can be expressed in terms of Bessel functions and modified functions. These two functions can also be expressed in terms of the Airy functions:

$$\begin{eqnarray}&&{f}_{J}(\zeta )=\displaystyle \frac{3}{2}\mathrm{Ai}(-\zeta )-\displaystyle \frac{\sqrt{3}}{2}\mathrm{Bi}(-\zeta )\end{eqnarray} \tag{ 24 }$$

$$\begin{eqnarray}&&{f}_{Y}(\zeta )=-\displaystyle \frac{\sqrt{3}}{2}\mathrm{Ai}(-\zeta )-\displaystyle \frac{3}{2}\mathrm{Bi}(-\zeta ).\end{eqnarray} \tag{ 25 }$$

Note that f_J is a real function that behaves like a sine function at ζ > 0, and f_Y is a real function that behaves like a negative cosine function (see Figure 4 of Paper I). Once f(ζ) is determined, g_ω(h; l) = g(ζ) can be determined from f_ω using the equation

$$\begin{eqnarray}&&g(\zeta )=f(\zeta )-{a}_{l}\displaystyle \frac{{df}(\zeta )}{d\zeta }.\end{eqnarray} \tag{ 26 }$$

The behavior of the functions f(ζ) and g(ζ) at three different regimes can be conveniently investigated using the asymptotic expressions

$$\begin{eqnarray}f(\zeta )\approx \left\{\begin{array}{ll}\sqrt{\tfrac{3}{\pi }}{\zeta }^{-\tfrac{1}{4}}{e}^{{{is}}_{l}\theta } & \mathrm{for}\zeta \gt 1\\ 0.7765\zeta +{s}_{l}(0.4483\zeta -1.2299)i & \mathrm{for}-1\lt \zeta \lt 1\\ \sqrt{\tfrac{3}{\pi }}| \zeta {| }^{-\tfrac{1}{4}}\exp (\tfrac{2}{3}| \zeta {| }^{\tfrac{3}{2}}){e}^{-(2\pi /3){s}_{l}i} & \mathrm{for}\zeta \lt -1\end{array}\right.\end{eqnarray} \tag{ 27 }$$

and

$$\begin{eqnarray}g(\zeta )\approx \left\{\begin{array}{ll}\sqrt{\tfrac{3}{\pi }}{\zeta }^{-\tfrac{1}{4}}{e}^{{{is}}_{l}\theta }\left(1-i| {a}_{l}| {\zeta }^{1/2})\right) & \mathrm{for}\zeta \gt 1\\ 0.7765(\zeta -{a}_{l})+{s}_{l}(0.4483(\zeta -{a}_{l})-1.2299)i & \mathrm{for}-1\lt \zeta \lt 1\\ \sqrt{\tfrac{3}{\pi }}| \zeta {| }^{-\tfrac{1}{4}}\exp (\tfrac{2}{3}| \zeta {| }^{\tfrac{3}{2}}){e}^{-(2\pi /3){s}_{l}i}(1+{a}_{l}| \zeta {| }^{1/2}) & \mathrm{for}\zeta \lt -1\end{array}\right.,\end{eqnarray} \tag{ 28 }$$

with

$$\begin{eqnarray}&&\theta \equiv \displaystyle \frac{2}{3}{\zeta }^{\tfrac{3}{2}}-\displaystyle \frac{5\pi }{12}.\end{eqnarray} \tag{ 29 }$$

We consider the atmosphere composed of L layers. At the interface h = h_l+1 between the lth layer and the l + 1th layer, both ψ_ω and φ_ω should be continuous:

$$\begin{eqnarray}&&\begin{array}{l}{C}_{l}\left({f}_{\omega }({h}_{l+1};l)+{r}_{\omega ,l}{f}_{\omega }^{* }({h}_{l+1};l)\right)\\ \ \ \ =\,{C}_{l+1}\left({f}_{\omega }({h}_{l+1};l+1)+{r}_{\omega ,l+1}{f}_{\omega }^{* }({h}_{l+1};l+1)\right)\end{array}\end{eqnarray} \tag{ 30 }$$

$$\begin{eqnarray}&&\begin{array}{l}{D}_{l}\left({g}_{\omega }({h}_{l+1};l)+{r}_{\omega ,l}{g}_{\omega }^{* }({h}_{l+1};l)\right)\\ \ \ =\,{D}_{l+1}\left({g}_{\omega }({h}_{l+1};l+1)+{r}_{\omega ,l+1}{g}_{\omega }^{* }({h}_{l+1};l+1)\right)\end{array}\end{eqnarray} \tag{ 31 }$$

for l = 0, .., L − 2. We also impose the top boundary condition that only progressive waves exist in the top layer indexed by l = L − 1:

$$\begin{eqnarray}&&{r}_{\omega ,L-1}=0,\end{eqnarray} \tag{ 32 }$$

and the bottom boundary condition that the velocity at the bottom is given by

$$\begin{eqnarray}&&{\varphi }_{\omega }({h}_{0};0)={D}_{0}\left({g}_{\omega }({h}_{0};0)+{r}_{\omega ,0}{g}_{\omega }^{* }({h}_{0};0)\right),\end{eqnarray} \tag{ 33 }$$

where φ_ω (h₀; 0) is related to v (h₀, t; 0) by the equation

$$\begin{eqnarray}&&{\varphi }_{\omega }({h}_{0};0)=\displaystyle \frac{1}{2\pi }{p}_{0}^{\tfrac{1}{2}}{\int }_{-\infty }^{\infty }v({h}_{0},t;0)\exp (i\omega t){dt},\end{eqnarray} \tag{ 34 }$$

which has been derived from Equation (13).

In order to model a pulse of finite duration, we consider the propagation of the wave packet produced by the vertical displacement of the form

$$\begin{eqnarray}&&\xi ({h}_{0},t)=-{\xi }_{0}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{t-{t}_{0}}{{t}_{p}}\right)}^{2}\right],\end{eqnarray} \tag{ 35 }$$

corresponding to the velocity of a piston

$$\begin{eqnarray}&&v({h}_{0},t)=\displaystyle \frac{{\xi }_{0}}{{t}_{p}}\displaystyle \frac{t-{t}_{0}}{{t}_{p}}\exp \left[-\displaystyle \frac{1}{2}{\left(\displaystyle \frac{t-{t}_{0}}{{t}_{p}}\right)}^{2}\right]\end{eqnarray} \tag{ 36 }$$

at the low boundary. According to this velocity model, the piston is at rest at $t=-\infty$, is accelerated downward, and reaches the maximum displacement ξ₀ at t = t₀ with zero velocity. After this, the piston is accelerated upward, and reaches the original equilibrium position at $t=\infty$. The disturbance is characterized by the maximum displacement ξ₀ and the dynamic timescale of the pulse t_p.

By inserting this model into Equation (34), we obtain the expression of φ_ω (h₀; 0)

$$\begin{eqnarray}&&{\varphi }_{\omega }({h}_{0};0)=i\exp (i\omega {t}_{0})\displaystyle \frac{{\xi }_{0}{p}_{0}^{1/2}}{\sqrt{2\pi }}{t}_{p}\omega \exp \left(-\displaystyle \frac{1}{2}{\omega }^{2}{t}_{p}^{2}\right).\end{eqnarray} \tag{ 37 }$$

Figure 3 presents the plot of ${c}^{2}/{c}_{r}^{2}$ versus h/H_r constructed from the values of pressure and density in the mean umbral model of Maltby et al. (1986). Here we set the reference level of h = z = 0 at the surface where the optical depth at 500 nm is unity. The normalizing parameters c_r and H_r were calculated with pressure and density at this reference level.

Zoom In Zoom Out Reset image size

The figure clearly indicates that c² rapidly decreases with h in the photosphere, and then very slowly increases with h in the chromosphere. We find that this variation can be fairly well approximated by a two-layer model K₀ = −1.44 and K₁ = 0.04 with the interface of ${c}_{1}^{2}/{c}_{r}^{2}=0.8$ being at h₁/ H_r = 0.2, where K_l is the dimensionless gradient in the lth layer, defined by

$$\begin{eqnarray}&&{K}_{l}\equiv \displaystyle \frac{{H}_{r}}{{c}_{r}^{2}}\displaystyle \frac{{{dc}}^{2}}{{dh}}.\end{eqnarray} \tag{ 38 }$$

From Equations (7) and (38), we obtain the relationship

$$\begin{eqnarray}&&{\kappa }_{l}=\displaystyle \frac{{c}_{r}^{2}}{{c}_{l}^{2}}{K}_{l},\end{eqnarray} \tag{ 39 }$$

which indicates κ₁ = K₁/0.8 = 0.05.

The bottom of the photospheric layer is regarded as a free parameter of the model. To specify it, we introduce the dimensionless parameter

$$\begin{eqnarray}&&\tau \equiv \displaystyle \frac{{c}_{0}^{2}}{{c}_{1}^{2}},\end{eqnarray} \tag{ 40 }$$

which is equal to the ratio of temperature T₀/T₁ when γ R is constant. Once τ is known, we can calculate c₀, and hence κ₀ using Equation (39). Then by inserting l = 0, c = c₁, and h = h₁ into Equation (8) and by solving it for h₀, we obtain

$$\begin{eqnarray}&&{h}_{0}={h}_{1}-\left(\displaystyle \frac{{c}_{1}^{2}}{{c}_{0}^{2}}-1\right)\displaystyle \frac{{H}_{r}}{{\kappa }_{0}}.\end{eqnarray} \tag{ 41 }$$

In addition, by inserting l = 0, h = z = 0 into Equation (10) and by solving it for z₀, we obtain the expression

$$\begin{eqnarray}&&{z}_{0}={h}_{0}\left(1-\displaystyle \frac{1}{2}{\kappa }_{0}\displaystyle \frac{{h}_{0}}{{H}_{r}}\right)\displaystyle \frac{{c}_{0}^{2}}{{c}_{r}^{2}},\end{eqnarray} \tag{ 42 }$$

and by inserting l = 0, h = h₁, z = z₁ into Equation (10) and by solving it for z₁, we obtain the expression

$$\begin{eqnarray}&&{z}_{1}={z}_{0}+({h}_{1}-{h}_{0})\left(1+\displaystyle \frac{1}{2}{\kappa }_{0}\displaystyle \frac{{h}_{1}-{h}_{0}}{{H}_{r}}\right)\displaystyle \frac{{c}_{0}^{2}}{{c}_{r}^{2}}.\end{eqnarray} \tag{ 43 }$$

In our two-layer model, the chromospheric layer is a simple nonisothermal layer where waves propagate upward. There are no regressive waves in this layer. In Paper I we have considered the acoustic transmission efficiency of a simple nonisothermal layer where only progressive waves exist. The transmission efficiency T_ω of such a layer was defined by the ratio of the wave flux to the maximum value that can be realized when pressure and velocity fluctuations are perfectly correlated. Because the upper chromospheric layer indexed by l = 1 is a temperature-increasing layer where only progressive waves exist, its acoustic cutoff frequency is given by

$$\begin{eqnarray}&&{\omega }_{\mathrm{ac}}=\displaystyle \frac{{\omega }_{1}}{\sqrt{1-{\kappa }_{1}}},\end{eqnarray} \tag{ 44 }$$

which is equal to 1.03ω₁. According to Paper I, the transmission efficiency of this layer is given by

$$\begin{eqnarray}&&{T}_{\omega }=\displaystyle \frac{3}{\pi }\displaystyle \frac{| {a}_{1}| }{| {f}_{\omega }({h}_{1};1)| \cdot | {g}_{\omega }({h}_{1};1)| }.\end{eqnarray} \tag{ 45 }$$

Figure 4 presents the plot of T_ω versus ω. We find from the figure that there is a smooth transition of T_ω from zero to unity at frequencies around ω_ac and the waves transport energy even at frequencies lower than ω_ac. This is the property of the acoustic cutoff in a nonisothermal layer, which was reported in Paper I.

Zoom In Zoom Out Reset image size

Note that T_ω is the property of the layer. It is wholly determined by the parameters ω₁ and κ₁ of the upper layer, and is not at all affected by the parameters of the lower layer. In contrast, the reflection coefficient to be discussed in the following is the property of the interface, characterized by the difference in the parameters between the lower layer and the upper layer.

In the case of the two-layer model with L = 2, one can derive from Equations (30)–(32) the expression of the reflection coefficient r_ω,0 in the lower layer:

$$\begin{eqnarray}&&{r}_{\omega ,0}=\displaystyle \frac{{f}_{\omega }({h}_{1};1){g}_{\omega }({h}_{1};0)-{f}_{\omega }({h}_{1};0){g}_{\omega }({h}_{1};1)}{{f}_{\omega }^{* }({h}_{1};0){g}_{\omega }({h}_{1};1)-{f}_{\omega }({h}_{1};1){g}_{\omega }^{* }({h}_{1};0)},\end{eqnarray} \tag{ 46 }$$

with the magnitude and phase defined by

$$\begin{eqnarray}&&{r}_{\omega ,0}=| {r}_{\omega ,0}| {e}^{i{\theta }_{r}}.\end{eqnarray} \tag{ 47 }$$

Figure 5 presents the plots of r_ω,0 computed using Equation (46) as a function of ω. The plots indicate that $| {r}_{\omega ,0}|$ rapidly decreases with ω at frequencies around ω₁. We find from the figure that $| {r}_{\omega ,0}|$ is anticorrelated with T_ω. Note that r_ω,0 depends on the parameters of the lower layer but T_ω does not. Thus, the strong anticorrelation between $| {r}_{\omega ,0}|$ and T_ω may be understood only if they are very sensitive to the parameters that commonly affect both. The comparison of Figures 4 and 5 suggests that both the parameters are very much affected by the acoustic cutoff frequency of the chromospheric layer ω_ac , which is practically determined by the minimum temperature of the atmosphere.

Zoom In Zoom Out Reset image size

To understand how r_ω,0 depends on the parameters of the two layers, we examine Equation (46) in the special case of ω = ω₁. This choice leads to the following expressions: ζ(ω₁, h₁; 0) = ζ(ω₁, h₁; 1) = 0, ${a}_{0}=-{(| {\kappa }_{0}| \tau )}^{1/3}$, and ${a}_{1}=| {\kappa }_{1}{| }^{1/3}$. It then follows that f_ω1 (h₁; 0) = 1.2299i, f_ω1 (h₁; 1) = −1.2299i, ${g}_{{\omega }_{1}}({h}_{1};0)=1.2299i+(-0.7765+0.4483i){a}_{0}$, and ${g}_{{\omega }_{1}}({h}_{1};1)=-1.2299i+(-0.7765-0.4483i){a}_{1}$. After some algebra, we obtain the expression

$$\begin{eqnarray}&&{r}_{{\omega }_{1},0}=\displaystyle \frac{-(1-q)/(1+q)+0.5773i}{1+0.5773i},\end{eqnarray} \tag{ 48 }$$

with the parameter

$$\begin{eqnarray}&&q\equiv -\displaystyle \frac{{a}_{1}}{{a}_{0}}={\left(\displaystyle \frac{| {\kappa }_{1}| }{| {\kappa }_{0}| \tau }\right)}^{1/3}={\left(\displaystyle \frac{| {K}_{1}| }{| {K}_{0}| }\right)}^{1/3}.\end{eqnarray} \tag{ 49 }$$

Figure 6 presents the plots of the magnitude $| {r}_{\omega ,0}|$ and phase θ_r of c calculated with ω = ω₁ is calculated as functions of q using Equation (48). The figure clearly indicates that both $| {r}_{\omega ,0}|$ and θ_r decrease with q. The lower the value of q, the stronger the reflection. This means that for a fixed K₁, a high value of $| {K}_{0}|$ causes strong reflection at the interface. Namely, the more rapidly the temperature decreases with height in the photosphere, the more strongly are the waves reflected at the temperature minimum.

Zoom In Zoom Out Reset image size

An important point is that q only weakly depends on K₁/K₀. In the specific model we consider, K₀ = −1.44 and K₁ = 0.04, so that we have q = 0.30. Even if K₁/K₀ were to increase or decrease by a factor of 2, the range of q would be confined to the range from 0.24 to 0.38, in which case $| {r}_{{\omega }_{1},0}|$ would be still in the range from 0.63 to 0.73, little deviating from the value of 0.68 calculated with the specific model. Therefore we conclude that the reflection coefficient is only weakly affected by the specific structure of the atmospheric model. It is primarily determined by ω₁.

It is interesting that the reflection coefficient ${r}_{{\omega }_{1},0}$ is completely determined by one parameter q or by the ratio K₁/K₀. From this result we expect that r_ω,0 at other frequencies is also controlled by K₁/K₀. This can be explicitly demonstrated for ω close to ω₁ when the asymptotic formula for f(ζ) yields a correction to the reflection coefficient r_ω,0, which is on the order of $({\omega }^{2}/{\omega }_{1}^{2}-1)q$. Thus the dependence of r_ω,0 on ${a}_{0}^{2}$ and ${a}_{1}^{2}$ is indeed relatively weak. It is somewhat simpler to analyze the dependence of r_ω,0 on ω if it is expressed in terms of the function f(ζ) and its derivatives:

$$\begin{eqnarray}&&{r}_{\omega ,0}=-\displaystyle \frac{{f}_{\omega }({h}_{1};1){f}_{\omega }^{{\prime} }({h}_{1};0)+{{qf}}_{\omega }({h}_{1};0){f}_{\omega }^{{\prime} }({h}_{1};1)}{{f}_{\omega }({h}_{1};1){f}_{\omega }^{* ^{\prime} }({h}_{1};0)+{{qf}}_{\omega }^{* }({h}_{1};0){f}_{\omega }^{{\prime} }({h}_{1};1)}.\end{eqnarray} \tag{ 50 }$$

This form makes it obvious, for instance, that $| {r}_{\omega ,0}| =1$ when q = 0, as expected physically.

The reflection at the temperature minimum produces the regressive waves, and the interference between the regressive waves and the progressive waves at the bottom of the photospheric layer results in a resonance at a specific frequency. Thus the photospheric layer plays the role of a resonator. We assume that the bottom of the photospheric resonator is subject to a piston driving, and its velocity can be set to a prescribed complex function φ_ω (h₀; 0) as in Equation (37). Then from Equation (33), one can derive the amplitude coefficient D₀:

$$\begin{eqnarray}&&{D}_{0}={D}_{\mathrm{0,0}}\displaystyle \frac{1}{1+| {r}_{\omega ,0}| {e}^{i({\theta }_{r}-2{\theta }_{0})}},\end{eqnarray} \tag{ 51 }$$

where D_0,0 is the value of D₀ in the limiting case of no reflection given by

$$\begin{eqnarray}&&{D}_{\mathrm{0,0}}\equiv \displaystyle \frac{{\varphi }_{\omega }({h}_{0};0)}{{g}_{\omega }({h}_{0};0)},\end{eqnarray} \tag{ 52 }$$

and θ₀ is the phase angle of g_ω (h₀; 0) defined by

$$\begin{eqnarray}&&{g}_{\omega }({h}_{0};0)=| {g}_{\omega }({h}_{0};0)| {e}^{i{\theta }_{0}}.\end{eqnarray} \tag{ 53 }$$

Therefore

$$\begin{eqnarray}&&{\left|\displaystyle \frac{{D}_{0}}{{D}_{\mathrm{0,0}}}\right|}^{2}=\displaystyle \frac{1}{1+| {r}_{\omega ,0}{| }^{2}+2| {r}_{\omega ,0}| \cos ({\theta }_{r}-2{\theta }_{0})},\end{eqnarray} \tag{ 54 }$$

which indicates that if $| {r}_{\omega ,0}|$ is not zero and is constant over ω, the amplitude $| {D}_{0},/{D}_{\mathrm{0,0}}{| }^{2}$ has a local peak at the frequency ω_M satisfying the condition

$$\begin{eqnarray}&&{\theta }_{r}-2{\theta }_{0}=\pi ,\end{eqnarray} \tag{ 55 }$$

with which it follows cos(θ_r − 2θ₀) = −1. This peak represents the resonance of the waves resulting from the interference between the progressive waves and regressive waves. Note that because $| {r}_{\omega ,0}|$ in fact varies with ω, the resonance frequency of Equation (54), denoted by ω_R, will actually differ from ω_M.

Figure 7 presents the plots of $| {D}_{0},/{D}_{\mathrm{0,0}}{| }^{2}$ as a function of ω and indicates the values of ω_R and ω_M in the three different cases of τ. Each resonance frequency ω_R lies lower than the corresponding ω_M because $| {r}_{\omega ,0}|$ decreases with ω. The figure clearly shows that the higher value of τ corresponds to a lower value of ω_R and a sharper peak of $| {D}_{0},/{D}_{\mathrm{0,0}}{| }^{2}$. The values of ω_R/ω₁ are 1.32, 1.07, and 0.98 for τ = 2.0, 2.5, and 3.0, respectively. It also should be noted that the resonance frequency ω_R has a value close to ω₁ in all cases.

Zoom In Zoom Out Reset image size

The time-averaged energy of the waves ${ \mathcal E }(h;l)$ at a height h is given by

$$\begin{eqnarray}&&{E}(h;l)=\rho (h;l)\displaystyle \frac{1}{{t}_{\mathrm{wv}}}{\int }_{-\infty }^{\infty }{v}^{2}(h,t;l){dt},\end{eqnarray} \tag{ 56 }$$

where t_wv is the duration of the wave packet or the length of the time interval where the velocity fluctuation has an amplitude significantly above zero. In the same way, the time-averaged flux of the wave energy is given by

$$\begin{eqnarray}&&{ \mathcal F }\equiv \displaystyle \frac{1}{{t}_{\mathrm{wv}}}{\int }_{-\infty }^{\infty }v{\rm{\Delta }}p\,{dt}.\end{eqnarray} \tag{ 57 }$$

By applying Parseval's theorem and making use of the fact that the amplitude is an even function of ω, we obtain

$$\begin{eqnarray}&&{E}(h;l)=\rho \displaystyle \frac{1}{{t}_{\mathrm{wv}}}{\int }_{0}^{\infty }\displaystyle \frac{4\pi }{p}| {\varphi }_{\omega }(h;l){| }^{2}d\omega \equiv {\int }_{0}^{\infty }{E}_{\omega }(h;l)d\omega \end{eqnarray} \tag{ 58 }$$

and

$$\begin{eqnarray}&&{ \mathcal F }=\displaystyle \frac{1}{{t}_{\mathrm{wv}}}{\int }_{0}^{\infty }4\pi {\mathfrak{R}}\left[{\varphi }_{\omega }{\psi }_{\omega }^{* }\right]d\omega \equiv {\int }_{0}^{\infty }{F}_{\omega }(h;l)d\omega .\end{eqnarray} \tag{ 59 }$$

These expressions lead us to define the energy spectrum, velocity power spectrum, and flux spectrum as follows:

$$\begin{eqnarray}&&{E}_{\omega }(h;l)=\displaystyle \frac{4\pi \rho }{{{pt}}_{\mathrm{wv}}}| {\varphi }_{\omega }(h;l){| }^{2}\end{eqnarray} \tag{ 60 }$$

$$\begin{eqnarray}&&{P}_{\omega }(h;l)=\displaystyle \frac{{E}_{\omega }(h;l)}{\rho }\end{eqnarray} \tag{ 61 }$$

$$\begin{eqnarray}&&{F}_{\omega }=\displaystyle \frac{4\pi }{{t}_{\mathrm{wv}}}{\mathfrak{R}}\left[{\varphi }_{\omega }{\psi }_{\omega }^{* }\right].\end{eqnarray} \tag{ 62 }$$

By inserting Equations (14) and (15), and making use of Equation (16), we obtain the expressions

$$\begin{eqnarray}&&{t}_{\mathrm{wv}}{E}_{\omega }(h;l)=\displaystyle \frac{4\pi \rho }{p}| {D}_{l}{| }^{2}| {g}_{\omega }+{r}_{\omega ,l}{g}_{\omega }^{* }{| }^{2}\end{eqnarray} \tag{ 63 }$$

$$\begin{eqnarray}&&{t}_{\mathrm{wv}}{F}_{\omega }=\displaystyle \frac{24\omega }{g}| {D}_{l}{| }^{2}(1-| {r}_{\omega ,l}{| }^{2})| {a}_{l}| .\end{eqnarray} \tag{ 64 }$$

This formula for flux spectrum F_ω is derived making use of the relationship³

$$\begin{eqnarray}&&{f}_{J}{g}_{Y}-{f}_{Y}{g}_{J}=-\displaystyle \frac{3}{\pi }{a}_{l},\end{eqnarray} \tag{ 65 }$$

which results from the Wronskian formulae for the Airy functions.

Note that E_ω in Equation (63) depends on h and l. Figure 8 shows the plots of E_ω calculated using Equation (63) at different heights in the specific case of τ = 2.5. We find from Figure 8 that a slight enhancement of E_ω(h)/E_ω(h₀) over unity is detected around the resonance frequency ω_R near the bottom of the photospheric resonator. The strength of the energy enhancement increases with height, and reaches a peak at the top of the resonator h/H_r = 0.20. The strength decreases with height at the upper heights h/H_r > 0.2. Thus the photospheric resonator has a node at its bottom end (h/H_r = −0.63), and an antinode at its top (h/H_r = 0.20). Figure 8 also shows that the energy spectra are affected by the acoustic cutoff. At low frequencies ω, E_ω(h)/E_ω(h₀) decreases with height and becomes much smaller than unity at heights h/H_r > 0.20. As a consequence of this cutoff, the profile of E_ω(h)/E_ω(h₀) at a large height, i.e., h/H_r = 2.52, appears asymmetric over ω, rapidly increasing with ω for ω ≲ ω₁ and slowly decreasing with ω for ω > ω₁.

Zoom In Zoom Out Reset image size

Figure 9 presents the plots of F_ω calculated using Equation (64). Note that F_ω is constant over h in the lth layer. It is in fact independent of l as well, as expected from the principle of energy conservation. The reference flux spectrum F_ω,0 used in the plots is the flux spectrum of acoustic waves in a homogeneous medium:

$$\begin{eqnarray}&&{F}_{\omega ,0}\equiv {c}_{0}{E}_{\omega }({h}_{0};0).\end{eqnarray} \tag{ 66 }$$

Zoom In Zoom Out Reset image size

The plot of F_ω/F_ω,0 in the case of τ = 2.5 is qualitatively similar to the plot of E_ω(h)/E_ω(h₀) at the height (h−h₁)/H_r = 2.52 (see Figure 8). We also find that the plots of F_ω/F_ω,0 in Figure 9 look similar to those of $| {D}_{0}/{D}_{\mathrm{0,0}}{| }^{2}$ in Figure 7, respectively, except for the low frequencies. These similarities support the notion that the peak of F_ω/F_ω,0 as well as that of E_ω(h)/E_ω(h₀) results from the resonance in the photospheric layer and the effect of acoustic cutoff.

In this section, we fit the observed oscillations of the Na i line velocity presented in Figure 1. Note that the velocity data have been taken from a specific point inside the umbra of a sunspot—the point marked by P in Chae et al. (2017). This point is located in a region where numerous umbral dots exist. It is well known that umbral dots are hotter than the diffuse umbral background (Solanki 2003). Because the umbral model of Maltby et al. (1986) describes the darkest regions in large sunspot umbrae, we cannot directly use the structure of c² in this model. Instead, to set the atmospheric model of the umbral dot region, we adopt the structure of ${c}^{2}/{c}_{r}^{2}$ in the umbral model as presented in Figure 3 and multiply it by the value of c_r² that is suitable for the umbral dot region. The sound speed of the surface in the model is c_r = 6.45 km s⁻¹, which corresponds to the surface temperature of 4040 K, but we adopt c_r = 7.0 km s⁻¹, which corresponds to 4740 K for the umbral dot region, which is consistent with the reported observation that umbral dots are on average 700 K hotter than the umbral minimum (Tritschler & Schmidt 2002). In addition, we adopt the surface pressure p_r = 2.7 × 10⁵ dyn cm⁻² as in the umbral model. The specification of the two-layer model of the umbral dot region is thus completed by the adoption of the relationship between ${c}^{2}/{c}_{r}^{2}$ and h/H_r, and the choice of c_r and p_r. With this choice of c_r, we obtain ${c}_{1}=\sqrt{0.8}{c}_{r}=6.26$ km s⁻¹ and ω_ac = 0.038 rad s⁻¹. This angular frequency corresponds to the normal cutoff frequency of 5.97 mHz or to the cutoff period of 2.79 minutes.

The formation height of the Na i line should be known for the observed velocity oscillations to be compared with the theory. Chae et al. (2017) estimated the formation height of the line at 680 km by analyzing the phase differences among the different spectral lines by making use of the progressive wave solution in the piecewise isothermal atmosphere. We now understand that the waves may not be fully progressive in the low atmosphere because of the reflection near the temperature minimum, and hence the above estimate may deviate from the real formation height. Nevertheless, we adopt z_D = 680 km as an estimate of the formation height of the Na i line in the umbra because we have no other reference as yet. In this regard, this choice should be considered as a working hypothesis in the present work.

At a given height z_D in the two-layer atmosphere described in Section 4.1, one can calculate velocity v as a function of time t using the equations given in Section 3. The free parameters of the model to be specified are τ, ξ₀, t_p, t₀, and t_wv. The calculated velocity time series can be compared with the observation of the Na i line presented in Figure 1. We have determined the optimal values of τ, ξ₀, and t_p by fitting the observed global wavelet power spectrum with the model, and that of t₀ by comparing the velocity variation between the observation and the model. The global wavelet power spectra are the wavelet power spectra averaged over the time span from t₀ to t₀ + t_wv. We set t_wv to 20 minutes, which is used to obtain the observed global wavelet power shown in Figure 1.

Figure 10 shows the result of the model fitting of the observed global wavelet power spectrum. The fitting looks quite satisfactory. The determined parameters are τ = 2.77, ξ₀ = 7.4 km, and t_p = 66 s. The value of τ is used to determine the parameters at the lower boundary, namely c₀, z₀, h₀, and p₀. The values of the parameters are listed in Table 1.

Zoom In Zoom Out Reset image size

Table 1.  Parameters of the Theoretical Model Used to Interpret the Observations

Category	Parameter	Value	Unit
Input	cr	7.0	km s−1
	pr	105.43	dyn cm−2
	${c}_{1}^{2}/{c}_{r}^{2}$	0.8
	h1/Hr	0.2
	K0	−1.44
	K1	0.04
	zD	680	km
	twv	20.0	min
Fitting	τ	2.77
	ξ0	7.4	km
	tp	66	s
	t0	9.0	min
Calculated	Hr	213	km
	h0	−167	km
	h1	42	km
	hD	778	km
	z0	−276	km
	z1	40	km
	c0	10.4	km s−1
	c1	6.2	km s−1
	cD	6.8	km s−1
	κ0	−0.65
	κ1	0.05
	p0	106.11	dyn cm−2
	p1	105.26	dyn cm−2
	pD	102.26	dyn cm−2
	${\nu }_{1,\mathrm{ac}}$	5.97	mHz
	vrms	0.53	km s−1
	F	0.45	106 erg cm−2 s−1

Download table as:  ASCIITypeset image

The figure also presents the scaled power spectrum of velocity of the driving motion we used as input. The scaling was made by multiplying the square root of the density ratio between z₀ and z_D. If the medium were to allow only progressive waves with height-invariant group velocity, then the power spectrum of the velocity at z_D would be the same as the scaled input spectrum, which in fact is not the case, as can be seen from the figure. This discrepancy between the two power spectra is a clue to understanding how the medium affects the waves. First, we find that the velocity power at z_D is much stronger than the input power spectrum at frequencies ∼ω_ac/2π. This results from the resonance in the photosphere. We also find that the velocity power at z_D is significantly suppressed at low frequencies <ω_ac/2π, which is due to the acoustic cutoff. This implies that the velocity power at z_D is insensitive to the input power spectrum at frequencies ≪ω_ac/2π. Finally, the two power spectra are similar to each other in that they both decrease with frequency at high frequencies.

Figure 11 illustrates the result of such a model calculation compared with the observed velocity profile with the choice of t₀ in Table 1. We find that the model is very similar to the observation in period, phase, and amplitude modulation. This successful match supports the notion that each wave packet that is characterized by the amplitude modulation of 10 to 20 minutes results from a single event of wave excitation (Chae et al. 2017).

Zoom In Zoom Out Reset image size

Figure 12 presents the plot of the flux spectrum ${F}_{\nu }\equiv 2\pi {F}_{\omega }$ with F_ω calculated using Equation (62). The figure shows that the flux spectrum is sharply peaked at a frequency very close to the frequency ω₁/2π, and has a more enhanced wing in the higher frequency side than in the lower frequency side. The total flux integrated over frequency is found to be F = 0.45 × 10⁶ erg cm⁻² s⁻¹.

Zoom In Zoom Out Reset image size

We have also examined the case of z_D = 500 km instead of 680 km. We found that the fitting of the observed global wavelet power spectrum is as good as in the above case, with the new values τ = 2.82, ξ₀ = 8.7 km, t_p = 56 s, and t₀ = 9.1 min, which are slightly different from the above values. The most noticeable change occurs in the values p_D = 10^3.06 and F = 2.6 × 10⁶ erg cm⁻² s⁻¹. Thus we conclude that the uncertainty of z_D does not affect the general validity of our model, but is critical in the determination of F.

In this paper, we have presented the photospheric resonator model of linear acoustic waves in a nonthermal atmosphere, with the goal of modeling three-minute oscillations in the chromosphere of sunspot umbrae. Specifically, we developed a new solution that describes the propagation of acoustic waves in a model atmosphere composed of two layers: one temperature-decreasing layer, and one temperature-increasing layer. Although a realistic description may demand a more detailed multilayer modeling, this simplest two-layer model already reveals a number of observationally significant features.

Our results indicate that when a wave driver is located just below the surface, the photosphere acts as a leaky resonator, bounded by the driver at the bottom and by the temperature minimum region at the top. The temperature minimum region becomes the interface between the two layers, which partially transmits and partially reflects the waves incident from below. The three-minute oscillations observed in the Na D₂ line can be successfully explained as the progressive waves emitted by the photospheric resonator that is driven by a subphotospheric pulse of vertical displacement.

The most important outcome of our study is how the peak frequency of the velocity power spectrum ν_p is to be interpreted. We have shown that it is close to the Lamb frequency ω₁/2π of the temperature minimum region. In the specific observation we analyzed, we obtained

$$\begin{eqnarray}&&{\nu }_{p}=k\displaystyle \frac{{\omega }_{1}}{2\pi }\end{eqnarray} \tag{ 67 }$$

with a factor k that is close to unity. In our experiment, k = 1.05 is a good choice. Assuming that this relation holds generally, we can determine from the observed value of ν_p the temperature of the temperature minimum region T₁ using the formula

$$\begin{eqnarray}&&\begin{array}{l}{T}_{1}=3680\,{\rm{K}}\left[\displaystyle \frac{\gamma }{1.67}\right]\left[\displaystyle \frac{\mu }{1.25}\right]{\left[\displaystyle \frac{g}{{\rm{27,400}}\mathrm{cm}{{\rm{s}}}^{-2}}\right]}^{2}\\ \quad \times \,{\left[\displaystyle \frac{k}{1.05}\right]}^{2}{\left[\displaystyle \frac{{\nu }_{p}}{6.0\mathrm{mHz}}\right]}^{-2}.\end{array}\end{eqnarray} \tag{ 68 }$$

For instance, we have obtained ν_p = 6.1 mHz in an umbral region. The observation of Chae et al. (2018) produced an average peak period of 2.89 minutes or ν_p = 5.8 mHz from more Na i line data. By adopting γ = 1.67, μ = 1.25, and $g={\rm{27,400}}\,\mathrm{cm}\,{{\rm{s}}}^{-2}$, we obtain T₁ = 3560 and 3940 K. Thus our study provides an independent constraint on the models of the sunspot atmosphere (e.g., Maltby et al. 1986).

We have investigated the photospheric resonator model in the low atmosphere composed of the photosphere and the low chromosphere by adopting the assumption that only progressive waves exist in the chromospheric layer. There are two reasons for this approach. First, our linear theory cannot be used to investigate the behavior of waves at high altitudes. In this upper atmosphere the waves may become highly nonlinear and develop into shock waves, so the linear theory is not applicable (e.g., Chae et al. 2017; Litvinenko & Chae 2017). Second, we wished to exclude the chromospheric resonator resulting from the reflection from the chromosphere-corona transition layer to focus on the photospheric resonator. The success of the simple two-layer model in explaining the Na D₂ line observations supports the validity of this assumption and suggests that the regressive waves originating from the reflection at the chromosphere-corona transition region may have much less power than the progressive waves. This is in line with the observations (e.g., Cho et al. 2015).

Our model is directly applicable to sunspot umbrae where the magnetic field is very strong and vertical. Slow magnetoacoustic waves in such regions are practically acoustic waves that propagate along the field lines. The driving at the bottom of the atmosphere may be physically identified with the excitation of slow waves by the fast magnetoacoustic waves incident from below, which is often called the process of mode conversion (Cally 2001). Our model may be also applicable to the internetwork regions of the quiet Sun, a practically field-free medium, where acoustic waves are observationally manifest as the three-minute oscillations if the driving process occurs near the photosphere. Our model, however, may not be directly applicable to network regions and plage regions where plasma β may not be small in the low atmosphere. In these regions, a generalization of the model may be required to include the dependence on magnetic field strength.

The photospheric resonator model is not the only model that can explain the peak of the power spectrum around the acoustic cutoff frequency determined by the minimum temperature. The chromospheric resonator model as well as the propagating wave model also reproduce the peak.

The photospheric resonator model differs from the chromospheric resonator model in the number of peaks. The chromospheric resonator model predicts a few more peaks above the cutoff, closely spaced in frequency (1 mHz) in addition to the peak near the cutoff. We note that unless spatiotemporally averaged, the observed power spectrum of the chromospheric velocity oscillation is characterized by a single dominant peak around the acoustic cutoff of the temperature minimum, in agreement with our photospheric resonator model. The wavelet power spectrum given in Figure 1 is an example of spatially and temporally resolved power spectra. This power spectrum is characterized by a single peak. More wavelet power spectra of velocity determined from other spectral lines taken either from ground or from space are given in Figure 4 of Chae et al. (2017) and Figures 5 and 7 of Chae et al. (2018). These spectra indicate a successive occurrence of oscillation packets, each of which typically lasts 15 minutes. During the middle of each oscillation packet, the wavelet power spectrum usually has only one peak, which does not agree with the chromospheric resonator model. An exception occurs when the nonlinearity highly develops to be manifest as steepened oscillation profiles. In this case, the second harmonic and higher-order harmonics occur simultaneously with the fundamental peak.

The peak frequency is found to change slightly from packet to packet, so that the time-averaged wavelet power spectra (which are often called global wavelet power spectra and are equivalent to the Fourier power spectra) should be characterized by multiple peaks that are closely spaced in frequency. Because the temperature structure may differ from region to region, the power spectra that are further averaged over space may be found to have a larger number of peaks. This characteristic of the dependence on the spatiotemporal average was clearly illustrated by the observations of Christopoulou et al. (2003) as well.

The photospheric resonator model may explain the frequency variation from oscillation packet to oscillation packet. As clearly illustrated in Figure 11, the photospheric resonator model can successfully reproduce each oscillation packet. The frequency variation can be attributed either to the temporal change of temperature structure with a timescale of several tens of minutes, or to the transverse motion of fine structures of different temperature across the observing point. The first possibility has not been investigated yet, as far as we know, and the other possibility is quite likely because it is observationally well known that an umbra includes a number of fine inhomogeneous structures (such as umbral dots and umbral cores) that often move at significant speeds (e.g., Solanki 2003). Chae et al. (2018), for instance, reported a standard deviation of 0.3 minutes (0.6 mHz) and a mean of 2.89 minutes (5.8 mHz) in the peak period of the oscillation packet observed through Na i D2. With Equation (68), this variation of the peak frequency can be interpreted as the variation of the minimum temperature, with the mean being 3940 K and the standard deviation being 820 K. A temperature variation of this amount can be easily explained by the variation of the fine features.

The photospheric resonator model and the propagating wave model have subtle differences. In the chromosphere above the temperature minimum, there is no fundamental difference between the two models. In both, waves propagate without much reflection in the chromosphere, and their oscillation spectrum has a single peak. The two models differ from each other mainly in the photosphere, where the photospheric resonator is located. In the propagating wave model that does not take reflection at the temperature minimum into account, the wave power should be either constant or decrease with height. In contrast, in the photospheric resonator model, the wave power at frequencies around the cutoff is expected to increase with height because of the reflection at the temperature minimum (see Figure 8). In other words, the wave power can be accumulated and amplified in the photospheric layer. We can observationally distinguish between the propagating wave model and the photospheric resonator model when the waver power at frequencies around the cutoff is determined as a function of height.