Asteroseismic Signature of a Large Active Region

1. Introduction

Surface magnetic activity shifts the frequencies of the global modes of acoustic oscillation during solar and stellar activity cycles (e.g., Pallé et al., 1989; García et al., 2010; Santos et al., 2016; Kiefer et al., 2017). Asteroseismology can in turn inform us about the strength and the latitude distribution of a band of magnetic activity on the stellar surface (Gizon, 2002; Chaplin et al., 2007; Gizon et al., 2016). However, the signature of a single large active region in stellar p-mode oscillation power spectra has not been discussed so far in detail. A complication inherent to this problem comes from the fact that the perturbation associated with a rotating active region is not steady in the observer's frame. Yet, this problem is relevant, since large starspots are detected in photometric variability (e.g., Mosser et al., 2009) and using Doppler or Zeeman-Doppler imaging (e.g., Strassmeier, 2009).

We build on preliminary work by Gizon (1995, 1998) who considered the effect of a single sunspot in corotation on the low-degree modes of solar oscillation. The problem of unsteady perturbations has been considered in the past in different contexts. The interaction with acoustic modes with an inclined magnetic field with respect to the stellar rotation axis was first studied by Kurtz (1982) to explain the oscillation power spectra of roAp stars (see also Kurtz, 2008). In the oblique pulsator model the effect of the magnetic field dominates over rotation, and the pulsation axis is aligned with the magnetic axis of the star. Only modes of oscillations symmetric with respect to the magnetic field axis are excited. The oblique pulsator model was extended by Dziembowski and Goode (1985) to include the first-order effects of the Coriolis and Lorentz forces, and then by Shibahashi and Takata (1993) to account for the distortion in the eigenfunctions. In parallel Dziembowski and Goode (1984) and Gough and Taylor (1984) discussed the combined influence of rotation and of an inclined magnetic field in corotation on multiplets of solar acoustic oscillations. They explicitly mentioned that each multiplet consist of (2ℓ + 1)² components in the power spectrum (this was already hinted at by Dicke, 1982).

We focus on stars with a level of activity higher than the Sun, which may have active regions with larger surface coverage, and therefore better chances for detection. Following the same approach of Goode and Thompson (1992), we investigate the linear changes induced in a nℓ-multiplet by an unsteady perturbation that mimics an active region (AR) rotating with the star. In particular we study the power spectra of the dipole and quadrupole multiplets. For the active region we consider a simple two-parameter model, where near-surface sound-speed perturbations have a given amplitude and surface coverage.

As a complement to the linear analysis we also explore the non-linear regime of the active-region perturbation by means of 3D numerical simulations, by studying the combined effect of rotation and mode mixing on the observed power spectra using the wave propagation code GLASS (see Hanasoge and Duvall, 2007; Papini et al., 2015). We note that the non-linear regime was studied in the context of strong magnetic fields in roAp stars by, e.g., Cunha and Gough (2000); Bigot and Dziembowski (2002); Saio and Gautschy (2004), and Cunha (2006).

2. Methods

2.1. Signature of an Active Region in the Oscillation Power Spectrum: Linear Theory

2.1.1. Linear Problem in the Corotating Frame

The normal modes of oscillation of a spherically symmetric non-rotating star are identified by three integer numbers: the radial order n, the angular degree ℓ, and the azimuthal order m, with |m| ≤ ℓ. In the absence of attenuation the degenerate mode frequencies, ${\omega }_{n\ell }^{(0)}$, are real and the displacement eigenvectors ${\xi }_{n\ell m}^{(0)}(\text{r})exp(-\text{i}{\omega }_{n\ell }^{(0)}t)$ solve the linearized equation of motion

where ${{L}}^{(0)}$ is a linear spatial differential operator (see e.g., Unno et al., 1989). Hereafter superscripts “(0)” denote quantities associated with the non-rotating stellar model. In spherical polar coordinates r = (r, θ, ϕ), the displacement eigenfunctions can be written as

where $Y_l^m(\theta ,\varphi )$ are spherical harmonics and ξ_{r, nℓ} and ξ_{h, nℓ} are functions of radius only that can be calculated numerically for a given stellar model (see e.g., Aerts et al., 2010).

We now consider two perturbations, the first perturbation due to rotation (e.g., Hansen et al., 1977) and the second arising from the presence of an active region that rotates with the star (see e.g., Schunker et al., 2013, for the effect of a sunspot on high-degree modes). The latter perturbation is unsteady in an inertial frame of reference. Here we aim to study how these two effects may affect the fine structure of the modes within a fixed multiplet (n, ℓ). The two perturbations taken together completely remove the (2ℓ + 1)-fold degeneracy in m.

Provided that there is only one active region, it is much more convenient to first tackle the problem in a reference frame that is corotating with the active region, where both perturbations are steady (e.g., Dziembowski and Goode, 1984; Goode and Thompson, 1992). In Figure 1 we define three frames of reference, ${R}$, ${{R}}_{\beta }$, and ${\hat{{R}}}_{\beta }$, all three with the same origin at the center of the star. Frame ${R}$ is an inertial frame of reference, with polar axis z pointing toward the observer. We denote by β the colatitude of the active region in frame ${R}$. The other two frames are both corotating with the active region at the angular velocity Ω_β about the rotation axis of the star. The polar axis of ${{R}}_{\beta }$ is directed along the stellar rotation axis and is inclined by an angle i with respect to z. The polar axis of the frame ${\hat{{R}}}_{\beta }$ is inclined by the angle β with respect to the rotation axis. In ${\hat{{R}}}_{\beta }$ the center of the active region is at the north pole. Provided that the starspot has no proper motion, the angular velocity Ω_β is equal to the surface rotational angular velocity of the star at colatitude β. We call r = (r, θ, ϕ), r_β = (r, θ_β, ϕ_β), and ${\hat{\text{r}}}_{\beta }=(r,{\hat{\theta }}_{\beta },{\hat{\varphi }}_{\beta })$ the spherical-polar coordinates associated with ${R}$, ${{R}}_{\beta }$, and ${\hat{{R}}}_{\beta }$, respectively.

