A study of the rapid rotator ζ Aql: differential surface rotation? Open Access

ABSTRACT

We report new, extremely precise photopolarimetry of the rapidly-rotating A0 main-sequence star ζ Aql, covering the wavelength range ∼400–900 nm, which reveals a rotationally-induced signal. We model the polarimetry, together with the flux distribution and line profiles, in the framework of Roche geometry with ω-model gravity darkening, to establish the stellar parameters. An additional constraint is provided by TESS photometry, which shows variability with a period, P_phot, of 11.1 h. Modelling based on solid-body surface rotation gives rotation periods, P_rot, that are in only marginal agreement with this value. We compute new ester stellar-structure models to predict horizontal surface-velocity fields, which depart from solid-body rotation at only the ∼2 per cent level (consistent with a reasonably strong empirical upper limit on differential rotation derived from the line-profile analysis). These models bring the equatorial rotation period, P_rot(e), into agreement with P_phot, without requiring any ‘fine tuning’ (for the Gaia parallax). We confirm that surface abundances are significantly subsolar ([M/H] ≃ −0.5). The star’s basic parameters are established with reasonably good precision: |$M = 2.53\pm 0.16\, \mbox{M}_{\odot }$|⁠, log (L/L_⊙) = 1.72± 0.02, |$R_{\rm p}= 2.21\pm 0.02\, \mbox{R}_{\odot }$|⁠, T_eff = 9693 ± 50 K, |$i = 85{^{+5}_{-7}}^\circ$|⁠, and ω_e/ω_c = 0.95 ± 0.02. Comparison with single-star solar-abundance stellar-evolution models incorporating rotational effects shows excellent agreement (but somewhat poorer agreement for models at [M/H] ≃ −0.4).

1 INTRODUCTION

The discovery of rotationally-induced wavelength-dependent linear polarization in Regulus (α Leo; Cotton et al. 2017a) unveiled a new tool for investigating the properties of rapid rotators, opening the possibility of reasonably precise determinations of mass (and other parameters) in single stars. Subsequent studies have shown that high-quality photopolarimetry can afford powerful tests of stellar-atmosphere physics and of evolutionary models incorporating rotation (Bailey et al. 2020b; Lewis et al. 2022).

Here, we examine new results for ζ Aql (HD 177724, HR 7235, ‘Okab’), to explore further the diagnostic potential of very precise photopolarimetry. Selected basic data for the star are assembled in Table 1, including some of the observational material reviewed in Section 2. Section 3 outlines our modelling procedures. Results are given in Section 4, which examines the question of possible differential surface rotation. Section 5 confronts the inferred stellar parameters with stellar-evolution models. The determination of v_esin i, and other aspects of the line-profile analysis, are relegated to Appendix  B.

2 OBSERVATIONS

2.1 Photopolarimetry

The new observations reported here were obtained using HIPPI, the HIgh-Precision Polarimetric Instrument, and its successor, HIPPI-2; the design and operation of these instruments are described by Bailey et al. (2015, 2020a). We also made use of one previously published observation (Bailey, Lucas & Hough 2010), obtained using PlanetPol (Hough et al. 2006).

Details of observing runs and instrument configurations are given in Appendix  A (Table A1), and the individual observations of ζ Aql are listed in Table 2. As well as shot noise, the uncertainties reported therein on the normalized Stokes parameters q, u, and on the polarization p, include (in quadrature) a wavelength-dependent positioning error which results from inhomogeneities across the face of the ferroelectric liquid-crystal modulator, and which sets the accuracy limit for the instruments. For HIPPI-2, this limit ranges from 1.1 parts per million (ppm) for the reddest passbands, through 2.5 ppm in g′, to 13.7 ppm in the 425SP filter – similar to, but slightly better than, the corresponding figures for HIPPI (Bailey et al. 2020a). Position-angle (PA) calibration was performed using highly-polarized standard stars; typical PA zero-point uncertainties are ∼1 ^○ (i.e. are smaller than the statistical errors on our measurements).

Results are plotted in Fig. 1; the scatter in the observations slightly exceeds the quoted formal errors (though is still small compared to results from traditional polarimeters), as was also found in our study of θ Sco (Lewis et al. 2022). Although intrinsic low-level polarization variability cannot be ruled out, the scatter is consistent with a combination of instrument-configuration changes and imperfect characterization of low-polarization standard stars (used to correct for polarization arising in the telescope optics), which can lead to zero-point drifts of up to ∼10 ppm between runs (Bailey et al. 2021). These factors are accommodated in our modelling by use of bootstrapping in the error analysis (Section 4).

2.2 TESS photometry

The TESS satellite (Ricker et al. 2015) observed ζ Aql in sector 54 (2022 July–August). The light-curve is shown in Fig. 2; a periodic signal is immediately evident. A generalized Lomb–Scargle periodogram ( Ferraz-Mello 1981; Zechmeister & Kürster 2009) yields a fundamental frequency and corresponding semi-amplitude of

|$\phantom{X} \nu _0= 2.1587 (8)$| d⁻¹  [P_phot = 11.118(4) h]

