## Abstract

We have searched helioseismic data collected by the Birmingham Solar Oscillations Network (BiSON) for solar-cycle changes to those low-ℓ p-mode parameters that relate to the excitation and damping of the resonances. These data — collected between 1991 and 1997 — cover the complete declining phase of solar activity cycle 22 (up to and including the cycle 22/23 boundary). Over the range 2600≤ν≤3600 μHz, we uncover a mean 24±3 per cent increase in the frequency-domain linewidths; a mean decrease of 46±5 per cent decrease in the mode heights, and a mean decrease of 22±3 per cent in the modal velocity powers. The rate at which energy is supplied to the modes remains constant, at the level of precision of the observations (measured change 0±4 per cent). We use expressions derived from the equation of a damped harmonic oscillator to illustrate the diagnostic properties of the observables: these indicate that both the signs and relative sizes of the extracted variations can arise from changes solely to the net damping; the net forcing of the modes need not change.

The results possibly hint at the changes being maximal at frequencies near ∼3100 μHz. They might therefore suggest an origin for the observed variations that is peaked in the superadiabatic layer of the convection zone, which couples most strongly to the eigenfunctions of modes at the centre of the p-mode spectrum.

## Introduction

The study of those aspects of helioseismology that relate to the excitation and damping of the acoustic resonances of the Sun has recently attracted a great deal of interest. The quality of data being collected by both modern satellite and ground-based instrumentation (e.g., VIRGO, MDI and GOLF on *SOHO*, and the BiSON, GONG and LowL/ECHO networks) are now such that subtle, hitherto-unseen effects are becoming apparent. Perhaps the most conspicuous example is the presence of skew-symmetry in the observed resonant mode profiles. This arises as a natural consequence of the modes being excited by a source whose radial extent is small relative to that of the principal axis of the resonant cavity. Interference along different acoustic ray paths emanating from the source then gives rise to phase shifts in the re-constructed interference signal and asymmetry in the resonant peaks in the frequency domain. A careful study of the complex resonant forms that can arise — not only in the vicinity of the peaks of the resonances, but also in the intervening minima — uncovers information regarding the localization of the source relative to the outer boundary of the cavity and its nodes (see, e.g., Gabriel 1995, Abrams & Kumar 1996, Rast & Bogdan 1998, Nigam et al. 1998, Rosenthal 1998 and Kumar & Basu 1999a). Recent analyses of this type place the source in the outer few hundred kilometres of the convection zone (see, e.g., Chaplin & Appourchaux 1999, Chaplin et al. 1999, Kumar & Basu 1999b and Nigam & Kosovichev 1998, 1999).

Further, such studies can also reveal how the modes interact with their environment in the neighbourhood of where the excitation occurs. There is now an increasing body of evidence to indicate that observations of the evanescent tails of the bounded resonances are ‘contaminated ’ by a correlated component of the convective granulation (which is responsible for exciting the modes); it is actually a component of the convective overshoot present in the photosphere which imprints itself upon the oscillation observations, and in so doing modifies the phase properties of the resonant signal to such an extent that the shape of the observed profiles can be noticeably altered (Roxburgh & Vorontsov 1997; Nigam et al. 1998; Kumar & Basu 1999a, b; Nigam & Kosovichev 1999). The relative size disparity between the imprinted signal that acts when either Doppler velocity or intensity observations are made is now believed to explain the opposite sign of asymmetry that is usually revealed in the two types of data.

Theories of the excitation of p modes by stochastic turbulence indicate that the forcing is peaked in the outer layers of the convection zone (see, e.g., recent simulations by Stein & Nordlund 1998). The fact that analyses of the type discussed above, together with observations of individual ‘acoustic ’ events (Rimmele et al. 1995; Espagnet et al. 1996; Goode et al. 1998) also lead to similar conclusions lends credence to stochastic theories. They are also sustained by observations of the distribution of observed mode powers, which largely follow the expected Boltzmann distribution (e.g. Elsworth et al. 1995; Chaplin et al. 1997; Chang & Gough 1998), although some question marks remain concerning the origin of occasional very high-excitation events (Foglizzo et al. 1998).

In this paper we concern ourselves with another aspect of excitation and damping phenomenology, namely that of the variation of the mode linewidths, amplitudes, energies and energy supply rates, on the time-scale of the solar activity cycle. Evidence of changes to these parameters can in principle provide important clues to the physical origin of the solar activity cycle. Studies in the literature show clear evidence for a decrease in the amplitudes of the acoustic modes with increasing solar activity (e.g. Palle et al. 1990a; Anguera-Gubau et al. 1992; Elsworth et al. 1993); evidence for changes to the modal linewidths is sketchy at best, and in some cases contradictory (Palle et al. 1990b; Jefferies et al. 1990, 1991; Meunier 1997). This is perhaps not surprising, given that the task of extracting reliable estimates of the linewidths is fraught with difficulties. Further, the intrinsic precision with which the width and amplitude parameters can be measured is appreciably inferior to that achievable in the mode frequencies. As such, efforts to uncover variations in them demand the use of long, high-quality, well-filled data sets.