We consider the effects of rotation and the active region on the acoustic oscillations as small perturbations. In the frame ${{R}}_{\beta }$ each mode is identified by the index M, −ℓ ≤ M ≤ ℓ, and we expand the displacement eigenvectors and eigenfunctions as

and

where δξ_nℓM is orthogonal to each unperturbed eigenvector ${\xi }_{n\ell m}^{(0)}$ with the same ℓ and n (e.g., Gough and Thompson, 1990), and the coefficients $A_m^M$ are (real) amplitudes. We write the wave operator as

with

where ${{L}}_{\Omega }$ accounts for the effects of rotation and ${{L}}_{\text{AR}}$ for the effects of the active region.

To first order, the linearized equation of motion reduces to

We define the inner product between two vectors ξ(r_β) and η(r_β) on the Hilbert space of displacement vectors as

where ^* denotes the complex conjugate and V is the stellar volume. The unpertubed eigenmodes are normalized such that

We take the inner product of Equation (7) with ${\xi }_{n\ell m^{\prime }}^{(0)}$ to obtain

Because ${{L}}^{(0)}$ is Hermitian symmetric and $\langle {\xi }_{n\ell m^{\prime }}^{(0)},\delta {\xi }_{n\ell M}\rangle =0$, the second term on the left-hand side of the above equation vanishes

Introducing the perturbation matrix elements

where

Equation (9) becomes

To simplify the notation we dropped the indices nℓ on δω_M. In matrix form,

where $A^M={\left[A_{-\ell }^M{A}_{-\ell +1}^M\dots A_{\ell }^M\right]}^T$ is the vector of amplitudes. To find A^M and δω_M we have to solve the above eigenvalue problem, Equation (14).

The rotation perturbation matrix O^Ω is diagonal in the frame ${{R}}_{\beta }$. The active region perturbation matrix O^AR is not diagonal in ${{R}}_{\beta }$, but it can be obtained in terms of the diagonal perturbation matrix ${\widehat{O}}^{\text{AR}}$ expressed in the frame ${\hat{{R}}}_{\beta }$,

where the matrix R^(ℓ) performs a clockwise rotation of β about the y axis that transforms the frame ${\hat{{R}}}_{\beta }$ into the frame ${{R}}_{\beta }$. More explicitly, the elements of the rotation matrix are given by

where $r_{m{m}^{\prime }}^{(\ell )}(\beta )$ is given by Messiah (1960).

2.1.2. Frequency Splittings Due to Rotation

In the corotating frame ${{R}}_{\beta }$ the rotation perturbation matrix is diagonal:

with

where Ω(r, θ) is the internal angular velocity in an inertial frame. By construction, the angular velocity Ω_β of the frame ${{R}}_{\beta }$ is Ω_β = Ω(R, β), where R is the radius of the star. The first and second terms on the right-hand side of Equation (18) describe the effect of differential rotation, where the functions K_nℓm(r, θ) are the rotational sensitivity kernels (Hansen et al., 1977) and C_nℓ are the Ledoux constants (Ledoux, 1951) that account for the effect of the Coriolis force. The last term describes the quadrupole distortion of the star due to the centrifugal forces (e.g., Saio, 1981; Aerts et al., 2010) and is proportional to the ratio of the centrifugal to the gravitational forces at the surface η = Ω²R³/(GM), where M is the mass of the star and G is the universal constant of gravity. The term Q_2ℓm accounts for the quadrupolar component of the centrifugal distorsion

where $P_{\ell }^m(x)$ are the associated Legendre functions and P₂ is the Legendre polynomial of second order. The centrifugal term is very small in the case of slow rotators like the Sun (Dziembowski and Goode, 1992), however it increases rapidly with rotation and it is not negligible anymore for stars rotating a few times faster than the Sun (e.g., Gizon and Solanki, 2004).

2.1.3. Frequency Splittings Due to the Active Region

In this section we parametrize the effects of the active region on the oscillation frequencies in the corotating frame ${\hat{{R}}}_{\beta }$, where the active region is at the north pole.

Modeling the complex influence of surface magnetic fields on acoustic oscillations is challenging (e.g., Gizon et al., 2010; Schunker et al., 2013). Here we choose to drastically simplify the physics and to focus on the geometrical aspects of the problem.

Assuming that the area of the active region covers a polar cap with $0\le {\hat{\theta }}_{\beta }\le \alpha$ (see Figure 1) and that the structure of the active region is separable in r and ${\hat{\theta }}_{\beta }$, we can parametrize the perturbation matrix as follows:

where

is a geometrical weight factor that accounts for the surface coverage of the active region. For small values of α, $G_{\ell }^{m^{″}}(\alpha )$ decreases fast as |m″| increases (Figure 2). Since the value of $G_{\ell }^{m^{″}}(\alpha )$ does not depend on the sign of m″, the eigenvalues of O^AR are degenerate in |m| (see also, e.g., Kurtz and Shibahashi, 1986).