|$\phantom{X} a_0=56.5 (25)$| ppm

where parenthesized values are 1σ uncertainties on the least significant digits, generated from Monte Carlo simulations using a residual-permutation, or ‘prayer beads’, algorithm (adopted because the residuals are strongly correlated as a result of lower-frequency drifting). The signal appears to be rather simple; power at the first harmonic accounts almost entirely for departures from a pure sinusoid (Fig.   2).

2.3 Ancillary observational data

Some additional observational material is required for our analysis; furthermore, given the level of precision of the polarimetry, consideration needs to be given to potential contaminating sources (whether stellar or circumstellar).

2.3.1 Parallax

Parallaxes, required principally to establish the stellar radius, are available from both the Hipparcos and Gaia missions (van Leeuwen 2007; EDR3, Gaia Collaboration 2021a, b). The Hipparcos result is 1.05 ± 0.38 mas larger than the EDR3 value (Table 1), a 2.6σ difference. Although this difference is of little consequence for most derived parameters (cf. Table 4), it proves to be of some significance for the interpretation of the photometric period (Section 4.1). We therefore performed calculations for both values (taking the distance, ∼26 pc, to be simply the inverse of the parallax¹).

2.3.2 Interstellar polarization and extinction

Cotton et al. (2017a) developed a model of interstellar polarization from which we expect p(λ_max) ≃ 30 ppm for ζ Aql, where λ_max is the wavelength of maximum linear polarization (typically ∼0.5 |$\mu$|m). Observations of stars over a range of distances in the general direction of ζ Aql  support this estimate, with an upper limit of ∼50–60 ppm (Fig.   3). The interstellar polarization is directly estimated as part of the modelling (Section 3.1), and yields p(λ_max) ∼17–24 ppm (in position angle θ_i = 55 ± 2^○), in good accord with the Cotton et al. model.

Such small polarizations imply very little foreground dust and interstellar reddening. Serkowski, Mathewson & Ford (1975) found E(B − V) ≳ p/9 per cent for nearby stars, suggesting a barely non-zero extinction. For the purposes of correcting the observed flux distribution we adopt |$E(B-V) = 0{${_{.}^{\rm m}}}005$| as a suitably small, if arbitrary, round-number estimate; our results are very insensitive to the exact value (Table 4).

2.3.3 Projected rotation velocity

There is, additionally, some sensitivity to the surface-rotation profile (i.e. the variation, or otherwise, of ω with colatitude θ). For example, models generated with the ester stellar-structure code, discussed in Appendix B5, have a differential rotation profile that results in v_esin i values ∼4 km s⁻¹ smaller than does solid-body surface rotation.

Our modelling takes these factors fully into account; the overall range of acceptable v_esin i values is ∼300–325 km s⁻¹, with v_esin i = 306 km s⁻¹ for our final preferred model (Section 4).

2.3.4 Companion stars

The Washington Double Star Catalog (WDS; Mason et al. 2001) lists four visual companions to ζ Aql A (=WDS J19054+1352A); of these, only the B component, at separation ρ = 7″ (De Rosa et al. 2014; Gaia Collaboration 2021a), is close enough to potentially affect the observations discussed in this paper. (It is also the only physical companion, according to Gaia astrometry.)

The B component was discovered by Burnham (1874), who described it as ‘not fainter than...11 mag’. Wallenquist (1947) reported a visual² magnitude difference of |$\Delta {v} = 8{${_{.}^{\rm m}}}45$|⁠, while |$\Delta {G}=7{${_{.}^{\rm m}}}85$| (Gaia Collaboration 2021a) and |$\Delta {K} = 4{${_{.}^{\rm m}}}87$| (De Rosa et al. 2014). The colours and absolute magnitudes are consistent with an early-M dwarf companion, which would contribute <1 per cent of the flux at λ<1|$\mu$|m (<0.1 per cent at λ<0.6 |$\mu$|m), rising to ∼3 per cent only for λ ≳ 4 |$\mu$|m. The B component is therefore of no importance for the observations and analysis reported here.

As pointed out to us by our referee, the Gaia image parameters can provide additional information on potential close companions. The Renormalised Unit Weight Error (RUWE) for the ζ Aql A astrometric solution is 2.5, which initially appears to be rather large compared to the value of ∼1 expected for well-behaved solutions of single stars. However, we find that this is value is actually typical of very bright stars; the 548 stars in DR3 with G ≤ 4.0 have a median RUWE of 2.73. The ipd_gof_harmonic_amplitude parameter is a measure of image asymmetry; its value of 0.09 is fully consistent with a circular image (Fabricius et al. 2021). Peterson et al. (2006) also imply that companions with ΔR ≲ 8 within |$0{_{.}^{\prime\prime}} 5$| of ζ Aql A are ruled out by their interferometric observations.

2.3.5 Is ζ Aql A a spectroscopic binary?