Here, we present and discuss in detail what we believe has been a successful attempt to uncover changes to the linewidths of low-angular-degree (low-ℓ) p modes over the falling phase of solar activity cycle 22. To do this, we have used several high-quality data sets collected by the Birmingham Solar-Oscillations Network (BiSON). We also confirm previous reports of a reduction in the strength of the modes at high activity; in addition, our results are not inconsistent — at the level of precision of the data — with the premise that the net rate at which energy has been pumped into the modes has remained constant over the time frame covered by the data.

While we have confidence in the general conclusions drawn in the paper, we nevertheless stress that the analysis is difficult: this is because of the presence of bias in estimates of the p-mode parameters which varies according to the nature of the temporal window functions of the individual data sets, and the varying quality of the data across those sets. We therefore devote much of the paper to a discussion of the origin of these effects, and attempts to allow for them in the analysis of the BiSON data. However, we begin below by drawing the analogy of a damped oscillator in order to derive some simple, formalistic expressions for the various parameters under consideration; we then use these to illustrate the diagnostic properties of the observed solar p-mode parameters, before presenting and discussing in some detail the analysis procedures and results.

## Excitation and damping phenomenology

### Analogy of a damped oscillator

A useful analogy for the p modes is a forced, damped harmonic oscillator, of the form:

In the above, *x*(*t*) is the displacement of the oscillator, ω_{0} is its natural angular frequency, η is the damping constant, and *f*(*t*) is the forcing function. The Fourier transform of the oscillator equation gives the shape of the expected power spectrum of the velocity signal in the vicinity of the resonance (ω≈ω_{0}), i.e.

*F*(ω) — the power spectrum of

*f*(

*t*) — is a slowly varying function of ω, and η≪ω

_{0}. This is essentially a Lorentzian profile, with a radian fwhm of 2η. The peak height is

The total velocity power of a given mode is proportional to width times peak height, i.e.,

To a good approximation, the solar acoustic spectrum should be represented by an ensemble of such oscillators. The energy (kinetic plus potential) of a mode with associated mass *M* is given by

The rate at which energy is supplied to the modes, d*E*/d*t*, is readily derived by again having recourse to the oscillator analogy. The energy of a damped oscillator can be expressed as

*A*is the amplitude of the motion. It follows that and, taking derivatives,

If we combine equations (3), (4) and (5), we have

### What do changes to the mode parameters tell us?

With reference to the above equations, we see that the various parameters provide us with different phenomenological information. For example, measurements of the modal linewidths, Δν, might be expected to give a direct measure of the damping rate, *η;* the velocity power, *V*^{2}, in the modes reflects the balance between the excitation — as expressed by *F* (ω) — and the damping, while the energy supply rate, , provides information regarding the forcing function, *F* (ω). The use of Δν and in any analysis is particularly appealing, since these measures are independent of one another. So, given possible variations in η and *F*(ω), what changes to the ‘observables ’ (Δν, *H*, *V*^{2}, *E* and ) might we expect to see?

In the following discussion we have assumed that ω and *M* are constant at each (ℓ,*n*). Changes to the p-mode eigenfrequencies over the solar cycle are well documented (e.g. Elsworth et al. 1990; Libbrecht & Woodard 1990, 1991; Elsworth et al. 1994; Regulo et al. 1994), with the fractional magnitude of the variations increasing to δω_{0}ω_{0}∼1×10^{−4} at frequencies of ∼4000 μHz. [At frequencies above this, there is strong evidence to suggest that the sense of the trend reverses and becomes somewhat larger in magnitude at progressively higher *n* (e.g. Ronan, Cadora & LaBonte 1994; Elsworth et al. 1990; Libbrecht & Woodard 1990, 1991; Elsworth et al. 1994; Regulo et al. 1994; Chaplin et al. 1998a). Here, we aim to study changes to the modes at frequencies below 4000 μHz only.] The choice of whether or not to disregard changes to ω_{0} in our analyses clearly depends upon the inherent precision that we can achieve in the various ‘observables ’, given the quality of the 8-month BiSON spectra. While the frequencies of the modes can be determined to typical accuracies of a few parts in 10^{5} in such data, the fractional precision obtainable in *H* and Δν (which represent the two fundamental observables under consideration here) is of the order of 10 to 20 per cent or worse. We can therefore safely assume that δω_{0}ω_{0}≈0. The issue of whether *M* varies over the solar cycle depends upon possible changes to the modal eigenfunctions. Theoretical calculations — based upon variations in the granular size and mixing-length parameter — indicate that fractional changes of the order of ∼5 per cent are possible (Houdek, private communication). However, as we shall see, this is at or below the level of those changes to the ‘observables ’ that we are sensitive to, given the inherent quality of the data (and only then as a result of averaging changes over the main part of the spectrum to reduce uncertainties). Therefore, since we would not be able to detect variations in *M* of the size discussed above, we assume *δM**M*≈0. We note also that in the light of the equations in Section 2.1, fractional variations in *V*^{2} and *E* are identical (we list those for *V*^{2} only). So