The parameter ε_nℓ is a measure of the relative magnitude of the active region perturbation and contains all the physics. From studies of local helioseismology (e.g., Gizon et al., 2009; Moradi et al., 2010; Schunker et al., 2013), it is known that the net effect of an active region is to increase the frequencies of acoustic modes, i.e., waves propagate faster in magnetic regions and thus ε_nℓ is positive. Since the active region introduces a perturbation that is strongly localized near the surface, the value of ε_nℓ increases with radial order n. A proper calculation of the value of ε_nℓ goes beyond the scope of the present study. Instead, we parametrize the active region perturbation as an increase in sound speed near the surface. Following Papini et al. (2015), we write

where ρ₀(r) is the density of the unperturbed stellar background and Δc²(r) is the radial change in the squared sound speed. Here the integral is over the radial extent of the active region. In section 2.2 we will further specify the effective sound-speed perturbation.

Using Equation (15), the perturbation matrix elements in frame ${{R}}_{\beta }$ are

We note that for ℓ ≤ 3 and α ≲ 30°, we have $G_{\ell }^0(\alpha )\gg G_{\ell }^{m^{″}}(\alpha )$, where |m″| ≠ 0, and therefore the dominant eigenvalue is $O_{00}^{\text{AR}}={\omega }_{n\ell }^{(0)}{\epsilon }_{n\ell }{G}_{\ell }^0(\alpha )$, see Equation (20).

2.2. Numerical Setup for the Linear Problem

We now introduce the internal rotation model and the active region parameters, used to illustrate the theory. We consider a star with a rotation period of 8 days, about one third the rotation period of the Sun. This choice of rotation period ensures that the azimuthal modes in a multiplet are well separated in frequency space. For the internal rotation profile, we take

which is a scaled model of solar differential rotation as in Gizon and Solanki (2004). The centrifugal term has a significant effect. It shifts the m = 0 mode and introduces an asymmetry in the shifts for positive and negative azimuthal orders m, with a maximum frequency shift of more than 100 nHz in the case of a multiplet near 3 mHz. Therefore this term must be included when performing the analysis.

From observations of p-mode frequency changes during the solar cycle, Libbrecht and Woodard (1990) showed that the (positive) frequency shifts are almost independent of ℓ and increase with frequency, thus indicating that the effects of magnetic activity on acoustic oscillations are confined to the surface. Assuming that the perturbation covers two pressure scale heights below the photosphere and setting Δc²/c² ≃ 10% there, Equation (22) gives ε_nℓ ≃ 0.003 at frequencies near 3 mHz. The surface coverage of a stellar active region, as inferred from Doppler imaging, ranges from a percent up to 10% (Strassmeier, 2009). Here we consider two different surface coverages of either 4% or 7%, corresponding to cosα = 0.92 or 0.86 (α ≃ 23° or 30°). Finally, we consider an active region at a colatitude of either β = 20° (near the pole) or 80° (near the equator).

In the following section we focus on two multiplets with (ℓ, n) = (1, 18) and (2, 18). For each of these modes we solve the eigenvalue problem (14) by means of Jacobi's method (Press et al., 1992). For the calculation of the unperturbed eigenmodes $({\omega }_{n\ell }^{(0)},{\xi }_{n\ell m}^{(0)})$ we use the ADIPLS software package (Christensen-Dalsgaard, 2008) and solar Model S as the reference structure model (Christensen-Dalsgaard et al., 1996). The unperturbed frequencies ${\omega }_{n\ell }^{(0)}/2\pi$ of the dipole and quadrupole modes are 2695.40 and 2756.95 μHz, respectively. We note that the choice of reference solar model is unimportant for the present study.

2.3. Power Spectrum in the Observer's Frame: (2ℓ + 1)² Peaks

Given particular values for α, β, and ε_nℓ the eigenvalue problem (14) is fully specified and can be solved. In this section we use the solutions (3) and (4) to build a synthetic power spectrum in the observer's frame, in order to relate the results to observations.

We need to find an expression that connects the eigenmodes to the observed intensity fluctuations. For the sake of simplicity, we assume that the variation I(θ_β, ϕ_β, t) induced by the acoustic oscillations in the emergent photospheric intensity is proportional to the Eulerian pressure perturbation p of the acoustic wavefield (e.g., Toutain and Gouttebroze, 1993), as measured at the stellar surface r = R. The pressure perturbation p is related to the wavefield displacement ξ through the linearized adiabatic equation

where c₀ and P₀ are, respectively, the sound speed and the pressure of the unperturbed stellar model. The pressure perturbation p_nℓM(r, θ_β, ϕ_β) of the mode M is

where ξ_nℓM is given by Equation (3). To leading order, we have

with

where we dropped the subscripts nℓ. Acoustic oscillations in stars are stochastically excited and damped by turbulent convection, therefore I(θ_β, ϕ_β, t) is a realization of a random process. Since the perturbation is steady in the frame ${{R}}_{\beta }$, this random process is stationary in that frame. An expression for the intensity fluctuations I(θ_β, ϕ_β, ω) in Fourier space with the required statistical properties is

where the ${{N}}_M(\omega )$ are independent complex Gaussian random variables, with zero mean and unit variance:

In Equations (27) and (28) we only consider the positive-frequency part of the spectrum (ω and ω′ > 0) since I(θ_β, ϕ_β, t) is real. The negative-frequency part is related to the positive part by $I({\theta }_{\beta },{\varphi }_{\beta },-\omega )=I^{*}({\theta }_{\beta },{\varphi }_{\beta },\omega )$. The function L_M(ω) is a Lorentzian

appropriate for describing the power spectrum of an exponentially damped oscillator with full width at half maximum (FWHM) Γ (see e.g., Anderson et al., 1990).