At the time of writing, the WDS carries an unattributed note that ‘A is a spectroscopic binary’. We have been unable to find any documented source of that report. We therefore examined the sixty-three good-quality high-resolution spectra discussed in Appendix  B, obtained between 2005 May and 2014 June, which sample time-scales of minutes, hours, days, and years. Simple visual inspection reveals no obvious line-profile or radial-velocity variations, and cross-correlation velocity measurements of the Ca ii K line yield an r.m.s. dispersion of only 3.6 km s⁻¹ (cp. the resolution element, ∼4.6 km s⁻¹, and v_esin i, ∼300 km s⁻¹). We proceed on the assumption that if ζ Aql A is indeed a spectroscopic binary, then that is of no consequence for our analysis.

2.3.6 Exozodiacal-dust emission

Absil et al. (2008) reported a K-band excess of (1.69 ± 0.27) per cent of the photospheric flux,³ based on differences between short-baseline interferometric visibilities observed with CHARA and those predicted from a colour-based surface-brightness estimate of the angular diameter. Using a different detector (though the same instrument and methodology), Nuñez et al. (2017) obtained a ΔK of (1.23 ± 0.38) per cent, attributing this excess to hot (∼1000 K) exozodiacal dust.

The flux contribution from dust emission of this nature is again too small to have any direct consequences for our study. In principle, if aligned grains were involved, they might influence the photopolarimetry; however, current evidence suggests that the grains responsible for exozodiacal emission are too small to produce polarization at optical wavelengths (Marshall et al. 2016).

2.3.7 Interferometry

Finally, long-baseline optical interferometry can provide a useful check on our results. We found two reports in the literature: Peterson et al. (2006) briefly summarize an otherwise unpublished detailed analysis of NPOI measurements, while Boyajian et al. (2012) give a mean angular diameter from CHARA data. Results are included in Table 1; the angular diameters from the two studies differ by ∼7 per cent (∼3.5σ). Since the NPOI angular diameter (⁠|$\sim \, VR$| passband) is larger than the CHARA value (⁠|$\sim K$|⁠), the discrepancy cannot be attributed to the possible exozodiacal-dust emission discussed in Section 2.3.6.

3 MODELLING

The principal motivation for our modelling is to determine values for basic stellar parameters from the photopolarimetry. We conducted our analysis in the framework of standard Roche geometry (e.g. Collins 1963), with ω-model gravity darkening (Espinosa Lara & Rieutord 2011).

While the photometric variability suggests the possibility of some departure from axial surface-brightness symmetry, it is at a very low level (and in any case, we have no way of characterizing it in an appropriate manner). Because our polarimetric observations were taken at arbitrary phases, we do not anticipate systematic effects, and the bootstrap error analysis accommodates any increase in observational scatter that may arise.

3.1 Overview

Photospheric polarization was computed using the code described by Cotton et al. (2017a) and Bailey et al. (2020b), which calculates local surface intensities using the synspec spectral-synthesis program (Hubeny, Stefl & Harmanec 1985; Hubeny 2012), modified for fully polarized radiative transfer using the vlidort package (Spurr 2006). The underpinning atmosphere models are custom atlas9 line-blanketed LTE calculations (Castelli & Kurucz 2003).

The model polarization depends principally on four quantities (or equivalent surrogates): a reference temperature and gravity (e.g. polar values T_p, g_p); the inclination of the rotation axis to the line of sight, i; and the rotation parameter, ω_e/ω_c.

We can reduce this four-dimensional parameter dependency to a two-dimensional grid by exploiting complementary observations which independently constrain the temperature and gravity (Section 3.2), allowing us to compute predicted polarizations as functions of i and ω_e/ω_c alone (at the self-consistent T_p, g_p values).

3.2 Reducing the 4D dependency

3.2.1 Temperature

The method underpinning our temperature determinations is closely akin to the Infra-Red Flux Method of Blackwell & Shallis (1977), the basic principle being that the ratio of the observed fluxes in two suitable regions is a measure of temperature. Ideally, ‘suitable’ regions should have strongly different temperature dependences (e.g. be on either side of the peak of the flux distribution), and should record a significant part of the total luminosity.

We used V photometry and the UV flux from observations made with the International Ultraviolet Explorer (IUE), as summarized in Table 1. We investigated the use of longer-wavelength (RIJ) photometry in place of V, but found that this did not afford any useful gain in sensitivity (in part because of increasing observational uncertainties).

We used all available archival IUE observations obtained through its spectrographs’ large apertures at low resolution (resolving power R ≃ 350), which provide the most reliable flux measurements. The nine short-wavelength (∼115–198 nm) and nine long-wavelength (∼185–335 nm) flux-calibrated spectra were combined with weights proportional to exposure times, excluding image LWR 7295, which has a notably poor signal:noise ratio. All fluxes were corrected for the adopted interstellar extinction using a Seaton (1979) curve, with A_V/E(B − V) = 3.1. The UV flux accounts for ∼25–30 per cent of the luminosity, and the V band for ∼10 per cent (cf. Fig. 8).

3.2.2 Metallicity, radius