These expressions indicate that changes to the observed velocity power and energy will result from changes to the damping rate, regardless of variations in the forcing; the energy supply rate remains fixed, since it is independent of the damping. Changes to the forcing alone will give rise to variations in the modal powers, energies and energy supply rates, while the linewidths remain unaltered.

We now seek to extract — and then interpret — variations in the observables introduced above from an analysis of several BiSON power spectra which span the maximum to minimum phase of a solar cycle.

## Analysis

### Data

We have analysed 10 contiguous 8-month power spectra, generated from data collected by BiSON over the period 1991 January 7 to 1997 September 1. These data cover the complete declining phase of solar activity cycle 22, from an activity maximum through to the minimum at the cycle 22/23 boundary. The duty cycles of the time series used to construct the spectra span the range ∼50 to ∼80 per cent. The fractional fill of each data set, the resulting interaction of its associated window function with the acoustic mode profiles in the frequency domain, and the inherent quality of the data, all have a direct impact upon the magnitude of the parameters extracted by the mode-fitting procedures. This demands great care in the analysis, as we shall discuss at length below.

### Extraction of basic mode parameters

Adjacent ℓ=2 and 0, and ℓ=3 and 1 modes were fitted in pairs in each power spectrum as rotationally split, skew-symmetric multiplets with an associated flat background offset. The first, temporal sidebands were also included in the model. The power, χ, in each azimuthal mode component was modelled according to (Nigam & Kosovichev 1998):

where and ν_{0}is the frequency of the Lorentzian component, Δν its width, and

*H*its height. The parameter α is essentially the ‘fractional ’ asymmetry of the mode. This formalism reduces to a pure Lorentzian for α→0. Note that it is an approximation to a function obtained by solving a suitable, bounded wave equation, and is valid only in the vicinity of the resonance.

We used a multidimensional direction-set minimization algorithm (Press et al. 1992) to perform the fitting, by maximizing an appropriate log-likelihood function consistent with a χ^{2} distribution with 2 d.o.f.. Since there are data set gaps present in the time domain, the bins in the frequency domain are not statistically independent. We have performed extensive Monte Carlo simulations with realistic artificial data — which consequently allowed us also to study the effects of correlations between overlapping components at higher frequencies — in order to verify that even though our assumption of statistically independent χ^{2} is not strictly correct, the results extracted are not noticeably affected (see Chaplin et al. 1998a).

For frequencies up to ∼3600 μHz, the following parameters were varied iteratively until they converged to their best-fitting values:

- (i)
a frequency for each mode;

- (ii)
a single skewness parameter (the

*m*components of the modes in the pair are all assumed to have the same skewness); - (iii)
a single parameter describing a symmetric rotational splitting pattern for each mode (not applicable for radial modes);

- (iv)
a single width for both modes (i.e., each

*m*component was assumed to have the same width in both multiplets constituting the low-ℓ pair); - (v)
a single height — that of the outer, sectoral

*m*components — for each mode; the relative*m*-component height ratios were assumed to take fixed values (Christensen-Dalsgaard 1989); - (vi)
a single relative sideband height for each mode; at low frequencies, the relative sideband heights were fixed according to the characteristics of the window function, and

- (vii)
a flat, background offset for the whole fit.

The natural logarithm of the height, width and background terms were varied — not the straightforward parameter values themselves — in order to give a quasi-normal fitting distribution. Estimates of the formal uncertainties in each parameter could then be determined from the inverse of the Hessian matrix describing the hypersurface of the fit.

At frequencies above 3600 μHz the fractional sideband heights were fixed, according to the characteristics of the window function and the rotational splitting pattern was fixed for all modes, with a constant synodic spacing (appropriate to Δ|*m*|=1) of 400 nHz. (For a more in-depth discussion of the fitting of low-ℓ data, see Toutain & Appourchaux 1994; Elsworth et al. 1990; Libbrecht & Woodard 1990, 1991; Elsworth et al. 1994; Regulo et al. 1994,and Chaplin et al. 1998a.)

In this paper we are concerned chiefly with studying variations to the fitted heights and widths, and to other combinations of both parameters. Since the fitting returns best-estimates of the natural logarithms of each, where we define

andwe choose to develop the following analysis in terms of the logarithmic values. Extracted variations in *w*(ℓ, *n*) and *h*(ℓ, *n*) (or other additive combinations of the two) therefore correspond to fractional variations in Δν(ℓ, *n*) and *H*(ℓ, *n*) (or other multiplicative combinations of the two) respectively.

The logarithm of the velocity power (velocity amplitude squared) of each mode is given by