We transform Equation (27) back to the time domain by inverse Fourier transformation to obtain I(θ_β, ϕ_β, t). The intensity I(θ, ϕ, t) as seen by the observer in the inertial frame ${R}$ is obtained by applying a passive rotation of Euler angles (0, −i, Ω_βt) to express $Y_{\ell }^m({\theta }_{\beta },{\varphi }_{\beta })$ in terms of θ and ϕ (Messiah, 1960):

In the frequency domain, the intensity fluctuations become

where we used the property $r_{m^{\prime }m}^{(\ell )}(-i)=r_{m{m}^{\prime }}^{(\ell )}(i)$. Since L_M(ω − mΩ_β) peaks at frequency ω = ω_M+mΩ_β, the intensity spectrum observed in the inertial frame has (2ℓ + 1)² peaks, corresponding to all combinations of m and M.

To obtain the full-disk integrated intensity fluctuations I_obs(ω) we perform an integration over the visible disk of the star:

where W(θ) is the limb-darkening function. The components with m′ ≠ 0 vanish upon integration over ϕ, thus

with $B_m^M=A_m^M{V}_{\ell }{r}_{m0}^{(\ell )}(i)$ and

The matrix elements $r_{m0}^{(\ell )}(i)$ are written explicitly in terms of the associated Legendre polynomials (Messiah, 1960):

A realization of the power spectrum is given by

which depends on 2ℓ + 1 independent realizations of complex Gaussian random variables. The expectation value of P(ω) is

As mentioned earlier, the power spectrum displays (2ℓ + 1)² Lorentzian peaks, centered at frequencies

with peak power

Example power spectra for a dipole and a quadrupole multiplet are shown in Figure 3. The frequencies and amplitudes of the (2ℓ + 1)² peaks are obtained from Equations (37) and (38), using the limb-darkening function quoted by Pierce (2000) to calculate V_ℓ.

Figure 3 also displays the different contributions to the frequency splittings due to rotation and to the active region perturbation in the corotating frame ${{R}}_{\beta }$. For both multiplets, the M = 0 peaks are shifted by the largest amount, they are the most affected by the AR perturbation and in the frame ${\hat{{R}}}_{\beta }$ (see also Papini et al., 2015). This feature, which arises from geometrical considerations only, is preserved in the spectrum as seen in the observer's frame, where the M = 0 peaks are clearly visible. With increasing AR surface coverage the frequency shifts of the M ≠ 0 modes increase and the peaks get less clustered.

3. Results

3.1. Dipole and Quadrupole Power Spectra

In this section we describe the changes imprinted by a large active region in the spectrum of two multiplets with (ℓ, n) = (1, 18) and (2, 18), for a star rotating with a period of 8 days.

Figure 4 shows the results for the quadrupole multiplet, for four combinations of the values of α and β selected in section 2.2: the observed power spectra are plotted for two angles of observation, i = 30° and 80°, and are normalized with respect to V₂, i.e., with respect to the m = 0 peak of the pure rotational spectrum seen with i = 0. The corresponding theoretical Lorentzian envelope (solid line) has been calculated from Equation (36), by setting a value for the FWHM of Γ/2π = 1 μHz, typical for this multiplet in the Sun (see e.g., Chaplin et al., 2005). Due to the finite lifetime of the modes of oscillation, it is clear that is not possible to resolve all the (2ℓ + 1)² peaks, and an observer would identify not many more than (2ℓ + 1) peaks in a multiplet. In the cases shown here, it is possible to identify from 5 to 6 peaks for i = 80°, the additional peak coming from the uppermost shifted m = 2, M = 0 peak. We note that the Lorentzian envelope displays an asymmetric profile. Because of their large shift in frequency, the M = 0 peaks blend with peaks from different m-quintuplets. This blending increases with activity level. Figure 4A shows a case for which the (M, m) = (0, 0) and (M, m) = (1, 1) peaks have close frequencies and comparable amplitudes, they contribute equally to a single peak in the power spectrum.

The envelope of the power spectrum is very sensitive to the latitudinal position of the AR and to the inclination angle: in Figure 4C the power spectrum is near the standard rotationally-split spectrum, while the same configuration observed from a different inclination angle (Figure 4D) shows a more asymmetric profile with additional peaks. This is better seen in Figure 5, which shows contours of the acoustic power as a function of inclination angle, for both the ℓ = 1 and ℓ = 2 multiplets, for the same active region parameters as in Figure 4. For an active region at high latitude (middle panels of Figure 5), the central peak shows a significant shift and overlaps with the m = 1 peak. For a low-latitude active region (bottom panels) the envelope of the power spectrum displays more than 2ℓ + 1 peaks. A distinct feature is the presence of two peaks instead of one when observing an ℓ = 2 multiplet at zero inclination angle. The sensitivity of the spectrum to the colatitude of the AR, shown in Figure 6, is due in part to the variation with β of the non-diagonal elements of the rotation matrix R^(ℓ).

An observed power spectrum is, of course, much more difficult to interpret than its expectation value. The power spectrum in Figure 7 includes realization noise due to the stochastic nature of stellar oscillations and to additional shot noise. At each frequency the observed power is a realization of an exponential distribution (a chi-squared with two degrees of freedom) with standard deviation and mean equal to the expectation value of the power spectrum. Realization noise considerably degrades the spectrum, however in some cases it is still possible to distinguish between the pure rotational spectrum and a spectrum with the active region.

3.2. Asymptotics

3.2.1. Limit of Small Latitudinal Differential Rotation