Conversion of the observed flux ratio to a temperature is achieved by means of model-atmosphere flux distributions. This introduces a metallicity dependence, principally through line blanketing, whereby the observed UV flux is reproduced by lower-temperature models at lower metallicity. Gray et al. (2003) and Wu et al. (2011) report [M/H] = −0.68 ± 0.09, −0.52 ± 0.16, respectively, for ζ Aql; the synthetic spectra described in Appendix  B also indicate substantially subsolar metallicity (cf. Fig 6). Our analysis is not sensitive to the precise value (Table 4); we adopt [M/H] = −0.5 as a suitable round-number value (for calculation of both parameter grids and model polarizations). This results in effective temperatures ∼2 per cent lower than would be inferred from solar-abundance models.⁴

With the temperature established, the observed flux level directly yields the angular diameter, which can be converted into a stellar radius given the distance.

3.2.3 Rotational effects

Significant rotation introduces dependences on ω_e/ω_c and i to the modelled fluxes, through gravity darkening and aspect effects. However, for a given (or assumed) value of i, the equatorial rotation velocity, v_e, follows directly from the observed projected rotation velocity, v_esin i, thereby establishing ω_e (from the equatorial radius). For a given (or assumed) value of ω_e/ω_c, the corresponding mass can then be inferred (from equation 1). Consequently, for any specified i, ω_e/ω_c combination, there is a unique pair of temperature and polar-gravity values (and associated mass and radius) that reproduce the observed fluxes.

To determine those values in practice, we run a series of models to calculate fluxes at 100-K steps in T_eff, for specified i, ω_e/ω_c pairs [at given v_esin i, d, [M/H], and E(B − V)]. These models are computed in full limb- and gravity-darkened Roche geometry, using the code described by Howarth (2016), with atlas9 intensities from Howarth (2011). At each T_eff, we determine the polar radius that reproduces the observed flux (UV or V) using a simple interval-halving algorithm, in order to establish a locus of acceptable values in T_eff, log (g_p) space. The required final result is given by the intersection of the UV and V loci (which is well defined for this T_eff regime, confirming that these passbands are ‘suitable’ in the sense discussed in Section 3.2.1).

4 RESULTS

Fig. 4 shows selected results of the T_eff/log (g_p) modelling described in Section 3.2. Our first parameter grid was constructed with models based on the Hipparcos parallax (as it is more precise than the Gaia value) and the ‘Occam’s razor’ assumption of solid-body surface rotation (consistent with the quite strong upper limit on differential rotation obtained from the line-profile analysis reported in Section B3), using equation (B2a) to link v_esin i and ω_e/ω_c.

Detailed polarization models were constructed for each parameter-grid point, as set out in Section 3.1, with results passband-integrated for comparison with the observations. The resulting χ² map is shown in Fig. 5. Confidence intervals on this map were generated from bootstrapped data sets by evaluating best-fit inclinations for each ω_e/ω_c from the χ² maps, for 1000 samplings. The 1σ ranges give the corresponding Δχ², from which the confidence intervals may be inferred.

4.1 The rotation-period ‘problem’...

Table 3 shows that, for the initial parameter grid, the surface rotation periods, P_rot, along the χ² ‘valley’ of Fig. 5 are in the range 10.2–10.7 h (for i = 60 → 90^○). These values are sufficiently close to the TESS photometric period (P_phot = 11.1 h; Section 2.2) to suggest that the photometric period may very well be the rotation period.

In part to see if the model P_rot values could be straightforwardly reconciled with P_phot, we conducted parameter-sensitivity tests, and also recalculated the parameter grid with the Gaia parallax (which, being smaller, leads to slightly larger radii, hence larger values of P_rot for given v_e). Selected results are included in Tables 3 (Gaia parallax, columns 11, 12) and 4 (sensitivity tests).

The inferred effective temperature is, unsurprisingly, mildly sensitive (at the ∼100 K level) to [M/H], and to the value of the integrated UV flux, but other parameters (including P_rot) appear to be quite robustly determined, with changes that are generally within the spread of values resulting from the i, ω_e/ω_c indeterminancy. This robustness propagates into the polarization modelling; a solar-abundance test grid (with concomitant changes in all basic parameters) gives a χ² map that is practically indistinguishable from that shown in Fig. 5, confirming our expectation that the modelled polarization is primarily sensitive to i and to ω_e/ω_c, with relatively little dependence on other parameters (within reasonable bounds).

Although confirming the obvious result that smaller parallax (greater distance, hence larger radius) and smaller v_esin i should push P_rot to larger values, the model rotation periods in Tables 3 and 4 still all fall short of P_phot. The differences are smallest for the highest-inclination, Gaia-parallax models, with the extreme i = 90^○ model in Table 3 requiring a reduction of only ∼5 km s⁻¹ in v_esin i to force agreement. However, Fig. B1(f) indicates that even as small a change as this is barely compatible with the line-profile analysis (for a continuum placement chosen to give consistency with solid-body surface rotation, which already results in lower v_esin i values than would otherwise be found; Fig. B3).

We conclude that the combined hypotheses of both (i) solid-body surface rotation and (ii) P_phot = P_rot, while not completely ruled out, are at best only marginally consistent with the data.

4.2 ...and some possible resolutions