The factor *k*_{obs} consists of several terms (i.e., a suitable normalization factor to give the integrated power under each modal peak, a correction for the ℓ-dependent filter of the observations, and corrections for limb darkening and Doppler imaging; see Chaplin et al. 1998b for a full discussion). With reference to equations (4) and (5), the natural logarithms of the energy and energy supply rates are given by

*n*) is the logarithm of the associated modal mass. Here, our analysis aims to uncover temporal changes to the parameters: if we assume that both

*k*

_{obs}and ϱ(ℓ,

*n*) (see Section 2.2) take constant values across the 10 spectra at each (ℓ,

*n*), they can then both be ignored in any differential analysis.

The determination of the uncertainties on υ^{2}(ℓ, *n*), *e*(ℓ, *n*) and (ℓ, *n*) must be approached with a little caution, owing to the fact that there is a strong degree of correlation between the estimates of the heights and widths extracted by the fitting procedure. To illustrate, in Fig. 1 we show those values of *w*(ℓ, *n*) and *h*(ℓ, *n*) returned by fitting some 250 independent realizations of the same artificial low-ℓ mode pair: we plot the estimates returned for the radial mode. (Here, we constructed artificial ℓ=2, *n*=19, ℓ=0, *n*=20 pairs; see Section 3.4 for an in-depth discussion of the Monte Carlo simulations.) The correlation coefficient ρ(*w*, *h*) between the fitted parameters is of the order of −0.9, and this must be taken into account when constructing the errors in the velocity power and energy parameters from those in *w*(ℓ, *n*) and *h*(ℓ, *n*). For example, consider the uncertainty in υ^{2}(ℓ, *n*). If σ_{w(ℓ,n)} and σ_{h(ℓ,n)} are the uncertainties in *w*(ℓ, *n*) and *h*(ℓ, *n*) respectively, then the usual combination-of-errors formula must be modified according to

Similar corrections are required to the uncertainties in *e*(ℓ, *n*) and (ℓ, *n*).

### Solar-cycle variation analysis of extracted parameters

With estimates of *w*(ℓ, *n*) and *h*(ℓ, *n*) extracted by the fitting procedures from each of the 10 spectra, we then proceeded to analyse these data — and suitable linear combinations of them (cf. equations 8, 9 and 10) — for significant variations. The analysis for a chosen parameter (i.e., logarithm of width, height, velocity power, energy or energy supply rate) was made in the following manner.

At each (ℓ, *n*), we constructed a mean data set from a weighted average of the parameter estimates from all 10 spectra. In all cases, we used the formal uncertainties returned by the fitting procedure to fix the weights. Residuals were then calculated — with respect to the values in the average set — for each datum in each spectrum. Separate linear fits were then made, again at each (ℓ, *n*), between the 10 parameter residuals and the 10 corresponding values of the mean 10.7-cm radio flux (in rf units) observed over the duration of the time series used to generate each spectrum. (Note that 1 rf unit corresponds to 10^{−22} W m^{−2} Hz^{−1}.) The use of an activity-related measure, rather than the epoch of the observations, serves to provide a more ‘physical ’ parameter for the dependent variable in the fits. Here, we note that the precision in the measured parameters is such that the use of higher order functions to parametrize any shifts is inappropriate. The returned best-fitting gradients from the linear fits serve as a measure of the observed change in the parameter per unit change in activity (here, per unit rf).

### Bias in the analysis: verification with artificial data

#### Origin of bias

Bias in the estimates of the widths and heights arises chiefly from a combination of the effects of the window function, and the inherent quality of the data set.

The characteristics of the temporal window function generated by a multistation network can be quite complex. Its interaction in the frequency domain will give rise to sidebands, owing to the presence of a quasi-diurnal gap structure in the time domain, and a complex background in the vicinity of the resonance, owing to the quasi-random gap structure that results from inclement weather and occasional instrumental breakdown. These effects are illustrated in Fig. 2. Here, we used a randomly forced, damped-oscillator model — which gives a Lorentzian resonance in the frequency domain (see Chaplin et al. 1997) — to generate 1000 independent realizations of the same, underlying ℓ=2, *n*=18, ℓ=0, *n*=19 mode pair (with the basic limit characteristics chosen to match those obtained for real data). Each individual *m* component was generated in turn, and the separate time strings co-added to give a combined, mode-pair time series. We then obtained the data represented by the solid line in Fig. 2 by modulating each combined time series by a representative 8-month BiSON window function (with duty cycle 75 per cent), and adding the 1000 power spectra computed from each of the modulated sets of residuals. The dashed line illustrates the mean spectrum obtained without multiplying through by the window. (We have scaled the two mean data sets to the same peak power level in the plot.)

