© ESO, 2017

1. Introduction

Fleck & Schmitz (1991) were the first to point out that three-minute chromospheric oscillations naturally arise as a response of the solar atmosphere to a general velocity disturbance at its base. Kalkofen et al. (1994) presented solutions for both impulsive and continuous driving at the photospheric boundary, which illustrated how a velocity disturbance at the base of the chromosphere generates free atmospheric oscillations at the acoustic cutoff frequency. An important feature of the solutions is that the growth rate with height at the driving frequency, ω, is less than the growth rate of the oscillations excited at the acoustic cutoff frequency, ω₀, as long as ω<ω₀. Kalkofen et al. (1994) termed the former solutions evanescent because they grow with height more slowly than the free atmospheric oscillations. Consequently, the atmospheric response to a low-frequency driving is dominated at sufficiently large heights by an oscillation at the cutoff frequency, ω₀, corresponding to wave periods of about three minutes. (Non-evanescent waves at both frequencies are present if ω>ω₀.)

Independent studies of chromospheric oscillations indicated the presence of a strong three-minute peak in the data (e.g., O’Shea et al. 2002; Jefferies et al. 2006; Centeno et al. 2009). The excitation of the cutoff frequency mode provides a generally accepted explanation of the three-minute oscillations (e.g., Felipe et al. 2010; Chae & Goode 2015), although the oscillations were also modeled as standing acoustic waves trapped in a chromospheric cavity (e.g., Cally & Bogdan 1993), and as waves escaping a leaky resonator and propagating through a boundary (e.g., Zhugzhda 2008).

Theoretical studies investigated the atmospheric response to various types of continuous and pulsed driving at the lower boundary (e.g., Sutmann et al. 1998; Petukhov & Petukhov 2001). Excitation by a single pulse or by a periodic, low-frequency driver generally leads to a chromospheric response at the cutoff frequency, which decays with time. The observed persistent three-minute oscillations are probably excited by a more complicated time-dependent driving, for instance by a sequence of random pulses (Sutmann et al. 1998). Recently, Chae & Goode (2015) analyzed vertical propagation of a spectrum of acoustic waves with frequencies close to the cutoff frequency and argued that the power of the observed chromospheric oscillations at frequencies ω ≥ ω₀ is consistent with the acoustic energy flux in the photosphere, produced by a series of impulsive disturbances.

Traditionally, analytical studies of acoustic wave propagation in a gravitationally stratified atmosphere used linearized equations to describe the evolution of small perturbations. Here we use methods of nonlinear acoustics (e.g., Naugolnykh & Ostrovsky Naugolnykh & Ostrovsky 1998) to derive an analytical solution for dispersive nonlinear waves in the solar chromosphere.

The theoretical study is motivated by new opportunities for observing the three-minute oscillations, provided by instruments such as the HELioseismological Large Regions Interferometric DEvice (HELLRIDE) operating on the Vacuum Tower Telescope located at Tenerife (Wiśniewska et al. 2016) and the Fast Imaging Solar Spectrograph (FISS) installed at the 1.6-m New Solar Telescope (NST) of the Big Bear Observatory (Kwak et al. 2016). The theoretical results may also contribute to the analysis of the data from the Interferometric BIdimensional Spectropolarimeter (IBIS) installed at the Dunn Solar Telescope (e.g., Stangalini et al. 2011), the CRisp Imaging SpectroPolarimeter (CRISP) installed at the 1-m Swedish Solar Telescope (e.g., Stangalini et al. 2015), and the Interface Region Imaging Spectrograph (IRIS; e.g., Tian et al. 2014). These instruments can acquire high-cadence data in the upper solar chromosphere where, as we argue below, the nonlinear effects lead to observable consequences. The new analytical results may also help to interpret and guide numerical studies (e.g., Kalkofen et al. 2010; Fawzy & Musielak 2012).

2. Basic equations

The vertical propagation of plane adiabatic acoustic waves in the presence of gravity is governed by the momentum equation, $\mathit{\rho }\frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathit{\rho }\left(\frac{\mathit{\partial v}}{\mathit{\partial t}}\mathrm{+}\mathit{v}\frac{\mathit{\partial v}}{\mathit{\partial z}}\right)\mathrm{=}\mathrm{-}\frac{\mathit{\partial p}}{\mathit{\partial z}}\mathrm{-}\mathit{\rho g,}$(1)the continuity equation, $\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{t}}\mathrm{=}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{v}}\frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{t}}\mathrm{=}\mathrm{-}\mathit{\rho }\frac{\mathit{\partial v}}{\mathit{\partial z}}\mathit{,}$(2)and the entropy equation, $\frac{\mathrm{d}\mathit{p}}{\mathrm{d}\mathit{t}}\mathrm{=}{\mathit{c}}^{\mathrm{2}}\frac{\mathrm{d}\mathit{\rho }}{\mathrm{d}\mathit{t}}\mathrm{·}$(3)Here t is the time, z is the vertical coordinate, v(z,t) is the velocity along the z axis, ρ = ρ(v) and p = p(v) are the fluid density and pressure in the wave, $\mathit{c}\mathrm{=}\sqrt{\mathit{\gamma p}\mathit{/}\mathit{\rho }}$ is the local sound speed, and γ is the ratio of specific heats. The vector gravitational acceleration g = (0,0,−g) is directed downward.