The examples shown so far suggest that the power spectrum of a multiplet may be approximated by the sum of 2(2ℓ + 1) peaks, for two reasons. First, in the corotating frame the splitting due to rotation is small compared to the shift induced by the AR perturbation (see e.g., Figure 3). Second, the AR perturbation induces a shift that is largest for the M = 0 mode. Therefore, for each m, all the M ≠ 0 peaks are clustered near the pure rotational frequencies and appear as a single peak, while the M = 0 peaks are well separated.

Here we wish to find an approximation for the power spectrum of a multiplet in terms of 2(2ℓ + 1) Lorentzians only. In the following we assume the rotation perturbation to be small compared to the active region perturbation, and seek an approximate solution to the eigenvalue problem (14).

It is convenient to solve the eigenvalue problem in the frame where the dominant perturbation is diagonal, i.e., the frame ${\hat{{R}}}_{\beta }$ of the active region. The following results are similar to section 19.5 of Unno et al. (1989). We rewrite the elements of the full matrix $\hat{\text{O}}$ in ${\hat{{R}}}_{\beta }$ as

where the term with the sum in the right hand side correspond to ${\hat{O}}_{m^{″}{m}^{\prime }}^{\Omega }$. We look for solutions of the form

and

where $\left(\delta {\omega }_M^{(1)},{\hat{A}}^{M,(1)}\right)$ are the perturbations to the (partially degenerate) eigenvalues $\delta {\omega }_M^{\text{AR}}$ and eigenvectors ${\hat{A}}^{M,AR}$ of ${\hat{\text{O}}}^{\text{AR}}$, with ${\hat{A}}_m^{M,\text{AR}}={\delta }_{mM}$. The eigenvector perturbation is orthogonal to the reference eigenspace (see e.g., Messiah, 1960, p. 687),

which gives ${\hat{A}}_{\pm M}^{M,(1)}=0.$ The eigenvectors A^M in the frame ${{R}}_{\beta }$ are given by

via the rotation matrix R^(ℓ) (see Equation 16).

To first order, the eigenvalue problem (Equation 14) becomes

where I is the identity matrix. To calculate ${\hat{A}}^{M,(1)}$ and $\delta {\omega }_M^{(1)}$ we multiply the above equation on the left by the transpose of ${\hat{A}}^{m^{\prime },\mathrm{\text{AR}}}$ to obtain:

We find the perturbed eigenvalues by setting m′ = M in the above equation

The non-zero elements of A^{M, (1)} are obtained from the m′ ≠ ±M components of Equation (45):

The explicit expressions for the eigenvalues and eigenvectors of the combined perturbations are

and

The explicit expression for the amplitudes $A_m^M$ in the frame ${{R}}_{\beta }$ is then.

3.2.2. Neglecting Latitudinal Differential Rotation

If we neglect differential rotation then Equation (18) reduces to

Then, by using the identities (Unno et al., 1989; Gough and Thompson, 1990)

and

the frequency shifts (Equation 48) simplify to

Note that for moderately fast rotating stars the Coriolis term in the above equation is much smaller than the centrifugal distortion term (e.g., Gizon and Solanki, 2004).

Due to the clustering of the M ≠ 0 peaks, the power spectrum in the observer's frame can be approximately modeled by 2(2ℓ + 1) Lorentzians, some of which may overlap. Half of these correspond to the peaks with M = 0. The remaining 2ℓ + 1 (approximate) Lorentzians are obtained by summing over the M ≠ 0 peaks; their mean frequency shifts are given by a power weighted average of the M ≠ 0 frequency shifts $\delta {\omega }_M^m=\delta {\omega }_M+m{\Omega }_{\beta }$. We denote by 〈δω〉_m this average:

with

The corresponding averaged power amplitudes 〈P〉_m are (see Equation 38)

Figure 8 shows how good is the 2(2ℓ + 1)-Lorentzian model in reproducing the expected power spectrum for the case of Figure 7. The vertical red lines are for the (2ℓ + 1) peaks with M = 0, while the vertical black lines refer to the (2ℓ + 1) peaks with power 〈P〉_m and frequency shifts 〈δω〉_m. The envelope of power of the 2(2ℓ + 1)-model compares well with the envelope obtained by summing over all the (2ℓ + 1)² peaks.

Figure 9 shows the frequency shifts 〈δω〉_m (black lines) and the shifts of the M = 0 peaks (red lines) calculated from the solution of Equation (14) (solid lines) and from the approximate solutions (dashed lines), as a function of stellar rotation rate and for an active region with a surface coverage of 4% (α = 23°) and β = 80°. The linear approximation successfully returns the frequency shifts of the 2(2ℓ + 1) Lorentzians, even for moderately high rotation rates.

3.2.3. Limit of Small-Size Active Region (α ≲ 15°)

The above formulae for the frequency shifts simplify further when the active region has a small surface area. For small values of α, the integral in Equation (21) can be approximated,

such that $G_{\ell }^{m^{″}}(\alpha )$ becomes

As shown in Figure 2, the above approximation is very good for α ≲ 15° and ℓ ≤ 2. Up to order α², the active region induces a shift only in the frequency of the M = 0 mode

Perturbed frequency shifts and amplitudes are obtained by following the calculation described in section 3.2.1, but taking into account the fact that the eigenfrequencies $\delta {\omega }_M^{\text{AR}}$ are now degenerate for M ≠ 0. The frequencies in the observer's frame can then be approximated by

where

The peak amplitudes in the power spectrum, $P_m^M$, are given by Equation (38), where $A_m^M$ is replaced by

with