By adding many independent realizations of the same basic mode, the resulting summation should tend toward the underlying limit spectrum. As expected, the dashed line mimics the anticipated combined-Lorentzian profile. The solid line illustrates clearly the effects of the window function, in particular the complex background that can arise in the immediate vicinity of each resonance. While the fitting model takes care of the sidebands, it does not adequately reflect all the additional complex structure. As a result, those widths extracted by the fitting procedure overestimate the true widths. Further, because of the width-height anticorrelation discussed in Section 3.2 there is a corresponding tendency for the extracted heights to be underestimated.

For the levels of data fill we are concerned with here — which range from roughly 50 to 80 per cent in the 10 BiSON data sets — the precise size of the bias introduced does, unfortunately, depend upon the characteristics of the window function and not merely the associated duty cycle. At higher levels of fill, the corrupting effects of the window are somewhat less noticeable, and the resulting bias is a better-behaved, more-monotonic function of fill alone.

If we ignore for the moment variations in signal-to-noise ratio (S/N) that might arise from changes to the strength of the modes over the solar cycle, then the S/N can also vary as a result of differences in the inherent noise characteristics of the various data sets — resulting changes to the background noise level of the spectra have a more noticeable effect upon the S/N of the weaker, narrower, lower frequency modes, and variations in the properties of the window function — i.e., in a data set with a low duty cycle, a large amount of power can be redistributed (or aliased) from the main, resonance peak into the window-dependent background, thereby reducing the S/N in the mode. Provided reductions in S/N are sufficient to begin to ‘mask ’ the wings of the mode — and any window-dependent background structure (cf. Fig. 2) — there will be a tendency for the fitting procedure to underestimate the true linewidths.

Our analysis of the BiSON data entails studying 10 independent spectra. Each has been modulated by a distinct window function; further, the inherent average quality of the data varies to some extent across the 10 sets. By and large, those spectra generated from data collected during the early part of the time frame considered here, which correlate with higher levels of solar activity, are of poorer quality and lower fill.^{1} Consider then the bias that might be present in any search for temporal changes in the data. Naively we might expect the outcome of any analysis to uncover suggestions of an increase in the widths (and decrease in heights) with increasing activity where there may actually be none, owing to the window-dependent bias. However, this does not take into account the counteracting effect of the poorer S/N in those spectra that coincide with elevated activity levels — this will tend to reduce the linewidths. Clearly, we need to map out the interaction of these complex, competing effects upon the differential-parameter analysis procedure outlined in Section 3.3, and to do so we have performed an extensive series of carefully planned Monte Carlo simulations. We have adopted the strategy of constructing ‘complete ’ artificial power spectra. These are designed to mimic, as closely as possible, the characteristics of the 8-month spectra presented by the full-disc BiSON observations, as modulated by the (ℓ, *m*)-dependent filter of the BiSON instrumentation. These artificial data can be processed in the same manner as the real data, thereby providing as thorough a check as possible of the analysis procedures.

#### Construction of artificial data sets

We began the construction of each artificial set in the frequency domain. Each underlying limit spectrum consists of some 170 separate (ℓ, *n*) multiplets, covering the range 1000≤ν≤5200 μHz and 0≤ℓ≤5. All *m* components to which the full-disc observations are sensitive are represented explicitly in the spectrum as Lorentzian peaks. As previously noted, the observed acoustic resonance profiles are slightly skew-symmetric. While this is sufficient to shift noticeably the estimated frequencies depending upon the model assumed for the limit profiles (i.e., an asymmetric model or a Lorentzian profile), the impact upon estimates of the widths and heights is not of any concern in the analyses here, given the precision achievable in their determination. Therefore, even though we have used a Lorentzian model for the modes, we nevertheless still expect the quantitative results to be relevant to the problem of analysing the real solar data. The Nigam & Kosovichev (1998) formalism is valid only in the vicinity of the resonance it describes, and as such cannot be used to construct a ‘complete ’ spectrum, since it gives rise to unrealistically large power levels in the tails of each mode. [In the parlance of equation (7), the power level χ(ξ) tends to the constant offset *Hα*^{2} far from resonance, i.e., for ξ≫1.] A more realistic spectrum could be constructed by making use of the expression from which the Nigam & Kosovichev formalism is derived, i.e., the Green's function of a simple potential-well model which serves to describe the essential elements of the solar acoustic cavity. However, this is somewhat more computationally expensive than simply adopting a Lorentzian model for each mode. (Naturally, we chose to fit the BiSON spectra to the Nigam & Kosovichev formalism in order to extract more reliable estimates for the frequencies.)

We used a data base of mode-parameter estimates obtained from analyses of BiSON spectra in order to fix the characteristics of the modes across each spectrum (e.g., frequencies, widths, S/N, etc.). For all multiplets, we constructed a symmetric, rotationally induced splitting pattern with a synodic spacing between adjacent *m* components of 400 nHz. The *m*-component height ratios in multiplets with ℓ=2 were fixed at the values determined by Christensen-Dalsgaard (1989). Further, we also used these calculations to establish appropriate visibility levels for the barely detectable ℓ=4 and 5 modes. At the extreme ends of the modelled spectrum — where we do not have reliable, fitted estimates for the parameters — we made appropriate, low-order extrapolations of the parameter values from the adjacent frequency bands.