We consider nonlinear waves in an isothermal atmosphere. The background density ρ₀(z) and pressure p₀(z) are related by the hydrostatic equation, $\frac{\mathrm{d}{\mathit{p}}_{\mathrm{0}}}{\mathrm{d}\mathit{z}}\mathrm{+}{\mathit{\rho }}_{\mathrm{0}}\mathit{g}\mathrm{=}\mathrm{0}\mathit{,}$(4)and p₀/ρ₀ = const. Consequently, $\frac{\mathrm{d}{\mathit{\rho }}_{\mathrm{0}}}{\mathrm{d}\mathit{z}}\mathrm{=}\mathrm{-}\frac{\mathrm{1}}{\mathit{H}}{\mathit{\rho }}_{\mathrm{0}}\mathit{,}$(5)where the gravitational scale height $\mathit{H}\mathrm{=}\frac{{\mathit{p}}_{\mathrm{0}}}{{\mathit{\rho }}_{\mathrm{0}}\mathit{g}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\gamma g}}\mathrm{=}\mathrm{const}\mathrm{.}$(6)and ${\mathit{c}}_{\mathrm{0}}\mathrm{=}\sqrt{\mathit{\gamma }{\mathit{p}}_{\mathrm{0}}\mathit{/}{\mathit{\rho }}_{\mathrm{0}}}$.

In the remainder of the paper, we develop an analytical description of nonlinear acoustic waves, which should apply to the observed three-minute oscillations, caused by slow magnetoacoustic waves propagating along the chromospheric magnetic field (O’Shea et al. 2002). We focus on the propagation of finite-amplitude acoustic waves and neglect the effects caused by a temperature gradient (Routh & Musielak 2014) and a nonuniform background magnetic field (Afanasyev & Nakariakov 2015). Our approach is justified as long as the background temperature gradient and magnetic field are weak enough. Quantitatively, our model remains accurate as long as its characteristic length scale, that is, the height of wave breaking and shock formation, given by Eq. (22), is small in comparison with the characteristic lengths of a more complete solution incorporating the temperature gradient and magnetic field effects. Analysis of the linearized equations suggest that those effects become significant in a low-beta plasma above the transition region separating the chromosphere and corona (Routh & Musielak 2014; Afanasyev & Nakariakov 2015).

As an illustration of distinct parameter regimes, consider the plasma density n ≃ 10⁹ cm${-3}^{}$ and temperature T ≃ 5 × 10⁵ K, corresponding to the coronal values at height h ≃ 2500 km above the photosphere (e.g., level 1 of model C in Vernazza et al. 1981). Neglecting the neutral component of the gas at that temperature, we see that the plasma beta β = 8πknT/B²< 1 if the magnetic field B exceeds a couple of Gauss. Conversely, in a much cooler and denser chromosphere, we expect our nonlinear solution – and its key prediction of an observable second harmonic at twice the driving wave frequency in an upper chromosphere – to be sufficiently reliable to justify the neglect of those additional effects.