Figure 10 shows the frequency error in μHz introduced by the small-α approximation. Dark blue shades indicate the regions in parameter space (α and Ω_β) where the approximation is very good. For an active region with a surface coverage below 2%, i.e., α <15°, the frequency shift ${\overline{\delta \omega }}_{M=0}$ is within ~ 0.1μHz of δω_{M = 0}, even for fast rotation rates.

3.3. Correlations in Frequency Space

Due to the fact that the active region perturbation is unsteady in the observer's frame, the intensity fluctuations (Equation 32) are not statistically independent in frequency space. Given two frequencies ω and ω′, the intensity covariance is

To simplify the analysis we assumed that $k=({\omega }^{\prime }-\omega )/{\Omega }_{\beta }$ is an integer (this is not a weakness of the theory though). The above expression vanishes unless |k| ≤ 2ℓ. The quantities I_obs(ω) and $I_{\mathrm{\text{obs}}}({\omega }^{\prime })$ are correlated for frequency separations $\Delta \omega ={\omega }^{\prime }-\omega =(m^{\prime }-m){\Omega }_{\beta }$.

The power spectrum is also correlated for $\Delta \omega =(m-m^{\prime }){\Omega }_{\beta }$. Using the formulae given in appendix C of Fournier et al. (2014), we find

where the intensity covariance is given above. We note that the values of the intensity and power spectrum covariances (Equations 65 and 66) depend on the parameters of the model. However the existence of a correlation is a general feature, which arises from the fact that the active region rotates with the star. Remarkably, the expectation value of the power spectrum ${P}$(ω) (see Equation 36) is as if all the terms in Equation (32) were statistically independent.

4. Discussion

4.1. Non-linear Frequency Shifts and Amplitudes From Numerical Simulations

In order to compare with the perturbation theory, we study the non-linear regime by means of numerical simulations. We use the GLASS wave propagation code, with the same numerical setup employed by Papini et al. (2015).

4.1.1. Rotation in Post-processing

Running different numerical simulations for different values of β and different perturbation amplitudes is computationally expensive. Instead we performed simulations for a 3D polar perturbation to the sound speed and for a star with no rotation, that is equivalent to solve the numerical problem in the reference frame ${\hat{{R}}}_{\beta }$. We introduced the effect of rotation later in processing the output. This approach has the advantage that, for a given amplitude of the AR perturbation, we only need to run one simulation in order to calculate the power spectrum for any given value of β and rotation period. However we can only reproduce solid body rotation, and it is not possible to include the effects of the centrifugal distortion and of the Coriolis force, therefore in each nℓ-multiplet we expect to find only (2ℓ + 1)(ℓ + 1) peaks. Nonetheless the results are useful for exploring the non-linear regime. We note that, as a consequence of neglecting these rotational effects, M can be identified with the azimuthal degree of the spherical harmonics $Y_{\ell }^M({\hat{\theta }}_{\beta },{\hat{\varphi }}_{\beta })$ in the frame ${\hat{{R}}}_{\beta }$ (see section 2.1.1).

4.1.2. Sound-Speed Perturbation

As was done earlier, we approximate the perturbation introduced by the starspot by a local increase in sound speed. We consider separable perturbations in the square sound speed of the form

where ϵ > 0 is a positive amplitude, f is a radial profile, and g a latitudinal profile. Explicitly, g = 1/2 + cos(κθ)/2 is a raised cosine for 0 ≤ κθ < π and is zero otherwise, where π/(2κ) = 0.65 rad = 37.5°. The function $f=exp(-|r-r_c{|}^2/2{\sigma }^2)\left[1/2+cos(|r-r_c|/\sigma )/2\right]$ is a Gaussian function centered at radius r_c = 0.9985 R_⊙ with dispersion σ = 0.004 R_⊙, multiplied by a raised cosine. This functional form is the same as the one described by Papini et al. (2015). The perturbation is thus placed along the polar axis at a depth of 4 Mm, with a surface coverage of 12%.

4.2. Synthetic Power Spectrum

As in section 2.3, we assume that the intensity fluctuations are proportional to the Eulerian pressure perturbation measured at r₀ = R + 200 km above the surface (see Papini et al., 2015). To calculate the expectation value of the observed power spectrum, we performed a first set of simulations for which all the modes were excited with the same phase at the initial time.

The approximation that the changes in the eigenfunction of a mode M are limited to the same angular degree ℓ (Equation 3) does not hold in the non-linear case: the horizontal shape of an eigenfunction is a combination of spherical harmonics of different ℓ values. However, since the perturbation is axisymmetric (see Equation 67), the eigenfunction for a mode M is a combination of spherical harmonics $Y_{\ell }^M({\hat{\theta }}_{\beta },{\hat{\varphi }}_{\beta })$ of different angular degree ℓ and same M. Therefore the intensity fluctuations $I_M({\hat{\theta }}_{\beta },{\hat{\varphi }}_{\beta },t)$ due to all the modes with the same M take the form

where p_ℓM(r₀, t) are the coefficients of the spherical harmonic decomposition of the pressure wavefield $p(r_0,{\hat{\theta }}_{\beta },{\hat{\varphi }}_{\beta },t)$ returned by the GLASS code in the frame ${\hat{{R}}}_{\beta }$. The index ℓ_max is set by the spectral resolution of the spherical harmonic transform.

To obtain the full-disk integrated intensity in the observer's frame we express each $Y_{\ell }^M({\hat{\theta }}_{\beta },{\hat{\varphi }}_{\beta })$ in terms of the spherical harmonics in the frame ${R}$

by means of two consecutive rotations of the Euler angles (0, −i, Ω_βt) and (0, β, 0). Combining this equation with Equation (68) and integrating over the visible disk, we obtain the full-disk integrated intensity of each mM-component:

We then perform a Fourier transform to calculate the intensity I_mM(ω) in the frequency domain. Finally, we derive the expression for the expectation value of the power spectrum

that is analogous to Equation (36), but for the entire wavefield.

For the non-linear study we chose three perturbation amplitudes ϵ = 0.1, 0.2, and 0.3 which, for the multiplet (ℓ, n) = (2, 18), correspond to ϵ_nℓ ≃ 0.005, 0.010, and 0.015, i.e., roughly twice to six times the value used in the linear analysis. The simulations run for 80 days (stellar time), in order to reach an accuracy of ~ 0.14 μHz in the frequency domain. The wavefield computed by GLASS includes some numerical damping that increases with frequency with an exponential dependence. We took advantage of this damping and selected two ℓ = 2 multiplets: one with n = 18 and a FWHM comparable to the solar value, the other with n = 12 and a FWHM small enough to resolve all the mM peaks in the multiplet. Model S is convectively stabilized (Papini et al., 2014), which implies that the unperturbed frequencies of the quadrupole modes are 1970.50 μHz for n = 12 and 2783.62 μHz for n = 18. Figure 11 shows the observed power spectra of the two selected multiplets in the case β = 80° and i = 80°, normalized with respect to the highest M ≠ 0 peak. In the case ϵ = 0.1, the simulated power spectrum (black curve) of the (ℓ, n) = (2, 12) multiplet (top left panel) is well reproduced by the power spectrum computed with linear theory (red dot-dashed curves), except for the peaks corresponding to M = 0 (vertical red curves show the M = 0 peaks from linear theory), which are less shifted in frequency and have smaller amplitudes than predicted. For the (ℓ, n) = (2, 18) multiplet (top right panel) the non-linear effects are less visible, due to the overlapping Lorentzian profiles. A blue curve, displaying the contribution to the power spectrum of the M = 0 peaks, shows that also for this multiplet the M = 0 component of the power spectrum deviates from the linear behavior, both in frequency and amplitude. This plot also shows an example of the combined action of mode mixing and mode visibility, which almost suppresses the m = 0, M = 0 peak located at a frequency shift of ~ 3.1 μHz.

With increasing ϵ (middle and bottom panels) the interaction of the wavefield with the active region becomes strongly non-linear, and for a perturbation with ϵ = 0.3 results in a massive distortion of the power spectrum with respect to that one predicted by linear theory. Here two different behaviors are evident: in the (ℓ, n) = (2, 12) multiplet the m = 0, M = 0 peaks are almost suppressed, while the peak with m = 0, M = 0 in the (ℓ, n) = (2, 18) multiplet, that was suppressed in the case ϵ = 0.1, increases in amplitude as ϵ increases (middle and bottom right panels). Moreover the peaks with M = 1 start to deviate from the linear prediction. This is in agreement with what found in the non-rotating case by Papini et al. (2015), who also showed that second-order perturbations would correct most of the differences.

4.3. Correlations in Synthetic Power Spectra

In section 3.3 we showed that in presence of an active region rotating with the star, the power spectrum of a multiplet is correlated at frequency separations that are multiple of the rotational frequency of the active region. However, for small perturbations this correlation is too weak to be observed. Here we ask whether such a correlation can be measured in presence of a perturbation of moderate amplitude. For that purpose, we ran a second set of simulations, in which the acoustic waves were excited by applying a random forcing function at 150 km below the surface at each time step, as described by Hanasoge and Duvall (2007). The duration of the simulation is 80 days (stellar time). In the frequency domain, the observed intensity is

where I_mM(ω) is obtained by Fourier transformation of Equation (70) using the numerical realizations of p_ℓM(r₀, t).

The autocorrelation of the intensity spectrum is

where [ω_nℓ] denotes an appropriate frequency interval of size ~ 20 μHz containing the multiplet nℓ but excluding the nearby l = 0 mode. The average is performed over all ℓ = 2 multiplets with n ranging from 15 to 25.

Figure 12 shows the real part of C_I for the special case of an inclination angle i = 80° and a starspot at colatitude β = 80° with ϵ = 0.1. A correlation is clearly visible at frequency separations Δω = 2Ω_β and 4Ω_β. This suggests that the frequency-domain autocorrelation function could be used as a diagnostic tool to identify unsteady perturbations in the time series of stellar oscillation. We note that the imaginary part of C_I contains no visible signal above the noise level.

4.4. Toward a Physical Model for Mode Interaction With an Active Region

In this paper we replaced the active region by a localized increase in sound-speed near the stellar surface. We focused on the geometrical aspects of the problem rather than on the physics. One may ask, however, what would be the difference in the obtained results if we had instead considered a realistic model for the magnetic active region. Although we will not answer this question here, it is worth listing some of the steps involved.

A typical solar active region consists of a pair of sunspots surrounded by plage with strong vertical field. Local helioseismology of the visible disk and the far side indicates that p modes are strongly scattered by both sunspots and extended plage (see e.g., Gizon et al., 2009, and references therein). Some studies of the interaction of high-degree p modes with magnetostatic sunspots (e.g., Moradi et al., 2010; Cameron et al., 2011) have been carried out using MHD wave propagation codes (Cameron et al., 2008; Felipe et al., 2016). Other studies are based on numerical magneto-convective simulations (e.g., Rempel et al., 2009; DeGrave et al., 2014). The main conclusion of these simulations is that the interaction takes place in the top few hundred kilometers below the surface, where the direct effects of the magnetic field and the indirect effects due to changes in thermodynamic structure with respect to the reference atmosphere (e.g., the Wilson depression) are large. Wave simulations indicate that the outgoing p modes are phase shifted with respect to the incoming p modes in such a way that the effective wave speed is increased, as observed. The physical interaction involves the conversion of p modes into fast and slow magnetoacoustic modes in the sunspot (e.g., Khomenko and Collados, 2006; Cameron et al., 2008). A fraction of the incoming p-mode energy is tunneled downward in the form of slow MHD waves, leading to absorption (Braun, 1995; Zhao and Chou, 2016). See also, e.g., Saio and Gautschy (2004) and Cunha (2006) for mode conversion calculations in roAp stars.