Once a limit spectrum had been constructed, its real and imaginary parts were generated by applying Gaussian noise to the limit values (now in amplitude, rather than power). The corresponding time series was recovered by taking the inverse Fourier transform of the complex spectrum, and then modulated by one of the 8-month BiSON window functions. Additional normally distributed noise was also added in the time domain in order to give realistic signal-to-noise ratios at low frequencies. The sample standard deviation of the noise source was chosen to match the inherent quality of data associated with the selected BiSON window. The power spectrum of the time series could then be fitted, as per the real data (albeit with a Lorentzian model assumed for each peak).

In order to mimic fully the analysis procedure outlined in Section 3.3, each ‘subset ’ of Monte Carlo simulations required the construction of 10 complete artificial spectra, each modulated by one of the 10 BiSON window functions. Once generated, each 10-spectrum subset was then subjected to the ‘variation ’ analyses outlined in Section 3.3. We also ‘programmed in ’ changes to the widths and heights to some of the 10-spectrum subsets in order to test the response of the analysis to those variations. By performing an extensive series of such simulations, we were then in a position to understand, and hence to discriminate against or correct for, systematic bias in the analysis procedure.

#### Results of analysing artificial data sets

Here, we begin by detailing the results of performing an extensive series of independent, 10-spectrum-set simulations, for which the limit (i.e., input) widths and heights were unchanged at each (ℓ,*n*) in each subset. We recall the discussion in Section 3.4.1: any ‘inferred ’ variations to the parameters extracted by the analysis procedure must then arise from a combination of the effects of the individual window functions, and the varying inherent quality of the data (as reflected by the size of the white-noise component added in the time domain to each data string). The mean bias obtained from averaging over 30, 10-spectrum-set simulations is shown in Fig. 3, each datum being an average over adjacent 0≤ℓ≤3. The change-per-unit-rf coefficients extracted by the analysis (cf. Section 3.3) have been suitably multiplied to give the ‘implied ’ variation from solar minimum to maximum on the ordinate of the plot (under the pretence that our artificial data were collected at the different levels of activity reflected in the real data). The plotted error bars are the 1σ uncertainties in the mean values. (Here, we note that we have verified that the formal uncertainties in the raw heights and widths extracted by the mode-fitting procedure provide reliable estimates of the true precision in the data. As such, the error bars plotted in Figs 3 and 4 — which come from analyses of artificial data — and those in Figs 5 and 6— which come from analyses of real data — should be regarded as reasonable indicators of the true fluctuations in the displayed parameters.) The bias present in Fig. 3 can be understood as follows.

Consider first the large offsets at low frequencies. As indicated above, those spectra generated at epochs with high activity are also of generally poorer quality, owing to a combination of reduced fractional-fill and lower intrinsic data quality. Since the artificial data are designed to mimic the observed data, these trends will also be present in each artificial 10-spectrum subset. At low *n*, the S/N in the observed peaks is small — as a result, the fitting tends to underestimate the mode widths, since the wings of the resonant profiles are masked by background noise, and as such they take on a ‘spiky ’ appearance. Since this effect will be more marked in the poorer quality spectra — which we ‘pretend ’ have been collected at higher levels of solar activity — there is a downward, negative trend present in the width-variation plot at low *n*. Since, on average, underestimated widths imply overestimated heights (see discussion of anticorrelation in Section 3.2), there is an upward, positive trend present in the height plot. Over the central part of the p-mode spectrum, the bias in the inferred values is far less pronounced. Here, the competing effects of the pure window-background-induced bias and the varying intrinsic data quality largely offset each other (see Section 3.4).

Next, we tested the ability of the analysis procedure to uncover ‘programmed ’ changes to the artificial widths and heights. We constrained our choice of variation scenarios to test on the basis of results obtained from an initial analysis of the BiSON data, and also the diagnostic analyses made with the oscillator equations in Section 2.2. Here, we present results from the following: we imposed a minimum-to-maximum increase of 28 per cent at all *n* in the artificial widths [i.e., *δw*(ℓ,*n*)=0.28], and a corresponding decrease in all heights of 56 per cent. (Note that the changes were programmed as a linear function of the ‘activity ’ level for each spectrum.) This corresponds roughly to the variation uncovered from a first-cut analysis of the BiSON spectra. However, there is a suggestion in the BiSON results that the extracted variations are peaked at the centre of the 5-min envelope. In order to test the extent to which this might be due to bias in the analysis at low and high *n*— where the S/N is lower — we chose to test the flat, *n*-invariant change indicated above.