On differentiating Eq. (1) with respect to time and simplifying the result using Eqs. (2) and (3), we get $\frac{{\mathrm{d}}^{\mathrm{2}}\mathit{v}}{\mathrm{d}{\mathit{t}}^{\mathrm{2}}}\mathrm{=}{\mathit{c}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\mathit{\gamma g}\frac{\mathit{\partial v}}{\mathit{\partial z}}\mathrm{-}\mathit{\gamma }\frac{\mathit{\partial v}}{\mathit{\partial z}}\frac{\mathrm{d}\mathit{v}}{\mathrm{d}\mathit{t}}\mathrm{·}$(7)It is straightforward to verify that Eq. (7) contains two familiar limiting cases. First, in the small perturbation theory, linearization yields $\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{=}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\mathit{\gamma g}\frac{\mathit{\partial v}}{\mathit{\partial z}}\mathrm{·}$(8)The well-known solution (in complex form) for a small-amplitude, vertically propagating wave is given by $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{~}\exp \left(\frac{\mathit{z}}{\mathrm{2}\mathit{H}}\mathrm{+}\mathrm{i}\mathit{kz}\mathrm{-}\mathrm{i}\mathit{\omega t}\right)\mathit{.}$(9)Here the frequency ω and wave number k are connected by the dispersion relation ${\mathit{\omega }}^{\mathrm{2}}\mathrm{=}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}{\mathit{k}}^{\mathrm{2}}\mathrm{+}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}\mathit{,}$(10)and the lower cutoff frequency is defined as follows: ${\mathit{\omega }}_{\mathrm{0}}\mathrm{=}\frac{{\mathit{c}}_{\mathrm{0}}}{\mathrm{2}\mathit{H}}\mathrm{·}$(11)Second, in the absence of gravity (g = 0), a simple (or Riemann) wave $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}\mathit{F}\mathrm{(}\mathit{z}\mathrm{-}\mathit{ut}\mathrm{)}$(12)is an exact solution of the nonlinear Eq. (7). Here F is an arbitrary function and $\mathit{u}\mathrm{=}\mathrm{\pm }\mathit{c}\mathrm{+}\mathit{v}\mathrm{=}\mathrm{\pm }{\mathit{c}}_{\mathrm{0}}\mathrm{+}\frac{\mathit{\gamma }\mathrm{+}\mathrm{1}}{\mathrm{2}}\mathit{v}$(13)is the speed of a point in the wave profile (e.g., Lighthill 1978). The two signs correspond to waves propagating in the positive and negative directions along the z axis.

3. Nondispersive short-wavelength waves

A simple-wave solution v(z,t) = F(z−ut) cannot exactly satisfy Eq. (7) if the gravitational acceleration g ≠ 0. To understand how acoustic waves of finite amplitude are modified by gravity, it is instructive to use a heuristic argument (Lighthill 1978) first; this argument is valid if an effective wavelength λ ~ |v/ (∂v/∂z)| is short compared with a distance over which ρ₀(z) varies significantly, that is, if λ/H ≪ 1. We break down the problem of wave propagation through a stratified medium into successive small sections of length Δz ≪ H. Within each section, the change in the background density ρ₀ is small, and the wave motion is governed by the Riemann solution. Consequently, we have v ~ F(z−ut) within each section, where the speed u is given by Eq. (13), and a proportionality constant is determined by the entry conditions. Neglecting errors that might accumulate after a number of successive sections, we approximately take v(z,t) = φ(z)F(z−ut) to describe a disturbance moving in the positive direction along the z axis. Here φ(z) is a function, defined by the density profile ρ₀(z), which varies on a scale ≫Δz. Consistency with Eq. (9) of the small perturbation theory dictates that φ(z) = exp(z/ 2H). It follows that $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{\approx }{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\mathit{F}\left[\mathit{t}\mathrm{-}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{\gamma }\mathrm{+}\mathrm{1}}{\mathrm{2}}\frac{\mathit{vz}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}\right]$(14)for an upward-propagating wave, where we redefined the arbitrary function F and expanded c₀/u, assuming v/c₀ to be sufficiently small.