There are several possible circumstances that could address this modest discrepancy between P_rot and P_phot. For example, reducing the Gaia parallax by ∼2σ would increase the distance, model radius, and hence P_rot, by a further ∼2 per cent (although the parallax would then be ∼11σ from the Hipparcos value). Such ‘fine tuning’ cannot be excluded; nevertheless, we should also not discount the possibility that P_phot need not necessarily be identical to the solid-body rotation period.

We note, for example, that both g- and r-mode pulsation periods can be in the same range as the rotation period. The attraction of pulsation as the origin of photometric variability is that it does not require time-variable magnetic fields to be invoked as a mechanism to generate rotational modulation through starspots (in a broad sense). Where magnetic fields have been directly detected in OBA stars (which are thought to have primarily radiative envelopes), they appear to be fossil remnants, an interpretation consistent with associated highly reproducible, strictly periodic photometric and spectroscopic variability (e.g. Hubrig & Schöller 2021). This behaviour contrasts with much of the low-amplitude variability revealed by space photometry and widely attributed to rotational modulation (e.g. Balona & Abedigamba 2016).

Saio et al. (2018) found a low-frequency ‘hump’ in time-series power spectra at frequencies just below the rotation frequency, resulting from r-mode pulsation with azimuthal wavenumber |m| = 1. Subsequently, Lee & Saio (2020) and Lee (2021, 2022) have shown that overstable convective modes in the core of an early-type star can couple with g modes in the radiative envelope, provided the core rotates slightly faster than the envelope. These modes are therefore candidates for the processes underlying the photometric variability.

However, the r-mode modelling suggests a more complex frequency spectrum than is exhibited by ζ Aql A, while the periods driven by overstable convective modes should necessarily be shorter than the surface rotation periods (whereas, we find to the contrary, that P_phot ≳ P_rot). While pulsational variability close to the rotation frequency remains a possibility, arguments in favour of that hypothesis do not appear compelling at present.

4.3 Differential rotation?

If, instead, the photometric variability is attributed to rotational modulation of surface-brightness inhomogeneities then differential surface rotation offers a straightforward solution to the discrepancy between the modelled (equatorial) rotation period, P_rot(e), and P_phot: spots could simply be at latitudes rotating with a different angular speed.

As described in Appendix B2, an ad hoc parametrization of differential rotation was considered as part of our initial examination of rotational velocities. That analysis of the line-profile shape found no direct evidence for substantial differential rotation. However, a purely empirical approach cannot rule out lower-level differential rotation, and necessarily introduces an arbitrary element (i.e., an ad hoc characterization of the dependence of ω on colatitude θ).

For a second phase of the analysis, we therefore pursued a direct physical approach, using surface-rotation profiles, ω(θ), generated from ester 2D stellar-structure models, as elaborated in Appendix B5 (and illustrated in Fig. B2).

We first repeated the exercise of establishing the relationship between v_esin i and ω_e/ω_c for these differentially-rotating models (cf. Appendices   B4,B5), then generated new parameter grids (for both Hipparcos and Gaia parallaxes). Selected results are incorporated into Figs 4 and 5. Fig. 5 shows that only the combination of Gaia parallax and ester rotation profiles gives agreement between P_rot(e) and P_phot at i, ω_e/ω_c values that are consistent with the polarimetry, without requiring any fine tuning.

This result arises because ester-type differential rotation has greatest angular velocity at temperate latitudes (Fig. B2), leading to projected equatorial rotation velocities that are generally ∼1–2 per cent smaller than rigid-rotator results (Appendix B5). The relationships between v_esin i and ω_e/ω_c are essentially independent of parallax, but with lower v_e values favoured by lower ω_e/ω_c and by higher inclinations. The combination of these effects, along with their non-linear dependences on ω_e/ω_c, means that, within the χ² valley of Fig. 5, P_rot(e) matches P_phot only for Gaia+ester models, at high inclinations (i ≳ 80^○, such that v_e ≃ v_esin i) for the ω_e/ω_c values indicated by the photopolarimetry.

Given the sensitivity of P_rot(e) to rather small changes in v_esin i, this is clearly not a unique result (nor a strong validation) of the ester rotation profile. For example, an arbitrary surface-rotation law of the form of equation B1, with α ∼ −0.03, would give a similar outcome. Nevertheless, the success of the (almost) ‘no free parameters’ ester models in bringing about agreement between P_rot(e) and P_phot is encouraging and suggestive.

A further caveat is that Fig. 5 is not fully self-consistent, as the various P_phot = P_rot(e) lines are calculated under slightly different sets of assumptions, while the χ² map is specifically for solid-body surface rotation and the Hipparcos parallax. We have not addressed this minor inconsistency, for several reasons:

• We consider the existing calculations already to be sufficient to demonstrate that the modelled P_rot(e) and observed P_phot can readily be reconciled; while other stellar parameters are insensitive to details of the rotation profile and parallax.

• We know the χ² map is very insensitive to most input parameters (including small changes in parallax), excepting the i, ω_e/ω_c pair. While any sensitivity to differential surface rotation is less clear, it is unlikely to be large, given the small departures from solid-body rotation implied by the ester profiles (and the empirical limits on strongly differential rotation).