The mean results obtained from averaging over the resulting 30, 10-spectrum-set simulations are displayed in Fig. 4 (in the same manner as Fig. 3). The plotted error bars are the 1σ uncertainties in the mean values. The dotted lines indicate the programmed shifts. Over the central part of the spectrum, the analysis does a reasonable job of extracting the inserted changes; however, this is not so at low frequencies. This merely reflects the bias inherent in the procedure, as discussed above, and displayed in Fig. 3; again, it arises principally as a result of the poorer quality spectra ‘coinciding ’, by and large, with higher levels of activity.

We have corrected the raw, plotted data in Fig. 4 according to the inherent bias implied by the artificial data in Fig. 3— the corrected values have been plotted with triangular symbols in both panels. The corrected coefficients are little changed at frequencies greater than ∼2600 μHz. However, as noted above, they do nevertheless mimic the programmed variations reasonably satisfactorily. If we average the corrected shifts over the range 2600≤ν≤3600 μHz, we find inferred mean width and height variations of 25 and −51 per cent respectively (cf. the programmed changes of 28 and −56 per cent). The inferred width and height changes are therefore both slightly underestimated. However, the magnitude of the discrepancy is still only of the order of the uncertainty in the mean-shift estimate that would be expected from the analysis of a single, 10-spectrum subset (e.g., the actual BiSON set). We nevertheless note the presence of residual bias in the results.

The disagreement that remains at low *n* is still very marked. As such, this clearly indicates that our ability to draw reliable conclusions concerning changes over this part of the spectrum is limited. Since we will adopt a similar correction procedure for the BiSON data, the results to be uncovered should therefore be regarded with a little caution, in particular those at low *n*. We discuss these issues, in the context of the BiSON analysis, in the next section.

## Results: analysis of bison data

In Fig. 5 we show the results of analysing the 10 BiSON spectra. The upper two panels display the raw, absolute shifts — from solar minimum to maximum — of the natural logarithms of the mode linewidths (left-hand panel) and heights (right-hand panel), as a function of frequency. (Each plotted datum is an average over adjacent 0≤ℓ≤3.) The dotted line in each panel is a smooth fit through the appropriate Monte Carlo simulation results shown in Fig. 3; it serves to provide an estimate of the inherent bias present in the raw, extracted BiSON shifts plotted in the figure. The lower two panels show the same raw-shift data, but now averaged for each spectrum over 2600≤ν≤3600 μHz and plotted as a function of the mean 10.7-cm flux measured over the duration of each time series. All plotted error bars correspond to 1σ uncertainties.

We have corrected the raw shifts on the basis of the results from Fig. 3, and display the corrected values in Fig. 6. We also show the implied, corrected shifts for the modal velocity powers and energy supply rates; the corrections for υ^{2}(ℓ, *n*) and (ℓ, *n*) were obtained from suitable linear combinations of those established for *w*(ℓ, *n*) and *h*(ℓ, *n*). Again, the plotted errors are the associated 1σ uncertainties. In Table 1 we give the mean, implied variations (corrected) obtained by averaging over all modes in the central frequency range 2600≤ν≤3600 μHz. With reference to Fig. 3, the level of bias for these modes is modest; further, given the results in Fig. 4 (and associated discussion in Section 3.4.3), we should have reasonable confidence that the variations uncovered over this range should reflect the true changes taking place. These data have therefore been highlighted in Fig. 6 with the triangular symbols.

With reference to the results presented in Table 1, the extracted increase in the widths clearly implies a net increase in the damping rate, η, with activity. If we recall equation (6), the results for the energy supply rate are not inconsistent with the notion that the net forcing, *F*(ω), has remained constant over the solar cycle (given the precision in the data); remember that the energy supply rate is invariant of η. Further, if we accept the applicability of the oscillator equations in Section 2.1, then they indicate that the relative sizes and signs of the changes in Table 1 can arise from variations solely to η. Here, we note that indications of similar linewidth changes are suggested by the comparison of the size of the local dip at ≈3000 μHz returned by fits to VIRGO/LOI data collected at different epochs (Appourchaux 1998).

The issue of the possible presence of a frequency dependence in the uncovered changes is somewhat less clear. If we take for the moment the highlighted points in each of the panels in Fig. 6— which cover the central part of the p-mode spectrum — these data are suggestive of the variations being peaked near ∼3100 μHz. However, given the size of the uncertainties in the data, it is difficult to ascribe any satisfactory degree of significance to such a hypothesis. For example, the ratio of external and internal errors in the means of the shifts plotted in Fig. 6 are all close to unity, indicating that an appropriate model for each set of points is a simple, mean offset. If we now include the data outside this range, indications of ‘curvature ’ in the shifts are less pronounced.

## Conclusion and discussion