To obtain a more accurate solution for finite values of the parameter λ/H, we use the results from nonlinear acoustics (e.g., Naugolnykh & Ostrovsky 1998). The method of multiple scales yields a first-order uniform expansion for finite-amplitude plane waves propagating in an inhomogeneous medium. The resulting expression, given by Eq. (32) in Nayfeh (1975), describes nonlinear adiabatic waves in an inhomogeneous medium in a hydrostatic, but not necessarily isothermal, equilibrium for which ρ₀(z) and c₀(z) are assumed to be known. On introducing the new variables $\begin{array}{ccc}& & \mathit{\tau }\mathrm{=}\mathit{t}\mathrm{-}\mathrm{\int }\frac{\mathrm{d}\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\mathit{,}\\ & & \mathit{\xi }\mathrm{=}\mathrm{\int }\frac{\mathrm{d}\mathit{z}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\sqrt{{\mathit{c}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}}\mathit{,}\\ & & \mathit{w}\mathrm{=}\sqrt{{\mathit{c}}_{\mathrm{0}}{\mathit{\rho }}_{\mathrm{0}}}\mathit{v,}\end{array}$the solution is expressed in our notation as follows: $\mathit{w}\mathrm{\left(}\mathit{\xi ,\tau }\mathrm{\right)}\mathrm{\approx }\mathit{F}\left[\mathit{\tau }\mathrm{+}\frac{\mathit{\gamma }\mathrm{+}\mathrm{1}}{\mathrm{2}}\mathit{\xi w}\right]\mathit{.}$(18)In the particular case of an isothermal background, we have ρ₀(z) = ρ₀(0)exp(−z/H), c₀ = const and $\mathit{H}\mathrm{=}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\mathit{/}\mathit{\gamma g}\mathrm{=}\mathrm{const}$. It is then straightforward to rewrite Eq. (18) in terms of the original variables. We obtain $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{\approx }{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\mathit{F}\left[\mathit{t}\mathrm{-}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\mathrm{+}\mathrm{(}\mathit{\gamma }\mathrm{+}\mathrm{1}\mathrm{)}\frac{\mathit{vH}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}\mathrm{1}\mathrm{-}{\mathrm{e}}^{\mathrm{-}\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\mathrm{)}\right]\mathit{,}$(19)which clearly reduces to Eq. (14) in the limit z/H ≪ 1. Below we model an evolving wave profile by assuming a simple harmonic perturbation at the boundary: $\mathit{v}\mathrm{\left(}\mathrm{0}\mathit{,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{0}}\sin \mathrm{\left(}\mathit{\omega t}\mathrm{\right)}\mathit{,}$(20)where ω>ω₀ for a nonevanescent propagating wave. The boundary condition yields F(t) = v₀sin(ωt).

It is worth stressing that, although the analytical multiple-scale solution formally requires the wavelength λ to be a small parameter, in practice the solution remains qualitatively correct even for λ/H > 1. This advantageous feature is common to uniform expansions derived by the method of multiple scales. The solution, however, neglects the dispersive nature of the acoustic waves because it implies that ω ≈ c₀k and thus ω ≫ ω₀ in a small-amplitude limit. Therefore, the multiple-scale solution can describe only a shortwave regime λ ≪ 2πc₀/ω₀ = 4πH, corresponding to the wave frequencies significantly exceeding the cutoff frequency ω₀. In the solar chromosphere, H ≈ 100 km, and so we expect the analytical solution to provide a reasonably accurate description for wavelengths up to a few hundred kilometres.

We can use the approximate wave solution for v(z,t) to illustrate how the gravitational field influences the formation of a shock. Wave breaking occurs at a height z_wb, defined by $\mathrm{0}\mathrm{=}\frac{\mathit{\partial \tau }}{\mathit{\partial u}}\mathrm{=}\frac{\mathit{\partial \xi }}{\mathit{\partial u}}\mathrm{=}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\tau }}{\mathit{\partial }{\mathit{u}}^{\mathrm{2}}}\mathrm{=}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\xi }}{\mathit{\partial }{\mathit{u}}^{\mathrm{2}}}$(21)(Ostrovskii 1963). On using Eqs. (19) and (20), we get ${\mathit{z}}_{\mathrm{wb}}\mathrm{=}\mathrm{2}\mathit{H}\ln \left(\mathrm{1}\mathrm{+}\frac{\mathit{\gamma }}{\mathit{\gamma }\mathrm{+}\mathrm{1}}\frac{\mathit{g}}{\mathit{\omega }{\mathit{v}}_{\mathrm{0}}}\right)$(22)or alternatively ${\mathit{z}}_{\mathrm{wb}}\mathrm{=}\mathrm{2}\mathit{H}\ln \left(\mathrm{1}\mathrm{+}\frac{{\mathit{z}}_{\mathrm{wb}\mathit{,}\mathrm{0}}}{\mathrm{2}\mathit{H}}\right)\mathit{,}$(23)where ${\mathit{z}}_{\mathrm{wb}\mathit{,}\mathrm{0}}\mathrm{=}\frac{\mathrm{2}}{\mathit{\gamma }\mathrm{+}\mathrm{1}}\frac{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{\omega }{\mathit{v}}_{\mathrm{0}}}$(24)is the height of wave breaking in the absence of gravity (H = ∞). Clearly gravity acts to lower the height of shock formation. Physically, gravity modifies the background density profile, causing an increase of the amplitude of an upward-propagating wave, which in turn leads to a stronger deformation of an initial wave profile. Notably, z_wb → 0 as H → 0, confirming that the analytical multiple-scale solution for v(z,t) leads to a qualitatively correct description of nonlinear waves even if λ/H> 1.

Nonlinear oscillating systems generally exhibit oscillations with combination frequencies that are superposed on the normal oscillations of the system (Landau & Lifshitz 1969). Here we consider the basic nonlinear effect of a second harmonic excitation by nonlinear acoustic waves in the solar atmosphere, which may produce observable features at double frequency 2ω. Kalkofen et al. (1994) appear to have identified this effect in their numerical simulations of chromospheric oscillations.