• If the photometric variability does arise from starspots, we do not know their latitude(s). Although we have focused on reconciling the equatorial rotation period of the models with P_phot, the spot rotation periods could easily be 1–2 per cent different (in either direction) in the presence of differential surface rotation. As shown in Fig. 5, such small differences can have a relatively large effect on the location of the P_rot(e) locus in the i, ω_e/ω_c plane, which is likely to dominate the uncertainties.

5 DISCUSSION

Selected numerical results for the Gaia+ester parameter grid are given in Table 5 (other parameter grids give results intermediate between those in Tables 3 and 5). We take the ‘base’ model listed there, for which P_rot(e) = P_phot, as our adopted specific solution. If starspots do give rise to the photometric variability, then the combination of high axial inclination and ∼continuous variation implies that they must have an extensive distribution in longitude.

The line-profile modelling for the adopted solution is shown in Fig. 6, the predicted and observed polarizations in Fig. 7, and the flux distributions in Fig. 8. All these comparisons show satisfactory agreement between models and observations.

The mean angular diameters predicted by the tabulated models are in excellent accord with the interferometric value given by Boyajian et al. (2012; 0.895 ± 0.017 mas), but are inconsistent with the interim analysis reported by Peterson et al. (2006; cf. Table 1); their values of i = 90^○, ω_e/ω_c = 0.99, θ_* = 45^○ are also at odds with the photopolarimetric results (Table 3).⁵

5.1 Comparison with evolutionary models

We compare our empirical results with models of the evolution of rotating stars from Georgy et al. (2013) in Fig. 9.⁶ Their grids are for a range of ZAMS rotation rates, ω_e/ω_c(0), and include metallicities representative of solar and LMC abundances (Z = 0.006, corresponding to [M/H] ≃ −0.4). All our empirically determined masses fall within the range |$M=2.53^{+0.20}_{-0.13}\,\mbox{M}_{\odot }$|⁠, in excellent agreement with the evolutionary mass for the solar-abundance tracks. The LMC-abundance tracks suggest evolutionary masses ∼0.3 M_⊙ lower, barely consistent with empirical values.

Main-sequence A-type stars show a range of surface-abundance anomalies, usually involving selective metal enhancements (the Am, Ap, and HgMn stars); only stars in the λ Boo class are noted for their metal depletions. This group is also characterized by relatively rapid rotation. While ζ Aql is not a classic λ Boo star in terms of its spectral morphology (Gray, personal communication), its subsolar metallicity may arise through a similar mechanism, generally thought to involve photospheric accretion of depleted gas (e.g. Venn & Lambert 1990; Jermyn & Kama 2018). In that case, we would expect solar abundances to be more relevant to its evolution, as found in these comparisons

At these masses, the evolutionary tracks are not strongly sensitive to the precise value of ω_e/ω_c(0), although a high value is, of course, required for ζ Aql A. A ZAMS value close to ∼0.95 is consistent with observations (Fig. 9, right-hand panel).

6 SUMMARY AND CONCLUSIONS

We have presented new, very precise photopolarimetry of ζ Aql (Table   2). Modelling those observations, together with supplementary analyses of the flux distribution and rotational velocity, allows us to determine the locus of allowed combinations of i and ω_e/ω_c (Fig. 5). The polarimetry alone cannot break the degeneracy between these two parameters, but limits their values to i ≳ 60^○, ω_e/ω_c ≳ 0.93.

Periodic photometric variability, demonstrated here for the first time (from TESS observations), provides additional constraints under the plausible assumption that the newly established photometric period, P_phot = 11.12 h, can be identified with the rotation period. The rotation periods of models based on rigid-body surface rotation are only marginally consistent with P_phot, requiring extreme values and fine tuning of parameters to push v_e and/or the parallax to appropriately low values. However, model equatorial rotation periods are found to be in good agreement with P_phot for the combination of Gaia parallax and the differential surface rotation predicted by ester models.

The inferred physical parameters of ζ Aql A are quite insensitive to these issues, as demonstrated by the small range of solutions listed in Tables 3–5; our adopted specific characterization is given in column 3 of Table 5. Taking the full ranges of parameter values in that Table as a reasonably conservative estimate of the 1σ uncertainties, we find |$M = 2.53\pm 0.16\, \mbox{M}_{\odot }$|⁠, log (L/L_⊙) = 1.72± 0.02, |$R_{\rm p}= 2.21\pm 0.02\, \mbox{R}_{\odot }$|⁠, T_eff = 9693 ± 50 K, |$i = 85{^{+5}_{-7}}^\circ$|⁠, and ω_e/ω_c = 0.95 ± 0.02.

Comparison of our results with grids of single-star evolution calculations shows excellent agreement for solar-abundance models, but poorer agreement with models at lower metallicities that approximately match the depleted surface abundances. This suggests that the observed photospheric depletions may not be global, but instead confined only to the surface layers.

ACKNOWLEDGEMENTS