We have analysed 10 contiguous 8-month BiSON power spectra in an effort to search for changes to the excitation and damping of low-ℓ solar p modes over the falling phase of solar-activity cycle 22. For resonant modes in the range 2600≤ν≤3600 μHz (and 0≤ℓ≤3), our analysis uncovers a mean 24±3 per cent increase in the modal linewidths from activity minimum to maximum; a corresponding decrease of 22±3 per cent is observed in the modal velocity powers (mode amplitude squared). The net rate at which energy is supplied to the modes remains constant, at the level of precision of the observations (measured change 0±4 per cent). We use expressions derived from the equation of a damped harmonic oscillator to illustrate the diagnostic properties of the observables. These indicate that both the signs and relative sizes of the extracted variations could arise from changes solely to the damping; the net forcing of the modes need not change. Our results also hint at the changes possibly being maximal at frequencies near the centre of the 5-min envelope. If the true variations are indeed characterized in this manner, then they suggest an origin that is peaked in the outer, superadiabatic skin of the convection zone. The thermal adjustment time of this layer corresponds to a temporal frequency of approximately 3000 μHz (Balmforth 1992), and as such the eigenfunctions of those modes with similar resonant frequencies couple strongly to it. This also gives rise to a local minimum in the damping rates in this part of the spectrum.

The direct impact — via the Lorentz force — of a horizontal magnetic field structure on the predominantly vertically oriented motions of the modes at the solar surface has been addressed by, amongst others, Gough & Thompson (1988), Goldreich, Murray & Willette (1991) and Hindman, Jain & Zweibel (1997). Hindman et al. considered the suppression of p-mode power by a uniform field with an rms size of a few hundred gauss, and found reductions of up to 50 per cent in the presence of a field of rms strength ∼600 G.

Goldreich et al. (1991) measured the mean, global radial rms field strength from the central portion of whole-disc Kitt Peak magnetograms, and found that this varied from ≈120 to ≈190 G between solar minimum and maximum.^{2}

Any rough calculation based upon these values must be approached with some caution. This is because the magnetograms only measure the line-of-sight component of the surface field, not the horizontal component we require here; assumptions must therefore be made regarding the relative sizes of the vertical and horizontal fields. If we assume that any horizontal field is of similar size to the vertical component determined by Goldreich and colleagues — and this is very questionable — both the absolute sizes, and changes in rms strength between solar minimum and maximum, appear to be insufficient to explain the changes found in the BiSON data, i.e., the direct modification of the p-mode eigenfunctions by a magnetic field probably cannot alone account for the observations. We also note that Hindman, Jain & Zweibel's analyses infer maximal power suppression at frequencies of ≈4000 μHz; this does not appear to be matched by the BiSON results (although we again note that the issue of frequency dependence is uncertain owing to the size of the uncertainties on the determined shifts).

The Wilson depression of larger scale magnetic structures (vertical cross-section greater than ≈600 km; see, e.g., Stix 1991) gives rise to a localized patch over the area in question where the effective observation height in a spectral line is lowered by as much as a few-hundred kilometers (Solanki, Walther & Livingston 1993). In order to explain our uncovered reduction in p-mode power of ∼22 per cent, a disc-averaged observation-height depression difference of roughly 50 km would be required between the extremes of the activity cycle. Because of the low filling factor of strong, active regions (e.g. Montesinos & Jordan 1993) the depression change expected — of perhaps a kilometre at most — is inadequate to explain the observed power suppression; further, the action of this phenomenon does not impact on the damping of the modes.

Another possible explanation is based on variations in the convective properties near the solar surface which arise from the influence of magnetic structures. Balmforth, Gough & Merryfield (1996) concluded that purely thermal structural changes in the convection zone could not account for observed changes to the modal eigenfrequencies and the solar luminosity over the activity cycle. Houdek & Gough (1999) have expanded on the work of Balmforth et al., and shown that changing the horizontal length-scale of the convective eddies affects the p-mode damping rates. They note that this is an intriguing result, given reports of solar-cycle variations in the horizontal size of the granules (e.g. Berrilli et al. 1999). Clearly, this work needs to be tied in with the observed changes in the BiSON data (Houdek et al., in preparation).

The results presented in this paper provide another constraint on the theory; clearly, it would be desirable to take into account, in some form, the impact of the surface magnetic morphology on the interactive processes occurring between the acoustic eigenmodes and the convection. Since a cogent theory for turbulence does not yet exist, this is clearly not a rudimentary problem.

### Acknowledgments

We thank all those who are — or have been — associated with BiSON, in particular H. K. Williams, J. Litherland and J. Allison in Birmingham, and P. Fourie at SAAO; and our hosts, R. Stobie (SAAO); the Carnegie Institution of Washington; the Australia Telescope National Facility (CSIRO); E. J. Rhodes (Mt. Wilson, California), and members (past and present) of the IAC, Tenerife. WJC thanks T. Appourchaux, D. Gough and G. Houdek for stimulating discussions, and M. Gabriel for his constructive remarks as referee. BiSON is funded by the UK Particle Physics and Astronomy Research Council. WJC acknowledges the support of an ESA Internal Fellowship.

## References

^{3}G (e.g. Wiehr 1979; Stenflo 1989), or the surface-averaged value of the field. The latter measure is smaller than the rms field, because flux breaks through the surface in spatially thin concentrations.