The second harmonic generation is demonstrated most easily for moderately nonlinear, short-wavelength waves, in which case we can replace v(z,t) on the right-hand side of Eq. (19) by a linear solution from Eq. (9) in the limit ω ≫ ω₀, $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{\approx }{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}{\mathit{v}}_{\mathrm{0}}\sin [\mathit{\omega t}\mathrm{-}\frac{\mathit{\omega z}}{{\mathit{c}}_{\mathrm{0}}}\begin{array}{c}\end{array}\\ \begin{array}{c}\end{array} \mathrm{+}\mathrm{(}\mathit{\gamma }\mathrm{+}\mathrm{1}\mathrm{)}\frac{\mathit{\omega H}{\mathit{v}}_{\mathrm{0}}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\sin \left(\mathit{\omega t}\mathrm{-}\frac{\mathit{\omega z}}{{\mathit{c}}_{\mathrm{0}}}\right)]\mathrm{·}\end{array}& & \end{array}$(25)Now v(z,t) can be written as its Fourier series $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{\approx }{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}{\mathit{v}}_{\mathrm{0}}\sum _{\mathit{n}\mathrm{=}\mathrm{1}}^{\mathrm{\infty }}{\mathit{b}}_{\mathit{n}}\sin \mathit{n}\left(\mathit{\omega t}\mathrm{-}\frac{\mathit{\omega z}}{{\mathit{c}}_{\mathrm{0}}}\right)\mathit{,}$(26)where $\begin{array}{ccc}& & {\mathit{b}}_{\mathrm{1}}\mathrm{=}\frac{\mathrm{2}}{\mathit{a}}{\mathit{J}}_{\mathrm{1}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{-}\mathrm{2}{\mathit{J}}_{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\mathrm{1}\mathrm{-}\frac{\mathrm{3}}{\mathrm{8}}{\mathit{a}}^{\mathrm{2}}\mathrm{+}\mathit{...}\mathit{,}\\ & & {\mathit{b}}_{\mathrm{2}}\mathrm{=}\frac{\mathrm{4}}{\mathit{a}}{\mathit{J}}_{\mathrm{2}}\mathrm{\left(}\mathit{a}\mathrm{\right)}\mathrm{=}\frac{\mathrm{1}}{\mathrm{2}}\mathit{a}\mathrm{-}\frac{\mathrm{1}}{\mathrm{24}}{\mathit{a}}^{\mathrm{3}}\mathrm{+}\mathit{...}\mathit{,}\end{array}$and we used the Taylor expansions of the Bessel functions J₁ and J₂ in terms of a dimensionless velocity amplitude $\mathit{a}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{+}\mathrm{1}\mathrm{)}\frac{\mathit{\omega H}{\mathit{v}}_{\mathrm{0}}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{(}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\mathrm{-}{\mathrm{1}}^{\mathrm{)}}\mathit{.}$(29)The wave-breaking condition, specified by Eq. (22), is equivalent to a = 1.

For a small nonlinearity parameter a, the terms of order a² and higher can be neglected, and so the first two terms in the Fourier series provide an accurate approximation of the complete solution, $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{\approx }{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}{\mathit{v}}_{\mathrm{0}}\left[\sin \left(\mathit{\omega t}\mathrm{-}\frac{\mathit{\omega z}}{{\mathit{c}}_{\mathrm{0}}}\right)\mathrm{+}\frac{\mathit{a}}{\mathrm{2}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}\frac{\mathit{\omega z}}{{\mathit{c}}_{\mathrm{0}}}\right)\right]\mathit{,}$(30)which is just the sum of the linear solution and a second-harmonic correction. This expression also follows directly from Eq. (25) on replacing its right-hand side by a linear Taylor approximation with respect to the dimensionless amplitude a.

4. Dispersive long-wavelength waves

As noted above, a nondispersive solution for acoustic wave propagation in the presence of gravity can only be justified in a limiting case of short wavelengths. In practice, the dispersive nature of the waves is essential in a long-wave regime λ > 2πc₀/ω₀, which corresponds to 0 < ω/ω₀−1 ≪ 1. As an illustration, suppose that the gravitational scale height in a lower chromosphere is H = 100 km, g = 0.27 km s^-2, and so c₀ ≈ 6.7 km s^-1. On calculating the cutoff frequency ω₀ from Eq. (11) and using ω = 2π/T with T = 180 s, we get ω/ω₀ ≈ 1.1. Clearly Eq. (30) is not applicable in this case.

In order to derive a dispersive-wave solution, we begin by noting that moderately nonlinear waves can be described analytically if v/c₀ is sufficiently small. We obtain a second-order approximation by retaining only linear and quadratic terms with respect to v in Eq. (7) as follows: $\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{+}\mathit{\gamma g}\frac{\mathit{\partial v}}{\mathit{\partial z}}\mathrm{=}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\left(\frac{{\mathit{p}}_{\mathrm{1}}}{{\mathit{p}}_{\mathrm{0}}}\mathrm{-}\frac{{\mathit{\rho }}_{\mathrm{1}}}{{\mathit{\rho }}_{\mathrm{0}}}\right)\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{-}\mathrm{2}\mathit{v}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{v}}{\mathit{\partial z\partial t}}\mathrm{-}\mathrm{(}\mathrm{1}\mathrm{+}\mathit{\gamma }\mathrm{)}\frac{\mathit{\partial v}}{\mathit{\partial z}}\frac{\mathit{\partial v}}{\mathit{\partial t}}\mathit{,}$(31)where the quadratic terms are collected on the right-hand side of the equation (see also Petukhov & Petukhov 2001 and references therein). Because ρ₁ and p₁ are multiplied by a derivative of v, in this second-order approximation the density and pressure perturbations should be calculated from the linearized Eqs. (2) and (3), i.e., $\begin{array}{ccc}& & \\ & & \end{array}$As a concrete example, consider the chromospheric response to a source of monochromatic acoustic waves with frequency ω > ω₀ at the lower boundary z = 0. The boundary condition is given by Eq. (20). On neglecting the right-hand side of Eq. (31), we obtain a linear equation for small-amplitude waves, which yields the familiar linear solution ${\mathit{v}}_{\mathrm{lin}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{0}}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\sin \left(\mathit{\omega t}\mathrm{-}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\right)\mathit{,}$(34)or in complex form, ${\mathit{v}}_{\mathrm{lin}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{0}}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\mathrm{Im}\left[\exp \left(\mathrm{i}\mathit{\omega t}\mathrm{-}\mathrm{i}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\right)\right]\mathrm{·}$(35)In the next iteration, we solve Eq. (31) to obtain a correction δv to v_lin(z,t), such that v = v_lin + δv is correct to a second order in v₀. Equations (32) and (33) yield the first-order density and pressure perturbations, $\frac{{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathrm{=}\mathrm{Im}\left[\frac{{\mathit{\omega }}_{\mathrm{0}}\mathrm{+}\mathrm{i}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}}{\mathit{i\omega }{\mathit{c}}_{\mathrm{0}}}{\mathit{v}}_{\mathrm{lin}}\right]\mathit{,}$(36)$\frac{{\mathit{p}}_{\mathrm{1}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}}{{\mathit{p}}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathrm{=}\mathit{\gamma }\frac{{\mathit{\rho }}_{\mathrm{1}}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}}{{\mathit{\rho }}_{\mathrm{0}}\mathrm{\left(}\mathit{z}\mathrm{\right)}}\mathrm{+}\mathrm{Im}\left[\frac{\mathrm{2}\mathrm{(}\mathrm{1}\mathrm{-}\mathit{\gamma }\mathrm{)}{\mathit{\omega }}_{\mathrm{0}}}{\mathit{i\omega }{\mathit{c}}_{\mathrm{0}}}{\mathit{v}}_{\mathrm{lin}}\right]\mathit{,}$(37)where the complex form is used for compactness. On substituting v_lin, ρ₁, and p₁ into the right-hand side of Eq. (31), we obtain a linear equation for δv. The general solution is a sum of solutions of the corresponding homogeneous and inhomogeneous equations. We select a solution of the former equation, which describes an upward-propagating wave, whose amplitude increases with height as exp(z/ 2H), which is the same as the growth rate of the linear wave. The solution of the latter equation describes a wave, whose amplitude increases with height as exp(z/H). Given the weak dependence of Eq. (30) on γ, we assume that dispersive solutions also depend weakly on γ and simplify the algebra by performing the analysis for isothermal perturbations. Then Eq. (31) yields $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\delta v}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\delta v}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{+}\mathit{\gamma g}\frac{\mathit{\partial \delta v}}{\mathit{\partial z}}\mathrm{=}\frac{\mathrm{2}\mathit{\omega }{\mathit{v}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{0}}}\mathrm{Im}[\left(\mathrm{i}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}\right)\\ \begin{array}{c}\end{array}\mathrm{\times }\exp \left(\mathrm{i}\mathrm{2}\mathit{\omega t}\mathrm{-}\mathrm{i}\mathrm{2}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{z}}{\mathit{H}}\right)]\mathit{,}\end{array}& & \end{array}$(38)and the solution for δv is given by $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{\delta v}\mathrm{=}\frac{\mathit{\omega H}{\mathit{v}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}\mathrm{Im}[\exp \left(\mathrm{i}\mathrm{2}\mathit{\omega t}\mathrm{-}\mathrm{i}\mathrm{2}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{z}}{\mathit{H}}\right)\\ \begin{array}{c}\end{array}\mathrm{-}\exp \left(\mathrm{i}\mathrm{2}\mathit{\omega t}\mathrm{-}\mathrm{i}\sqrt{\mathrm{4}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\mathrm{+}\frac{\mathit{z}}{\mathrm{2}\mathit{H}}\right)]\mathit{,}\end{array}& & \end{array}$(39)where integration constants in the solution of the homogeneous equation are determined from the boundary condition $\mathit{\delta v}\mathrm{\left(}\mathrm{0}\mathit{,t}\mathrm{\right)}\mathrm{=}\mathrm{0}\mathit{,}$(40)which describes the absence of a second harmonic at z = 0. Of course the factor exp(z/H) in the solution indicates the breakdown of the formal expansion in the small parameter v₀/c₀ at large heights.

Next we use a simple matching argument to correct the nondispersive nonlinear solution given by Eq. (30). We modify Eq. (30), so that it reduces to the dispersive solution given by Eqs. (34) and (39) when ω ≈ ω₀ and γ = 1. The resulting approximate solution for an upward-propagating wave incorporates both the nonlinear and dispersive effects $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{lin}}\mathrm{+}\mathit{\delta v,}$(41)where v_lin is defined by Eq. (34), and the second-order correction δv is given by $\begin{array}{ccc}\hspace{0.17em}\begin{array}{c}\begin{array}{c}\end{array}\mathit{\delta v}\mathrm{=}\mathrm{(}\mathit{\gamma }\mathrm{+}\mathrm{1}\mathrm{)}\frac{\mathit{\omega H}{\mathit{v}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}}[{\mathrm{e}}^{\mathit{z}\mathit{/}\mathit{H}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}\sqrt{{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{{\mathit{c}}_{\mathrm{0}}}\right)\\ \begin{array}{c}\end{array}\mathrm{-}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\sin \mathrm{2}\left(\mathit{\omega t}\mathrm{-}\sqrt{\mathrm{4}{\mathit{\omega }}^{\mathrm{2}}\mathrm{-}{\mathit{\omega }}_{\mathrm{0}}^{\mathrm{2}}}\frac{\mathit{z}}{\mathrm{2}{\mathit{c}}_{\mathrm{0}}}\right)]\mathrm{·}\end{array}& & \end{array}$(42)The solution generalizes Eqs. (30) and (39) and therefore describes the second harmonic excitation in both the short-wavelength and long-wavelength limits.

In the limit ω → ω₀, the matching procedure that yields Eq. (42) leads to an error on the order of (γ−1)/2. Consider a second-order approximation in the limiting case ω = ω₀ for an arbitrary γ. Then Eq. (31) yields $\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\delta v}}{\mathit{\partial }{\mathit{t}}^{\mathrm{2}}}\mathrm{-}{\mathit{c}}_{\mathrm{0}}^{\mathrm{2}}\frac{{\mathit{\partial }}^{\mathrm{2}}\mathit{\delta v}}{\mathit{\partial }{\mathit{z}}^{\mathrm{2}}}\mathrm{+}\mathit{\gamma g}\frac{\mathit{\partial \delta v}}{\mathit{\partial z}}\mathrm{=}\mathrm{-}\frac{{\mathit{\omega }}_{\mathrm{0}}{\mathit{v}}_{\mathrm{0}}^{\mathrm{2}}}{\mathit{H}}\sin \mathrm{(}\mathrm{2}{\mathit{\omega }}_{\mathrm{0}}{\mathit{t}}^{\mathrm{)}}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathit{H}}\mathit{,}$(43)and the solution for δv is given by $\mathit{\delta v}\mathrm{=}\frac{{\mathit{v}}_{\mathrm{0}}^{\mathrm{2}}}{\mathrm{2}{\mathit{c}}_{\mathrm{0}}}\left[{\mathrm{e}}^{\mathit{z}\mathit{/}\mathit{H}}\sin \mathrm{(}\mathrm{2}{\mathit{\omega }}_{\mathrm{0}}{\mathit{t}}^{\mathrm{)}}\mathrm{-}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\sin \left(\mathrm{2}{\mathit{\omega }}_{\mathrm{0}}\mathit{t}\mathrm{-}\frac{\sqrt{\mathrm{3}}}{\mathrm{2}}\frac{\mathit{z}}{\mathit{H}}\right)\right]\mathit{,}$(44)which depends on γ only through the dependence on c₀. A more detailed matching of the solutions with ω ≈ ω₀ and ω ≫ ω₀ could be performed but would not lead to qualitatively different results in the second-order approximation.