This paper is based in large part on data obtained with the Anglo-Australian Telescope at Siding Spring Observatory; we acknowledge the traditional owners of the land on which the AAT stands, the Gamilaraay people, and we pay our respects to elders past and present. We made use of the Washington Double Star Catalog, maintained at the U.S. Naval Observatory, as well as observations made with the International Ultraviolet Explorer and TESS satellites, obtained from the MAST data archive at the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS 5-26555. Funding for the TESS mission is provided by the NASA Explorer Program. CFHT data were accessed by using the facilities of the Canadian Astronomy Data Centre, operated by the National Research Council of Canada with the support of the Canadian Space Agency. We also benefitted from NASA’s Astrophysics Data System bibliographic service, and the SIMBAD data base, operated at CDS, Strasbourg, France. We thank Nicholas Borsato, Dag Evensberget, Behrooz Karamiqucham, Jonathan Marshall, and Jinglin Zhao for contributions to observing runs, our anonymous referee for useful remarks, and Conny Aerts, Derek Buzasi, Richard Gray, and Michel Rieutord for helpful correspondence. DVC thanks the Friends of MIRA for their support.

DATA AVAILABILITY

The new polarization data used for this project are listed in Table 2. All other data are from publicly accessible archives.

Footnotes

REFERENCES

APPENDIX A: SUMMARY OF OBSERVING RUNS

Technical details of the observing runs leading to the results given in Table 2 are summarized in Table A1. HIPPI and HIPPI-2 are dual-beam photopolarimeters that use ferroelectric liquid crystals (FLCs) for primary modulation at 500 Hz, in order to overcome seeing noise and thereby to achieve high precision. The FLC used for all HIPPI/-2 observations of ζ Aql was manufactured by Boulder Nonlinear Systems; the performance of this unit has evolved over time, an issue addressed by the reduction pipeline used to process all the data from both instruments (Bailey et al. 2020a).

We made use of two types of Hamamatsu photomultiplier tube (PMT) as detectors; the blue-sensitive H10720-210, which we denote ‘B’, and red-sensitive H10720-20, ‘R’. Six broadband filters were employed: custom-built 425-nm and 500-nm ‘short-pass’ and 650-nm ‘long-pass’ filters (425SP, 500SP, and 650LP, respectively; Bailey et al. 2020a), SDSS g′ and r′, and Johnson V. The r′ filter was paired with both the R and B detectors; the 650LP filter only with R; and the remainder only with B.

Telescope optics introduce a small wavelength-dependent polarization, which was removed by reference to observations of low-polarization standard stars (cf. Table A1).

APPENDIX B: ROTATIONAL VELOCITY

The equatorial rotation velocity provides an important constraint on the stellar mass. A summary of published estimates of v_esin i for ζ Aql A is given in Table B1; Abt & Morrell (1995) provide the only primary measurement on the Slettebak et al. (1975) system that we have been able to locate. They note the absence of reliable fast-rotating calibrators in their data set, and mark their measurement as uncertain. We therefore undertook a new analysis in order to determine a modern precise value for v_esin i, using the Fourier-transform method (e.g. Carroll 1933; Smith & Gray 1976).

B1 Data

A search of online archives showed that observations obtained with the ESPaDOnS echelle spectropolarimeter at the Canada–France–Hawaii Telescope (CFHT) are of particularly good quality, and numerous. After rejecting a handful of relatively low-quality exposures, we corrected relevant sections of each of the remaining sixty-three spectra for weak telluric absorption (dividing individual spectra by a scaled telluric template constructed from the data merged in topocentric velocity space), then merged them (after correcting to heliocentric velocities). The resulting noise-weighted mean spectrum has a continuum signal:noise ratio of ≳3000 at ∼635 nm (as measured from residuals to low-order polynomial fits), at a resolving power of R ≃ 65 000. We identified no obvious line-profile or radial-velocity variability in these data (Section 2.3.5).

B2 Modelling: ad hoc characterization of differential rotation

Latitudinal differential rotation introduces changes to line-profile shapes, principally by modifying the Doppler redistribution of absorption arising at temperate latitudes (for given v_e and sin i; changes in gravity darkening introduce further, but secondary, effects). This in turn affects the Fourier transform of the profiles – notably, the separation of the first and second minima (e.g. Reiners, Schmitt & Kürster 2001; Reiners & Schmitt 2002).

We computed synthetic spectra (incorporating full Roche-model rotational effects) over a space intended to cover the likely range of parameter values at suitable sampling densities:

v_esin i in the range 270:330 km s⁻¹, at steps of 2 km s⁻¹;

sin i, 0.72:1.00 @ 0.02;

ω_e/ω_c, 0.85:0.99 @ 0.01; and

α, −0.24:+0.24 @ 0.03.

First results for regions around Ca ii K, Mg ii 448.1 nm, and Si ii 634.7 nm showed that significantly subsolar metallicity is required to match observed line strengths, in accord with reports by Gray et al. (2003) and Wu et al. (2011); we obtained reasonable agreement for [M/H] ≃ −0.5, and adopted that value. We then focused on the |$\rm{Si {\small II}}$| λ634.7 line profile for analysis, as it is one of the very few features not to show obvious blending in the spectra employed. (Although least-squares deconvolution is ostensibly capable of addressing the blending issue, and of improving the overall signal:noise, its underpinning principles do not hold when gravity darkening is significant, as is the case here. Moreover, systematics, rather than stochastic noise, prove to dominate uncertainties in the conclusions.)