Using 2D ray tracing, Liang et al. (2013) showed that high-degree helioseismic waveforms can be reproduced by increasing the effective wave speed by 10% in the sunspot. This provides some justification for the values that we have used in the present paper, although the extension to low-degree p modes has not been studied. Clearly, much additional work will be needed to determine the correct active-region perturbation amplitude from first principles. Until then, a simple calibration can be obtained from the observational study by Santos et al. (2016) who estimated empirically the contribution of sunspots to the low-degree p-mode frequency shifts associated with the solar cycle. By combining our Equation (20) with Equation (3) from Santos et al. (2016), we find ε_nℓ = −Δδ_ch/I_ℓ, where Δδ_ch is the integral phase difference introduced in the mode eigenfunction by a sunspot and I_ℓ is related to mode inertia. Using the value Δδ_ch = −0.44 estimated by Santos et al. (2016), we have ε_nℓ ~ 0.05 for a quadrupole mode, i.e., a larger value than proposed by Liang et al. (2013) and used in the present paper. This may suggest that for large active regions the spectra calculations may have to be carried out in the non-linear regime. However, only realistic numerical modeling would help settle this question. The full problem would also have to include multiple scattering by collections of flux tubes in plage (see e.g., Hanson and Cally, 2015).

5. Conclusions

In this paper we investigated the changes in global acoustic oscillations caused by a localized sound-speed perturbation on a rotating star mimicking a large active region, using both linear perturbation theory and 3D numerical simulations. In an inertial frame, the active region perturbation is unsteady.

We find that the power spectra of low-degree modes have a complex structure. The combined effects of the active region and differential rotation cause each nℓ-multiplet to appear as (2ℓ + 1)² peaks, each with a different amplitude. Most of the peaks are clustered near the classical rotationally-split frequencies, and only 2ℓ + 1 peaks (the M = 0 peaks, which correspond to the axisymmetric mode in the reference frame of the AR) are shifted to higher frequencies. This leads to an apparent asymmetry in the line profiles. However, due to the finite lifetime of acoustic oscillations, most of the peaks cannot be resolved. For solar-type stars, the results are not very sensitive to the choice of latitudinal differential rotation profile.

The structure of the power spectra is sensitive to the latitudinal position of the active region and to the inclination angle i of the stellar rotation axis. The latter plays a major role in determining the relative visibility of the individual peaks. We find that the envelope of the power spectrum becomes more complex as the latitude of the active region decreases. In practice, it would be very difficult to perform a fit of the (2ℓ + 1)² peaks in a multiplet, due to peak blending and noise. However, by neglecting differential rotation it is possible to derive a simplified formula that approximates the power spectrum of a multiplet to a sum of only 2(2ℓ + 1) Lorentzian profiles. For small-area active regions, the formula further simplifies and directly links the frequencies of the peaks in the power spectrum to the active region parameters. Such formula may find applications in the analysis of real asteroseismic observations.

Numerical simulations were performed to explore the non-linear regime of the perturbation. We find that the M = 0 peaks deviate from the linear behavior for active-region perturbation amplitudes ε_nℓ ≳ 0.005. Depending on each particular case, the amplitude of these peaks is either reduced or enhanced compared to first-order linear theory, due to mixing with modes with other values of ℓ and m (Papini et al., 2015).

We found that there are correlations in the power spectrum at frequency separations that are multiples of the active-region rotation rate. In the linear regime the correlation signal is too weak to be observed. However the numerical simulations show that for active-region perturbations of moderate amplitude, such a correlation might be detectable, provided that the frequency intervals are carefully selected to increase the signal to noise ratio.

We note that perturbation theory can easily be extended to compute the effect of multiple active regions provided that latitudinal differential rotation is small. The treatment of several active regions rotating at different rotation rates would require a different setup, since there is no frame in which these perturbations are steady. A numerical approach would be preferable in such a case.

The work presented in this paper uses simplified physics, but it should provide useful guidance to identify the seismic signature of a large active region in the power spectrum of stellar oscillations. Given that the values of ϵ_nℓ are uncertain, we believe that it is worth searching for low-degree multiplets consisting of 2(2ℓ + 1) components in available asteroseismic observations. Ideal targets are stars that are known to have high-quality oscillation power spectra (high SNR, narrow lineprofiles, clear rotational splitting, see e.g., Nielsen et al., 2014) and show evidence for long-lived starspots (e.g., Nielsen et al., 2013, 2019). The catalog of potential targets, currently limited to CoRoT and Kepler, will increase fast with TESS (Ricker et al., 2015) and PLATO (Rauer et al., 2014).

Author Contributions

LG provided the basic theory. EP provided the numerical applications. Both authors contributed to the analysis of the results and to the writing of the article.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This paper is a contribution to PLATO PSM activity Seismic diagnostics of stellar activity. We acknowledge financial support from the German Aerospace Center DLR and from the Max Planck Society. We thank Shravan Hanasoge for making the GLASS code available. Computational resources were provided by the German Data Center for SDO. This work appeared in part in the Ph.D. thesis of Papini (2016).

References