Finally, although the nonlinear solutions above describe upward-propagating waves, Eq. (31) with different initial conditions also admits both downward-propagating and standing-wave solutions, which may be of physical interest. In practice, the initial and boundary conditions are hard to specify for realistic photospheric driving. For instance, the simplest assumption that the atmosphere be initially at rest and unperturbed, v(z,0) = ρ₁(z,0) = p₁(z,0) = 0, leads to a standing-wave solution. As an illustration, consider the low-frequency limit ω = ω₀. The linear standing wave is described by ${\mathit{v}}_{\mathrm{lin}}\mathrm{=}{\mathit{v}}_{\mathrm{0}}{\mathrm{e}}^{\mathit{z}\mathit{/}\mathrm{2}\mathit{H}}\sin \mathrm{(}{\mathit{\omega }}_{\mathrm{0}}{\mathit{t}}^{\mathrm{)}}\mathit{.}$(45)As previously, we find a weakly nonlinear solution of Eq. (31) by substituting Eq. (45) into the right-hand side of Eq. (31) and solving the resulting linear equation, which leads to $\mathit{v}\mathrm{\left(}\mathit{z,t}\mathrm{\right)}\mathrm{=}{\mathit{v}}_{\mathrm{lin}}\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathrm{1}}\mathrm{+}\mathit{\delta }{\mathit{v}}_{\mathrm{2}}\mathit{.}$(46)Here δv₁ is a second-order correction to the amplitude of the oscillation at the driving frequency ω₀, and δv₂ is the second harmonic, which is of primary interest to us. On employing the zero initial and boundary conditions, after some algebra we obtain $\begin{array}{ccc}\mathit{\delta }{\mathit{v}}_{\mathrm{1}}& \mathrm{=}& \mathrm{(}\mathit{\gamma }\mathrm{-}\mathrm{1}\mathrm{)}\frac{{\mathit{v}}_{\mathrm{0}}^{\mathrm{2}}}{{\mathit{c}}_{\mathrm{0}}}\left[\exp \left(\frac{\mathit{z}}{\mathit{H}}\right)\mathrm{-}\exp \left(\frac{\mathit{z}}{\mathrm{2}\mathit{H}}\right)\right]\sin \mathrm{\left(}{\mathit{\omega }}_{\mathrm{0}}\mathit{t}\mathrm{\right)}\mathit{,}\\ \mathit{\delta }{\mathit{v}}_{\mathrm{2}}& \mathrm{=}& \end{array}$While both terms in Eq. (48) are nonevanescent, the spatially nonoscillating term dominates at larger heights.