For each model spectrum, the required values for T_eff and polar radius (a surrogate for polar gravity in these circumstances) were those that reproduce the observed V, UV fluxes, for the matching v_esin i, ω_e/ω_c, and inclination values, as discussed in Section 3.2, but for rigid rotation (regardless of the spectrum-synthesis value of α). This approximation, adopted for computational expedience, is of no consequence for the v_esin i analysis (as is also true for the adopted metallicity).

The comparison between observed and modelled λ634.7 nm profiles was conducted in Fourier space, using the rms differences between normalized transforms⁸ as a test statistic. Precise numerical results depend on the exact frequency (inverse velocity) interval chosen for the comparison, but our general conclusions are insensitive to this, for any reasonable values. We used the range (1.5–4.5) × 10⁻³km⁻¹ s, which encompasses the first two minima in the transform (Fig   6); including the third minimum does not materially change any conclusions, but starts to run into the noise. We found no evidence for any additional broadening processes (‘macroturbulence’) beyond the basic physical mechanisms integral to the modelling.

B3 Ad hoc characterization: an empirical limit on differential rotation

Some results of the initial analysis are summarized in Fig. B1. The basic empirical test for differential rotation is embodied in panel (a), where the minimum rms FT O−C for any (sin i, ω_e/ω_c) combination is shown as a function of v_esin i and α. This figure hints at possibly antisolar differential rotation (i.e. negative α), which would contrast with the handful of positive-α A-star detections reported in the literature (Reiners & Royer 2004; Ammler-von Eiff & Reiners 2012).

However, we find that quite small revisions to the adopted continuum normalization can introduce significant changes to the transform (cf. Dravins, Lindegren & Torkelsson 1990). We label our initial, ‘by eye’, continuum as ‘c0’ in the accompanying figures. Modifying the observed profile by division [resp., multiplication] with a cosine bell of half-width 320 km s⁻¹ and peak amplitude 0.1 per cent of this initial continuum gives a slightly deeper [shallower] line of slightly different shape arising from the slightly higher [lower] continuum, labelled c+ [c−]. The c+ continuum leads to the results shown in Fig. B1(b), which are entirely consistent with solid-body surface rotation.

Continuum uncertainties in the rectified observations are certainly possible at this level (if only because of unrecognized weak line blends). Furthermore, although the comparison model spectra can be rectified simply by division with the corresponding model continuum, in practice this does not lead to a result well suited to comparison to observations. A degree of subjectivity therefore also enters in rectifying the models (even though this was done in an automated procedure), accommodating further potential uncertainty. We conclude that a conservative interpretation of our results is that the initial line-profile analysis alone does not provide any compelling direct evidence for differential rotation in ζ Aql A, and constrains |α| to ≲ 0.05.

B4 Ad hoc characterization: ω_e/ω_c dependence

Fig. B1 also illustrates the sensitivity of v_esin i to other parameters of interest. The line profile offers no useful diagnostic potential for axial inclination, but there is a clear dependence of v_esin i on ω_e/ω_c (at any fixed α; e.g. panel f). This has a straightforward interpretation: as a consequence of gravity darkening, the high-velocity equatorial belt becomes less evident in the spectrum at high ω_e/ω_c, requiring an increase in v_esin i in order to fill in the extreme wings of the line profile (cf. Townsend, Owocki & Howarth 2004).⁹

(in km s⁻¹), where ϖ = ω_e/ω_c − 0.9, and the (a), (b) numerical values are from otherwise identical analyses based on Hipparcos and Gaia parallaxes, respectively (and confirm the expectation of negligible sensitivity of v_esin i to distance). Equation (B2a) is shown as a white line in Fig. B1(f), and represents our adopted characterization of v_esin i as a function of ω_e/ω_c for these α ≡ 0 models (valid over the range 0.85 ≤ ω_e/ω_c ≲ 0.99).

B5 ESTER modelling

As discussed in Section 4.1, an initial analysis based on assumed solid-body surface rotation and the Hipparcos parallax was challenged by disparities between rotation periods implied by the models and the TESS photometric period. We therefore examined the question of differential rotational further, under the constraint of theoretical models of differential surface rotation, rather than an arbitrary ad hoc formulation. To this end, we computed a series of structure models using the ester code¹⁰ (Espinosa Lara & Rieutord 2013; Rieutord, Espinosa Lara & Putigny 2016). Ester computes the stellar structure self-consistently with the radial and latitudinal differential rotation and the meridional circulation resulting from driving by the baroclinic torque (at solar abundance).

Initial parameter modelling indicated M ≃ 2.5 M_⊙, X_c ≃ 0.7, ω_k ≃ 0.8; ester surface-rotation results for this parameter set are shown in Fig. B2 (and were used for most of our subsequent modelling). The sensitivity to parameter variations from this baseline set is also illustrated. Even at rather rapid rotation, only quite modest departures from solid-body surface rotation are predicted, and they are insensitive to precise values of the free parameters (within the range of uncertainty of our results).