5. Discussion

The observed three-minute oscillations in the solar chromosphere, caused by slow magnetoacoustic waves, can be modeled under certain simplifying assumptions as vertically propagating acoustic waves. We used methods of nonlinear acoustics to describe vertical propagation of finite-amplitude acoustic waves in a gravitationally stratified atmosphere. Since dispersive effects are essential in the propagation of the low-frequency waves responsible for the chromospheric oscillations, we derived a nonlinear dispersive solution that is valid in a long-wavelength limit. We also matched the solution to a nondispersive solution that is valid in a short-wavelength limit. Therefore, some information can be obtained in an intermediate case by interpolating between the two extreme cases.

A key feature of the analytical solution is that it yields a simple description of the basic nonlinear effect of the second harmonic generation at frequency 2ω, where the wave driving frequency ω only slightly exceeds the chromospheric cutoff frequency ω₀ ≈ 0.03 s^-1. It appears that this effect was previously noted in a numerical study of chromospheric oscillations (Kalkofen et al. 1994). More general solutions incorporating the effects caused by a temperature gradient (Routh & Musielak 2014) or a nonuniform background magnetic field (Afanasyev & Nakariakov 2015) would be required to achieve a more accurate description of the waves in the transition region where the plasma beta decreases rapidly with height.

It is worth stressing that the nonlinearity of the system controls the evolution of a large-amplitude acoustic wave, and the strength of the second harmonic is quantified by the dimensionless nonlinearity parameter a, defined by Eq. (29), rather than just by a dimensionless wave amplitude v₀/c₀. Equations (22) and (29) show that wave breaking occurs at a height where a = 1. It follows from Eq. (30) that wave breaking occurs when the second-harmonic amplitude reaches a half of the wave amplitude at the driving frequency ω.

Our results suggest that the nonlinearity of acoustic waves is strong enough in the solar chromosphere to allow the experimental detection of a second harmonic signal. On setting a = 1 in Eq. (29), we arrive at a rough estimate of the height, where the strongest second harmonic signal should be expected, as follows: $\mathit{z}\mathrm{\approx }\mathrm{2}\mathit{H}\ln \mathrm{(}{\mathit{c}}_{\mathrm{0}}\mathit{/}{\mathit{v}}_{\mathrm{0}}\mathrm{)}\mathit{.}$(49)This simple analytical expression, derived in the nondispersive limit, remain essentially unaltered for the dispersive solution given by Eq. (42). If v₀ is assumed to be comparable with the photospheric convection speed, then c₀/v₀ ≈ 100, and so z ≈ 9H ≈ 900 km, implying that the second harmonic signal in chromospheric oscillations should be detectable in an upper chromosphere by an instrument such as the FISS NST at the Big Bear Observatory (Chae et al. 2013).

Work is currently underway to isolate the second harmonic signal in the FISS data (J. Chae et al., in prep.).

Acknowledgments

Y.L. is grateful to the hospitality of Seoul National University where this work began. The authors acknowledge the comments and suggestions by the referee, Marco Stangalini.

References