# The pulsating magnetosphere of the extremely slowly rotating magnetic β Cep star ξ1 CMa

M. Shultz, G. Wade, Th. Rivinius, C. Neiner, H. Henrichs, W. Marcolino, The Mimes Collaboration

► To cite this version:

M. Shultz, G. Wade, Th. Rivinius, C. Neiner, H. Henrichs, et al.. The pulsating magnetosphere of the extremely slowly rotating magnetic β Cep star ξ1 CMa. Monthly Notices of the Royal Astronomical Society, 2017, 471 (2), pp.2286-2310. 10.1093/mnras/stx1632. obspm-02190692

HAL Id: obspm-02190692

https://hal-obspm.ccsd.cnrs.fr/obspm-02190692v1

Submitted on 13 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# The pulsating magnetosphere of the extremely slowly rotating magnetic β Cep star ξ<fnref n="1" /> CMa

M. Shultz,<fnref n="1" />,<fnref n="2" />,<fnref n="3" />,<fnref n="4" /> G. A. Wade,<fnref n="2" /> Th. Rivinius,<fnref n="3" /> C. Neiner,<fnref n="5" /> H. Henrichs,<fnref n="6" />,<fnref n="7" />
W. Marcolino<fnref n="8" /> and the MiMeS Collaboration

Accepted 2017 June 26. Received 2017 June 21; in original form 2017 March 27

## ABSTRACT

ξ<fnref n="1" /> CMa is a monoperiodically pulsating, magnetic β Cep star with magnetospheric X-ray emission that, uniquely amongst magnetic stars, is clearly modulated with the star's pulsation period. The rotational period *P*<sub>rot</sub> has yet to be identified, with multiple competing claims in the literature. We present an analysis of a large ESPaDOnS data set with a 9 yr baseline. The longitudinal magnetic field ⟨*B*<sub>z</sub>⟩ shows a significant annual variation, suggesting that *P*<sub>rot</sub> is at least of the order of decades. The possibility that the star's H α emission originates around a classical Be companion star is explored and rejected based upon Very Large Telescope Interferometer AMBER and PIONIER interferometry, indicating that the emission must instead originate in the star's magnetosphere and should therefore also be modulated with *P*<sub>rot</sub>. Period analysis of H α equivalent widths measured from ESPaDOnS and CORALIE spectra indicates *P*<sub>rot</sub> > 30 yr. All evidence thus supports that ξ<fnref n="1" /> CMa is a very slowly rotating magnetic star hosting a dynamical magnetosphere. H α also shows evidence for modulation with the pulsation period, a phenomenon that we show cannot be explained by variability of the underlying photospheric line profile, i.e. it may reflect changes in the quantity and distribution of magnetically confined plasma in the circumstellar environment. In comparison to other magnetic stars with similar stellar properties, ξ<fnref n="1" /> CMa is by far the most slowly rotating magnetic B-type star, is the only slowly rotating B-type star with a magnetosphere detectable in H α (and thus, the coolest star with an optically detectable dynamical magnetosphere), and is the only known early-type magnetic star with H α emission modulated by both pulsation and rotation.

**Key words:** stars: early-type – stars: individual: ξ<fnref n="1" /> CMa – stars: magnetic field – stars: mass-loss – stars: oscillations.

## 1 INTRODUCTION

The bright (V = 4.3 mag), sharp-lined β Cep pulsator ξ<fnref n="1" /> CMa (HD 46328, B0.5 IV) was first reported to be magnetic by Hubrig et al. (2006), based upon FORS observations at the Very Large Telescope. Silvester et al. (2009) confirmed the detection using ESPaDOnS at the Canada–France–Hawaii Telescope (CFHT).

Despite intensive observation campaigns with different spectropolarimeters, there are conflicting claims regarding the star's rotational period: Hubrig et al. (2011) derived a rotational period of ~2.18 d based on FORS1/2 and SOFIN data, while Fournune-Ravard et al. (2011) found a longer period of ~4.27 d based on ESPaDOnS measurements. Both Hubrig et al. and Fournune-Ravard et al. agreed that the star was likely being viewed with a rotational pole close to the line of sight. However, Fournune-Ravard et al. noted that given the small variation in ⟨*B*<sub>z</sub>⟩, they could not rule out intrinsically slow rotation.

<footnote n="1">¹Department of Physics, Engineering Physics & Astronomy, Queen's University, Kingston, ON K7L 3N6, Canada</footnote>

<footnote n="2">²Department of Physics, Royal Military College of Canada, Kingston, ON K7K 7B4, Canada</footnote>

<footnote n="3">³ESO – European Organisation for Astronomical Research in the Southern Hemisphere, Casilla 19001, Santiago 19, Chile</footnote>

<footnote n="4">⁴Department of Physics and Astronomy, Uppsala University, Box 516, SE-75120 Uppsala, Sweden</footnote>

<footnote n="5">⁵LESIA, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Universités, UPMC Univ. Paris 06, Univ. Paris Diderot, Sorbonne Paris Cité, 5 place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote n="6">⁶Anton Pannekoek Institute for Astronomy, University of Amsterdam, Science Park 904, NL-1098 XH Amsterdam, the Netherlands</footnote>

<footnote n="7">⁷Department of Astrophysics/IMAPP, Radboud University, Nijmegen, PO Box 9010, NL-6500 GL Nijmegen, the Netherlands</footnote>

<footnote n="8">⁸Universidade Federal do Rio de Janeiro, Observatório do Valongo, Ladeira Pedro AntUnio, 43, CEP 20.080-090 Rio de Janeiro, Brazil</footnote>

<footnote>* E-mail: matthew.shultz@physics.uu.se</footnote>
Frost (1907) was the first to note $\xi^1$ CMa's variable radial velocity (RV). McNamara (1955) found a pulsation period of ~0.2 d, identifying the star as a $\beta$ Cephei variable. While many $\beta$ Cep stars are multi-periodic, $\xi^1$ CMa appears to be a monoperiodic pulsator (Saesen, Briquet & Aerts 2006) with an essentially constant period (Jerzykiewicz 1999). In their analysis of CORALIE high-resolution spectroscopy, Saesen et al. (2006) found a pulsation period of 0.209 5764(4) d (where the number in parentheses refers to the uncertainty in the final digit). Saesen et al. also reported that $\xi^1$ CMa is one of the few $\beta$ Cep stars that achieve supersonic pulsation velocities. Several spectroscopic mode identification methods revealed that the oscillation frequency most likely corresponds to either a radial $(l, m) = (0, 0)$ mode or a dipolar $(l, m) = (1, 0)$ mode, the latter viewed at small inclination. An $l=2$ mode was ruled out via modelling of the velocity moments. Photometric mode identification indicates that only the radial mode agrees with the frequency (Heynderickx, Waelkens & Smeyers 1994).

Like many magnetic early-type stars, $\xi^1$ CMa displays the emission signatures of a magnetized stellar wind in optical, ultraviolet and X-ray spectra. The star's X-ray emission spectrum is very hard, and the brightest amongst all magnetic $\beta$ Cep stars, as demonstrated with Einstein (Cassinelli et al. 1994), ROSAT (Cassinelli et al. 1994) and XMM-Newton (Oskinova et al. 2011). The star's X-ray emission has recently been shown to be modulated with the pulsation period (Oskinova et al. 2014), a so-far unique discovery amongst massive magnetic pulsators (Favata et al. 2009; Oskinova et al. 2015).<fnref n="1" /> Various wind-sensitive ultraviolet lines show strong emission (Schnerr et al. 2008), which is typical for magnetic stars. $\xi^1$ CMa has also been reported to display $H\alpha$ emission (Fournier-Ravard et al. 2011), although the properties of this emission have not yet been investigated in detail.

Our goals in this paper are, first, to determine the rotational period based upon an expanded ESPaDOnS data set, and secondly, to examine the star's $H\alpha$ emission properties.

The observations are presented in Section 2. In Section 3, we refine the pulsation period using RV measurements. In Section 4, we analyse the line broadening, determine the stellar parameters, and investigate the relationship between the star's pulsations and its effective temperature. In Section 5, we examine the possibility that the $H\alpha$ emission may be a consequence of an undetected binary companion, constraining the brightness of such a companion using both interferometry and RV measurements. The magnetic measurements and magnetic period analysis are presented in Section 6. We examine the long- and short-term variability of the $H\alpha$ emission in Section 7, and investigate the influence of pulsation on wind-sensitive UV resonance lines. Magnetic and magnetospheric parameters are determined in Section 8, Section 9 presents a discussion of the paper's results, and the conclusions are summarized in Section 10.

# 2 OBSERVATIONS

## 2.1 ESPaDOnS spectropolarimetry

Under the auspices of the Magnetism in Massive Stars (MiMeS) CFHT Large Program (Wade et al. 2016), 29 ESPaDOnS Stokes V spectra were acquired between 2008/01 and 2013/02. One additional observation was already published by Silvester et al. (2009). A further four ESPaDOnS observations were acquired by a PI programme in 2014, and another 22 by a separate PI programme in 2017.<fnref n="2" /> ESPaDOnS is a fibre-fed echelle spectropolarimeter, with a spectral resolution $\lambda/\Delta\lambda \sim 65\,000$, and a spectral range from 370 to 1050 nm over 40 spectral orders. Each observation consists of four polarimetric sub-exposures, between which the orientations of the instrument's Fresnel rhombs are changed, yielding four intensity (Stokes I) spectra, one circularly polarized (Stokes V) spectrum and two null polarization (N) spectra, the latter obtained in such a way as to cancel out the intrinsic polarization of the source. The majority of the data were acquired using a sub-exposure time of 60 s, with the exception of the first observation (75 s), and the data acquired in 2017, for which a 72 s sub-exposure time was used to compensate for degradation of the coating of CFHT's mirror. The data were reduced using CFHT's Upena reduction pipeline, which incorporates LIBRE-ESPRIT, a descendant of the ESPRIT code described by Donati et al. (1997). Wade et al. (2016) describe the reduction and analysis of MiMeS ESPaDOnS data in detail.

The log of ESPaDOnS observations is provided in an online appendix in Table B1. The quality of the data is excellent, with a median peak signal-to-noise ratio (S/N) per spectral pixel of 828 in the combined Stokes V spectrum. The two observations acquired on 2017/02/11 had a peak S/N below 500, and were discarded from the magnetic analysis. The log of sub-exposures, which were used for spectroscopic analysis, is provided in an online appendix in Table C1. Note that while there are 56 full polarization sequences, due to one incomplete polarization sequence there are 227 rather than 224 sub-exposures listed in online Table C1.

## 2.2 MuSiCoS spectropolarimetry

Three Stokes V spectra were obtained in 2000/02 with the MuSiCoS spectropolarimeter on the Bernard Lyot Telescope (TBL) at the Pic du Midi Observatory. This instrument, one of several similar fibre-fed échelle multi-site continuous spectroscopy (hence MuSiCoS) instruments (Baudrand & Bohm 1992) constructed at various observatories, was uniquely coupled to a polarimeter (Baudrand & Bohm 1992; Donati et al. 1999). It had a spectral resolution of 35 000 and covered the wavelength range 450–660 nm, across 40 spectral orders. As with ESPaDOnS, each polarimetric sequence consisted of four polarized sub-exposures, from which Stokes I and V spectra, as well as a diagnostic null N spectrum, were extracted. The log of MuSiCoS data is provided in an online appendix in Table B1. One of the MuSiCoS observations has a low S/N (below 100), and was discarded from the analysis. The data were reduced using ESPRIT (Donati et al. 1997).

Normal operation of the instrument was verified by observation of magnetic standard stars in the context of other observing programmes (Shorlin et al. 2002; Kochukhov et al. 2004; Wade et al. 2006). In addition to this, we utilize the MuSiCoS spectra of 36 Lyn presented by Wade et al. (2006), one of which was obtained

<footnote n="1">¹ Short-term X-ray variability has also been reported for the magnetic $\beta$ Cep star $\beta$ Cen (Cohen, Cassinelli & MacFarlane 1997; Alecian et al. 2011), and for the non-magnetic (Grunhut, private communication) $\beta$ Cep star $\beta$ Cru (Cohen et al. 2008). However, while subsequent analyses of the X-ray light curves of these stars have confirmed variability, they have not confirmed those variations to be coherent with the primary pulsation frequencies (Raassen et al. 2005; Oskinova et al. 2015); thus, $\xi^1$ CMa is the only $\beta$ Cep star for which X-ray variation is unambiguously coherent with the pulsation period.</footnote>

<footnote n="2">² Programme codes 14AC010 and 17AC16, PI: M. Shultz.</footnote>
during the same observing run as those of $\xi^1$ CMa, in order to compare them to observations of the same star acquired recently with Narval, a clone of ESPaDOnS that replaced MuSiCoS at TBL and that achieves essentially identical results (Wade et al. 2016). This comparison is provided in Appendix A.

## 2.3 CORALIE optical spectroscopy

A large data set (401 spectra) was obtained between 2000/02 and 2004/10 with the CORALIE fibre-fed echelle spectrograph installed at the Nasmyth focus of the Swiss 1.2 m Leonhard Euler Telescope at the European Southern Observatory's (ESO) La Silla facility (Queloz et al. 2000, 2001). The spectrograph has a spectral resolving power of ~100 000, and covers the wavelength range 387–680 nm across 68 spectral orders. The data were reduced with TACOS (Baranne et al. 1996). The first analysis of these data was presented by Saesen et al. (2006).

The log of CORALIE observations is provided in an online appendix in Table C2. The median peak S/N is 205. One observation, acquired on 10/12/2003, was discarded as it had an S/N of 16, too low for useful measurements.

## 2.4 IUE ultraviolet spectroscopy

$\xi^1$ CMa was observed numerous times with the International Ultraviolet Explorer (IUE). The IUE could operate in two modes, high-dispersion ($R \sim 2000$) or low-dispersion ($R \sim 300$), with two cameras, the Short Wavelength from 115 to 200 nm and the Long Wavelength from 185 to 330 nm. We retrieved the data from the MAST archive.<fnref n="3" /> The data were reduced with the New Spectral Image Processing System (NEWSIPS). There are two simultaneous low-resolution spectra obtained with the Short Wavelength Prime (SWP) and Long Wavelength Redundant (LWR) cameras, along with 13 high-resolution spectra obtained with the SWP camera. The high-resolution data are of uniform quality as measured by the S/N, which is approximately 20 in all cases. 12 of the spectra were obtained in close temporal proximity, covering a single pulsation cycle in approximately even phase intervals. The first observation was obtained 154 d previously. We used the absolute calibrated flux, discarded all pixels flagged as anomalous and merged the various spectral orders.

## 2.5 VLTI near-infrared interferometry

While the angular radii of $\xi^1$ CMa and its circumstellar material are certainly too small to be resolved interferometrically, the star is bright enough that the data can be used to search for a high-contrast binary companion.

We have acquired four Very Large Telescope Interferometer (VLTI) observations: low-resolution AMBER *H* and *K* photometry, a high-resolution AMBER NIR spectro-interferogram and one low-resolution *H*-band PIONIER observation. AMBER (Astronomical Multi-Beam Recombiner) offers three baselines and can operate in either low-resolution photometric mode or high-resolution (*R* ~ 12 000) spectro-interferometric mode (Petrov et al. 2007). PIONIER (Precision Integrated-Optics Near-infrared Imaging ExperIMENT) combines light from four telescopes, offering visibilities across six baselines together with four closure phase measurements (Le Bouquin et al. 2011). All data were obtained using the 1.8 m

Auxiliary Telescopes, which have longer available baselines, and hence better angular resolution than the 8 m Unit Telescopes. AMBER observations were obtained with baselines ranging from 80 to 129 m. PIONIER was configured in the large quadruplet, with a longest baseline of 140 m. The nearby star $\xi^2$ CMa (A0 III, *H* = 4.63) was used as a standard star for calibration.

The observing log is given in Table D1. As the execution time for a full measurement (~1 h) is a significant fraction of the pulsation period, pulsation phases are not given.

Owing to the relatively faint magnitude, the wavelength edges of the AMBER *H* and *K* profiles were particularly noisy. These data points were edited out by hand before commencing the analysis. No such procedure was necessary for the PIONIER observation.

# 3 PULSATION PERIOD

## 3.1 Radial velocities

We measured RVs from the ESPaDOnS and CORALIE spectra using the centre-of-gravity method (Aerts, Christensen-Dalsgaard & Kurtz 2010). In addition to using the same line (Si II 455.3 nm) used by Saesen et al. (2006), we have used N II 404.4 and 422.8 nm; N III 463.4 nm; O II 407.2, 407.9, 418.5 and 445.2 nm; Ne II 439.2 nm; Al III 451.3 nm; and Si IV 411.6 nm. For ESPaDOnS data, we used individual sub-exposures rather than the Stokes *I* profiles corresponding to the full polarization sequences: each polarization sequence encompasses about 2.3 per cent of a pulsation cycle (4 × 60 s sub-exposures + 3 × 60 s chip readouts), as compared to about 0.33 per cent for the sub-exposures (although the 2017 data had slightly longer sub-exposure times, this was compensated for by shorter readout times). The CORALIE data are not as uniform as the ESPaDOnS data in this regard: the median exposure time corresponds to about 2.8 per cent of a pulsation cycle, and some are up to about 10 per cent (all CORALIE measurements were retained; however, HJDs were calculated at the middle rather than the beginning of each exposure).

Since centre-of-gravity measurements can be biased by inclusion of too much continuum, it is important to choose integration limits with care. Thus RVs were measured iteratively. A first set of measurements was conducted using a wide integration range (±50 km s⁻¹ about the systemic velocity), chosen to encompass the full range of variation. These initial RVs were then used as the central velocities for a second set of measurements, with an integration range of ±30 km s⁻¹. These were in turn used a final time as the central velocities to refine the RVs with a third iteration using the same integration range. The error bar weighted mean RV across all lines was then taken as the final RV measurement for each observation, with the standard deviation across all lines as the uncertainty in the final RV. ESPaDOnS and CORALIE RV measurements are tabulated in an online appendix in Tables C1 and C2, respectively. We determined the systemic velocity from the mean RV across all observations to be 22.5 ± 1 km s⁻¹.

## 3.2 Frequency analysis

We analysed the RVs using PERIOD04 (Lenz & Breger 2005). Frequency spectra are shown in Fig. 1, where the top panels show the frequency spectra for the original data set, and the panels below show frequency spectra after pre-whitening with the most significant frequencies from previous frequency spectra. An S/N threshold of 4 was adopted as the minimum S/N for significance (Breger

<footnote n="3">³ Available at https://archive.stsci.edu/iue/</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Frequency spectra for RVs measured from CORALIE data (left), ESPaDOnS data (middle) and the combined (right) data set, for the RVs (top row), and after pre-whitening with significant frequencies (subsequent rows). Dotted light blue lines indicate $f_{\text{pul}} \sim 4.771\,53\,\text{d}^{-1}$ and its first two harmonics. In the bottom panels, the window function is overplotted in light blue. In the bottom-right panel, the white frequency spectrum shows the results of pre-whitening the combined data with a pulsation period increasing at a linear rate of 0.0096(5) s yr$^{-1}$. By contrast, using a constant period to pre-whiten the combined data set (black frequency spectrum) does not fully remove the peak at $f_{\text{pul}}$.</figcaption></figure>

et al. 1993; Kuschnig et al. 1997). The uncertainty in each frequency was determined using the formula from Bloomfield (1976), $\sigma_F = \sqrt{6\sigma_{\text{obs}}/(\pi\sqrt{N_{\text{obs}}A}\Delta T)}$, where $\sigma_{\text{obs}}$ is the mean uncertainty in the RV measurements, $N_{\text{obs}}$ is the number of measurements, A is the amplitude of the RV curve and $\Delta T$ is the time-span of observations.

Period analysis of the CORALIE data set (left-hand panels in Fig. 1) yielded the same results as those reported by Saesen et al. (2006), with significant frequencies at $f_{\text{pul,COR}} = 4.771\,5297(5)\,\text{d}^{-1}$ with an amplitude of $16.6\,\text{km s}^{-1}$, and at the harmonics $2f_{\text{pul,COR}}$ and $3f_{\text{pul,COR}}$ with amplitudes of 0.6 and $0.3\,\text{km s}^{-1}$, respectively. After pre-whitening with these frequencies, no significant frequencies remain (bottom-left panel of Fig. 1), and all peaks are at the $1\,\text{d}^{-1}$ aliasing of the spectral window. The same analysis of the ESPaDOnS data (middle panels in Fig. 1) finds maximum power at $f_{\text{pul,ESP}} = 4.771\,5007(4)\,\text{d}^{-1}$, and at $2f_{\text{pul,ESP}}$.

The combined data set (right-hand panels in Fig. 1) yields the strongest signal at $f_{\text{pul,Comb}} = 4.771\,5121(3)\,\text{d}^{-1}$, along with significant peaks at $2f_{\text{pul,Comb}}$ and $3f_{\text{pul,Comb}}$. In this case, pre-whitening does not completely remove the peak corresponding to $f_{\text{pul,Comb}}$. This could indicate that the pulsation frequency is not the same between the data sets. The difference in frequencies between the CORALIE and ESPaDOnS data sets, $2.9 \times 10^{-5}\,\text{d}^{-1}$, is about

<figure><img src="image_4.png" /><figcaption>**Figure 2.** Top panel: composite RVs (using $N_{\text{II}}$ 404.4 and 422.8 nm; $N_{\text{III}}$ 463.4 nm; $O_{\text{II}}$ 407.2, 407.9, 418.5 and 445.2 nm; $Ne_{\text{II}}$ 439.2 nm; $Al_{\text{III}}$ 451.3 nm; $Si_{\text{III}}$ 455.3 nm; and $Si_{\text{IV}}$ 411.6 nm) as a function of pulsation phase, using the non-linear ephemeris defined in equation (1). Open black circles are CORALIE measurements, and filled blue circles ESPaDOnS measurements. The solid red line indicates a third-order sinusoidal fit. Bottom panel: residual RVs. Dashed red lines indicate $\pm 1$ standard deviation, about 0.48 km s$^{-1}$, and the dotted red lines show the median RV uncertainty.</figcaption></figure>

30 times larger than the formal uncertainties in either frequency. This difference corresponds to an increase in the pulsation period of $0.1\,\text{s}$, approximately compatible with the increase of $0.063 \pm 0.009\,\text{s}$ expected from the constant period change of $+0.0037(5)\,\text{s yr}^{-1}$ reported by Jerzykiewicz (1999).

To explore this hypothesis, phases were calculated assuming a constant rate of period change $-0.02 < \dot{P} < +0.02\,\text{s yr}^{-1}$, with $f_0 = 4.771\,529(7)\,\text{d}^{-1}$ taken as the frequency inferred from the CORALIE data acquired in 2000, and allowed to vary within the uncertainty in this frequency. The phases $\phi$ were calculated as

$$ \phi = \frac{\text{HJD} - (T_0 + N P_0 + 0.5 \dot{P} N^2)}{P_0 + N \dot{P}}, \qquad (1) $$

where $P_0$ is the initial period, and $N$ is the number of pulsation cycles elapsed between the observation and the reference epoch $T_0 = 2451591.42576$, defined as the time of the first RV maximum one cycle before the first observation in the data set. The goodness-of-fit (GOF) $\chi^2$ statistic was calculated for each combination of $P_0$ and $\dot{P}$ using a third-order least-squares sinusoidal fit in order to account for the pulsation frequency and its first two harmonics. The minimum $\chi^2$ solution was found for $P_0 = 0.209\,5763(1)\,\text{d}$ and $\dot{P} = +0.0096(5)\,\text{s yr}^{-1}$, where the uncertainties were determined from the range over which $\chi^2$ does not change appreciably. The RVs are shown phased with equation (1) using these parameters in Fig. 2. The bottom panel of Fig. 2 shows the residual RVs after subtraction of the third-order sinusoidal fit. The standard deviation
of the residuals is 0.48 km s⁻¹, an improvement over the standard deviation of 0.94 km s⁻¹ obtained using the constant period determined from the full data set. The bottom-right panel of Fig. 1 shows the frequency spectrum obtained for the residual RVs in Fig. 2: in contrast to the results obtained via pre-whitening using a constant pulsation frequency, all power at *f*<sub>pul</sub> and its harmonics is removed, with no significant frequencies remaining.

### 3.3 Comparison to previous results

Our analysis of the CORALIE data set recovers the same pulsation frequency as that found by the analysis performed by Saesen et al. (2006) of the same data. However, we find *Ṗ* to be almost three times larger than the rate found by Jerzykiewicz (1999). Whether this reflects an acceleration in the rate of period change will need to be explored in the future when larger data sets with a longer temporal baseline are available. Neilson & Ignace (2015) noted that ξ¹ CMa is one of the few β Cep stars for which *Ṗ* is low enough to be consistent with stellar evolutionary models. The higher *Ṗ* suggested by our results brings ξ¹ CMa closer to the range observed for other β Cep stars, which is to say, slightly above the *Ṗ* predicted by evolutionary models for a star of ξ¹ CMa's luminosity.

## 4 STELLAR PARAMETERS

This section presents an analysis of the fundamental stellar parameters of ξ¹ CMa. The results are summarized in Table 1.

### 4.1 Line broadening

The high S/N, high spectral resolution and wide wavelength coverage of the ESPaDOnS data, combined with the sharp spectral lines of this star, present numerous opportunities for constraining the line broadening mechanisms. As with measurement of RVs, in order to minimize the impact of RV variation, Stokes *I* spectra from individual sub-exposures were used rather than the spectra computed from the combined exposures. We selected two lines for analysis: the Si III 455.3 nm line, which is sensitive to the pulsational properties

<table><caption>**Table 1.** Observational and physical properties. The third column gives the reference, where numerical references indicate the section within this paper where the parameter is determined.</caption><tr><td>v<sub>sys</sub> (km s<sup>-1</sup>)</td><td>22.5 ± 1</td><td>3</td></tr><tr><td>v<sub>sin i</sub> (km s<sup>-1</sup>)</td><td>≤8</td><td>4.1</td></tr><tr><td>v<sub>mac</sub> (km s<sup>-1</sup>)</td><td>8 ± 2</td><td>4.1</td></tr><tr><td>Sp. Type</td><td>B0.5 IV</td><td>Hubrig et al. (2006)</td></tr><tr><td>V (mag)</td><td>4.3</td><td>Ducati (2002)</td></tr><tr><td>π (mas)</td><td>2.36 ± 0.2</td><td>van Leeuwen (2007)</td></tr><tr><td>d (pc)</td><td>424 ± 35</td><td>van Leeuwen (2007)</td></tr><tr><td>BC (mag)</td><td>-2.61 ± 0.09</td><td>4.5</td></tr><tr><td>A<sub>V</sub> (mag)</td><td>&lt;0.06</td><td>4.5</td></tr><tr><td>m<sub>V</sub> (mag)</td><td>-3.86 ± 0.18</td><td>4.5</td></tr><tr><td>M<sub>bol</sub> (mag)</td><td>-6.47 ± 0.27</td><td>4.5</td></tr><tr><td>log (L/L<sub>☉</sub>)</td><td>4.49 ± 0.11</td><td>4.5</td></tr><tr><td>T<sub>eff</sub> (kK)</td><td>27 ± 1</td><td>4.3</td></tr><tr><td>log g</td><td>3.78 ± 0.07</td><td>4.2</td></tr><tr><td>R<sub>*</sub> (R<sub>☉</sub>)</td><td>7.9 ± 0.6</td><td>4.3</td></tr><tr><td>M<sub>*</sub> (M<sub>☉</sub>)</td><td>14.2 ± 0.4</td><td>4.6</td></tr><tr><td>Age (Myr)</td><td>11.1 ± 0.7</td><td>4.6</td></tr><tr><td>τ<sub>MS</sub></td><td>0.77 ± 0.04</td><td>4.6</td></tr><tr><td>R<sub>*</sub>/R<sub>ZAMS</sub></td><td>1.8 ± 0.1</td><td>4.6</td></tr></table>

of β Cep stars (Aerts & De Cat 2003), and for comparison N II 404.4 nm, a weaker line with lower sensitivity to pulsation.

We applied a model incorporating the projected rotation velocity vsin *i*, radial/tangential macroturbulence v<sub>mac</sub> and assumed radial pulsations, performing a GOF test on a grid spanning 0–20 km s⁻¹ in vsin *i* and 0–30 km s⁻¹ in v<sub>mac</sub>, in 1 km s⁻¹ increments. The disc integration model is essentially as described by Petit & Wade (2012), with the exception of radial pulsations. Pulsations were modelled as a uniform velocity component normal to the photosphere, with the local line profiles in each surface area element shifted by the projected line-of-sight component of the pulsation velocity. A limb darkening coefficient of ϵ = 0.36 was used, obtained from the tables calculated by Díaz-Cordovés, Claret & Giménez (1995) for a star with T<sub>eff</sub> =27 kK and log g = 3.75. Local profiles were broadened with a thermal velocity component of 4 km s⁻¹ for Si III, and 5.6 km s⁻¹ for N II, and the disc-integrated profile was convolved according to the resolving power of ESPaDOnS. Line strength was set by normalizing the equivalent widths (EWs) of the synthetic line profiles to match the EWs of the observed line profiles. In the course of this analysis, we found a Baade–Wesselink projection factor of measured-to-disc-centre RV of *p* = 1.45 ± 0.02, which is consistent with determinations for other β Cep variables (e.g. Nardetto et al. 2013).

Fig. 3 shows the resulting best-fitting models. Taking the mean and standard deviation across all observations, both Si III 455.3 nm and N II 404.4 nm yield vsin *i* = 5 ± 3 km s⁻¹. The two lines give different values for v<sub>mac</sub>, however: 19 ± 1 and 8 ± 2 km s⁻¹, respectively. The much lower value of v<sub>mac</sub> determined from N II 404.4 nm may indicate that the extended wings of Si III 455.3 nm, which require a higher value of v<sub>mac</sub> to fit, are at least partly a consequence of pulsations.

As demonstrated in Fig. 3, there is essentially no difference in quality of fit between models with non-zero vsin *i* and v<sub>mac</sub>, and zero vsin *i*; conversely, a much worse fit is obtained by setting v<sub>mac</sub> = 0 km s⁻¹. For Si III 455.3 nm, the χ² of the model without v<sub>mac</sub> is 13 times higher than the χ² of the model without vsin *i*, while the χ² of the model without vsin *i* is only 1.4 times higher than the model with both turbulent and rotational broadening. The difference is not as great for N II 404.4 nm, although the χ² of the model with rotational broadening only is still higher than either the model with turbulent broadening only or both rotational and turbulent broadening. Sundqvist et al. (2013) used magnetic O stars with rotational periods known to be extremely long (i.e. with vsin *i* ~0 km s⁻¹) and found that both the GOF test and the Fourier transform methods severely overestimate vsin *i* when v<sub>mac</sub> > vsin *i*, as is the case for both of the lines we have examined. We conclude that the star’s line profiles are consistent with vsin *i* = 0 km s⁻¹.

### 4.2 Surface gravity

The surface gravity log g was determined by fitting TLUSTY BSTAR2006 synthetic spectra (Lanz & Hubeny 2007) to the H β and H γ lines of the ESPaDOnS spectrum acquired at pulsation phases at which the RV was closest to the systemic velocity of 22.5 km s⁻¹ (i.e. at pulsation phases 0.25 or 0.75). As both of these lines are close to the edges of their respective orders, in order to avoid warping the line, the orders were merged from unnormalized spectra, then normalized using a linear fit to nearby continuum regions. Five surface gravities were tested between log g = 3.25 and 4.25, and a low-order polynomial was fitted to the resulting χ² in order to identify the lowest χ² solution. Uncertainties were determined by fitting models with T<sub>eff</sub> of 26, 27 and 28 kK, spanning the
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Line profile fitting to determine $v\sin i$ and $v_{\text{mac}}$, using Si II 455.3 nm (top) and N II 404.4 nm (bottom). Left: $\chi^2$ landscapes (high values are dark, low values are light, where the reduced $\chi^2$ goes from approximately 1 to $2 \times 10^5$), corresponding to the weighted mean $\chi^2$ across all observations. The best-fitting solutions for individual observations are indicated by open red circles. Right: line profile fits at pulsation phases close to the RV minima and maxima ($\phi = 0.498, 0.982$) and RV = $0 \text{ km s}^{-1}$ ($\phi = 0.279, 0.750$). Observations are shown by black circles. The solid red line shows the best overall fit for each observation. The dashed blue line shows the best-fitting model with $v\sin i = 0 \text{ km s}^{-1}$. The dash-dotted green line indicates the best-fitting model with $v_{\text{mac}} = 0 \text{ km s}^{-1}$. While the overall best fit has a non-zero $v\sin i$ of a few km s$^{-1}$, excluding $v_{\text{mac}}$ results in a noticeably worse fit, but excluding $v\sin i$ does not substantially worsen the fit. The weaker N II line yields a lower $v_{\text{mac}}$, but does not change the result that $v_{\text{mac}} > v\sin i$, or that equivalently good fits can be achieved with or without $v\sin i$.</figcaption></figure>

approximate range in the uncertainty in $T_{\text{eff}}$ (see below). Fig. 4 shows the best-fitting models, $\log g = 3.78 \pm 0.07$, compared to the observed H $\beta$ line. Fitting to H $\gamma$ yielded identical results.

## 4.3 Effective temperature

The effective temperature $T_{\text{eff}}$ was determined using photometry, spectrophotometry and spectroscopy. Using the Strömgren *uvbyβ* photometric indices obtained by Hauck & Mermilliod (1998), and the IDL program *UVBYBETA.PRO*<fnref n="4" /> that implements the calibrations determined by Moon (1985), yields $T_{\text{eff}}=26.2$ kK. Using the Johnson *UBVRIJHK* photometry collected by Ducati (2002) and the colour-temperature calibration presented by Worthey & Lee (2011) yields $T_{\text{eff}}=25 \pm 3$ kK, where the uncertainty was determined from the range of $T_{\text{eff}}$ across the (U–B), (V–K), (B–V), (V–R) and (J–K) colours.

The $T_{\text{eff}}$ was measured using spectrophotometric data by fitting synthetic BSTAR2006 spectra to the IUE low-dispersion spectra, Johnson *UBVRI* photometry and Strömgren *uvbyβ* photometry. The Johnson photometry was converted into absolute flux units using the calibration provided by Bessell, Castelli & Plez (1998), and the Strömgren photometry converted using the calibration determined by Gray (1998). The synthetic spectra were scaled by $1/d^2$ and $4\pi R^2$, where the distance $d$ was fixed by the Hipparcos parallax and the stellar radius $R_*$ was left as a free parameter. The resulting fit is shown in Fig. 5. The best-fitting model yields $T_{\text{eff}}=26.8 \pm 0.4$ kK and $R_* = 7.8$ R⊙. The uncertainty comes from the spread in

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Surface gravity determination using H $\beta$. Top: model fits for $T_{\text{eff}} = 26, 27$ and 28 kK (cyan, blue and pink, respectively) to the observed profile (black). The observation nearest $\phi = 0.75$ (where the RV=0 km s$^{-1}$) was used. Bottom: residuals. The horizontal lines indicate the mean flux uncertainty. As explored in Section 7.1, there is emission near the line core; this is more prominent in H $\alpha$. This region was not included in the fit.</figcaption></figure>

<footnote n="4">⁴ Available at http://idlastro.gsfc.nasa.gov/ftp/pro/astro/uvbybeta.pro</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 5.** Low-resolution *IUE* spectroscopy (solid lines) with the SWP (black) and LWR (blue) cameras, and photometry (filled symbols) compared to a synthetic photospheric spectrum (red dashed line). No dereddening correction has been applied to the data. Inset: comparison between synthetic spectra (dashed lines) and observed spectra (solid lines) dereddened using three values of $E(B-V)$: 0 (green), 0.02 (cyan) and 0.04 (pink). Both the observed spectrum and the best-fitting synthetic spectra have been smoothed with a boxcar algorithm. The dot-dashed (purple) line is the reddening curve in arbitrary units. Even mild dereddening produces a bulge in the vicinity of the 220 nm silicate absorption feature.</figcaption></figure>

$T_{\text{eff}}$ over the range of $\log g = 3.78 \pm 0.07$, and for two values of $E(B-V)$, 0 and 0.02. Fluxes were dereddened using the *IDL* routine *FM_UNRED*,<fnref n="5" /> which utilizes the Fitzpatrick & Massa (1990) extinction curve parametrization. As demonstrated by the inset in Fig. 5, there is a clear absence of interstellar silicate absorption near 220 nm. The dereddened spectra in the inset in Fig. 5 assume that the strength of the absorption bump is $c_3 = 2.1$, the minimum Galactic value (Clayton et al. 2003). Even with this low value of $c_3$, using $E(B-V) = 0.04$ results in an increase of the flux near the silicate absorption bump that is not observed. Using the average Galactic value of $c_3 = 3.23$ instead would require that $E(B-V)$ be even lower.

## 4.4 Comparison to previous results

Previous measurements of $v\sin i$ have in general yielded somewhat higher, non-zero results than those obtained here ($15 \pm 1.5$ km s$^{-1}$, Saesen et al. 2006; $9 \pm 2$ km s$^{-1}$, Lefever et al. 2010; $14 \pm 5$ km s$^{-1}$, Aerts et al. 2014). Saesen et al. did not include $v_{\text{mac}}$ in their profile fit, but did include pulsation and thermal broadening; Lefever et al. did not consider pulsations; and Aerts et al. included $v_{\text{mac}}$, but neither pulsation velocity nor thermal broadening. Our inclusion of thermal, macroturbulent, pulsational and rotational velocity fields likely explains why our best-fitting value of $v\sin i$, $5 \pm 3$ km s$^{-1}$, is somewhat lower than the values found previously, since all four sources of line broadening are of a similar magnitude in this star.

Stellar parameters for $\xi^1$ CMa have been determined numerous times using a variety of different spectral modelling methodologies, e.g. DETAIL/SURFACE (Morel, Hubrig & Briquet 2008), FASTWIND (Lefever et al. 2010) and Potsdam Wolf-Rayet (PoWR; Oskinova et al. 2011). Niemczura & Daszyńska-Daszkiewicz (2005) employed an algorithmic method to simultaneously determine reddening, $T_{\text{eff}}$, metallicities and surface gravities for β Cep stars observed with *IUE*, and for $\xi^1$ CMa found $T_{\text{eff}} = 24 \pm 1$ kK, $\log g = 3.89$ and $E(B-V) < 0.01$, i.e. their analysis favoured a cooler, slightly less evolved star, but with the same very

low level of reddening. Niemczura & Daszyńska-Daszkiewicz derived $\log g$ from photometric calibrations, rather than line profile fitting, which likely explains the different values of $\log g$. Morel et al. (2008) found $T_{\text{eff}} = 27.5$ kK and $\log g = 3.75$, Lefever et al. (2010) derived $T_{\text{eff}} = 27 \pm 1$ kK and $\log g = 3.80 \pm 0.15$, and Oskinova et al. (2011) found $T_{\text{eff}} = 27.5 \pm 0.5$ kK; all are in agreement with the results found here. The slightly higher $T_{\text{eff}}$ found in these works is likely a consequence of the higher value of reddening that, as is shown in Section 4.3 and Fig. 5, is not consistent with the absence of silicate absorption. While reddening has an effect on the spectrophotometric method applied in Section 4.3, we obtained essentially identical results using line strength ratios and ionization balances below in Section 4.5, which should not be affected by reddening.

## 4.5 Effective temperature variation

The photospheric temperature of a β Cephei star is variable, a property that can be explored both photometrically and spectroscopically. The only public two-colour photometry for $\xi^1$ CMa is the Tycho BV$_T$ observations. While pulsational modulation is visible in both passbands, the precision of the data is not sufficient to detect a coherent signal in the colours. This is not surprising, as B-V colours are poorly suited to determining $T_{\text{eff}}$ for such hot stars. We therefore attempted to constrain the $T_{\text{eff}}$ change via the ionization balances of various spectral lines, under the assumption that all EW changes are due to temperature variation.

We first measured EWs from the ESPaDOnS observations, and from appropriate synthetic spectra obtained from the TLUSTY BSTAR2006 library (Lanz & Hubeny 2007). We used the ESPaDOnS Stokes *I* spectra obtained from the combined sub-exposures, in order to maximize the S/N. The following spectral lines were used: He I 587.6 nm, He II 468.6 nm, C II 426.7 nm, C III 406.8 nm, N II 404.4 nm, N III 463.4 nm, O I 777.4 nm, O II 407.9 nm, Ne I 640.2 nm, Ne II 439.2 nm, Si II 413.1 nm, Si III 455.3 nm and Si IV 411.6 nm. Since the EW can depend on $\log g$ as well as $T_{\text{eff}}$, the model EW ratios were linearly interpolated to $\log g = 3.78$ between the results for $\log g = 3.75$ and 4.0. We then determined $T_{\text{eff}}$ by interpolating through the resulting grid to the

<footnote n="5">⁵ Available at http://idlastro.gsfc.nasa.gov/ftp/pro/astro/fm_unred.pro</footnote>
<figure><img src="image_8.png" /><figcaption>**Figure 6.** $T_{\text{eff}}$ determined from EW ratios of various ion pairs as a function of pulsation phase. While the mean $T_{\text{eff}}$ differs by ~2 kK between different ion pairs, a coherent and identical modulation with the pulsation period is seen in all cases. The semi-amplitude of the $T_{\text{eff}}$ variation, as determined from the least-squares sinusoidal fit to the data, is given at the bottom of each panel. The bottom-right panel shows the weighted mean $T_{\text{eff}}$ across all ion pairs, with the dashed blue line indicating the $T_{\text{eff}}$ determined from SED fitting (Fig. 5).</figcaption></figure>

EW ratios measured from the observed spectra. The results are illustrated in Fig. 6. While there is a considerable spread in values, from approximately 25 to 28 kK, the mean value of 26.8 kK agrees well with the spectrophotometric determination. A coherent variation with the pulsation period is seen for all ion strength ratios. Individual atomic species yield variations with semi-amplitudes up to $\Delta T_{\text{eff}} = \pm 500$ K. Taking the weighted mean across the effective temperatures from all ion strength ratios at each phase yields a variation of $\Delta T_{\text{eff}} = \pm 310 \pm 30$ K. If the zero-point of each $T_{\text{eff}}$ variation is first subtracted, the resulting mean $T_{\text{eff}}$ semi-amplitude is 330 $\pm$ 30 K, slightly larger but overlapping within the uncertainty.

To check the validity of the temperature variation, we used it to model the *Hipparcos* light curve, as demonstrated in Fig. 7. With the uncertainty in the pulsation period and rate of period change determined in Section 3, the accumulated uncertainty in pulsation

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Predicted variation in visual magnitude (red line) as determined from the variation in $R_*$ from RVs (Fig. 2) and in $T_{\text{eff}}$ from EW ratios (Fig. 6) as compared to the *Hipparcos* photometry (black circles). Pulsation phases were calculated using equation (1).</figcaption></figure>

phase for the *Hipparcos* photometry is between 0.009 and 0.012 cycles. We began by determining the mean radius $R_*$ from the mean $T_{\text{eff}}$ and the luminosity $\log L$. $\log L$ was obtained from the *Hipparcos* parallax distance $d = 424 \pm 35$ pc, the apparent $V$ magnitude 4.06 and the bolometric correction BC $= -2.61 \pm 0.09$ mag obtained by linear interpolation through the theoretical BSTAR2006 grid according to $T_{\text{eff}}$ and $\log g$ (Lanz & Hubeny 2007). The absolute $V$ magnitude $m_V = V - A_V - \mu = -3.86 \pm 0.18$ mag, where $\mu = 5\log d - 5 = 8.13 \pm 0.18$ mag is the distance modulus and $A_V = E(B-V)/3.1 < 0.06$ mag is the extinction. The bolometric magnitude is then $M_{\text{bol}} = m_V + \text{BC} = -6.47 \pm 0.27$ mag, yielding $\log(L/L_{\odot}) = (M_{\text{bol,}\odot} - M_{\text{bol}})/2.5 = 4.49 \pm 0.11$ where $M_{\text{bol,}\odot} = 4.74$ mag. This yields $R_* = \sqrt{(L/L_{\odot})/(T_{\text{eff}}/T_{\text{eff,}\odot})}^4 = 7.9 \pm 0.6$ R$_{\odot}$, where $T_{\text{eff,}\odot} = 5.78$ KK. This radius agrees well with the value found via spectrophotometric modelling, $R_* = 7.8$ R$_{\odot}$, in which the radius was left as a free parameter (Fig. 5).

Integrating the RV curve (Fig. 2) in order to obtain the absolute change in radius, and assuming radial pulsation, yields a relatively small change in radius of $\pm 6.58 \times 10^4$ km, corresponding to between 1.0 and 1.5 per cent of the stellar radius. The radius variation was then combined with the $T_{\text{eff}}$ variation from Fig. 6 to obtain $\log L$ and the BC at each phase. The apparent $V$ magnitude was then obtained by reversing the calculations in the previous paragraph. The solid red line in Fig. 7 shows the resulting model light curve compared to the *Hipparcos* photometry, where we made the assumption that $V$ and $H_p$ are approximately equivalent. The larger 500 K $T_{\text{eff}}$ semi-amplitude measured from individual pairs of ions is not consistent with the light curve, yielding a photometric variation larger than observed, with a semi-amplitude of 0.037 mag as compared to the observed semi-amplitude of ~0.021 mag. The apparent phase offset of ~0.05 cycles between the predicted and observed photometric extrema may suggest that $\dot{P}$ may not in fact be constant; alternatively, there may simply be too few measurements to constrain the $T_{\text{eff}}$ variation well enough to obtain a close fit to the photometric data, as the latter is clearly not a perfect sinusoid.

## 4.6 Mass and age

Fig. 8 shows $\xi^1$ CMa’s position on the Hertzsprung–Russell diagram (HRD, top) and the $T_{\text{eff}}-\log g$ diagram (bottom), where the mean $T_{\text{eff}}$ was used. Comparison to the evolutionary models calculated by Ekström et al. (2012, which assume an initial rotational velocity of 40 per cent of the critical velocity) indicates that the stellar mass $M_* = 14.2 \pm 0.4$ M$_{\odot}$, and that the absolute stellar age is $t = 11 \pm 0.7$ Myr. The Ekström et al. (2012) models do not include the effects of magnetic fields; however, grids of self-consistent
<figure><img src="image_10.png" /><figcaption>**Figure 8.** ξ¹ CMa's position on the HRD (top) and on the T<sub>eff</sub>-log g diagram (bottom), relative to other magnetic stars with similar stellar parameters. The solid line indicates the zero-age main sequence (ZAMS), and the dot-dashed line the terminal age main sequence. Dashed lines show evolutionary tracks from the (rotating, non-magnetic) models by Ekström et al. (2012), and dotted lines isochrones from the same models for log (t/yr) = 7.0 (left) and 7.1 (right). The star's fundamental parameters indicate M* = 14.2 ± 0.4 M<sub>☉</sub> and t = 11 ± 0.7 Myr. ξ¹ CMa is the most evolved magnetic B-type star above about 12 M<sub>☉</sub>.</figcaption></figure>

evolutionary models including these effects in a realistic fashion are
not yet available.

The positions of the star on the two diagrams are mutually con-
sistent. The other known magnetic B-type stars with similar stellar
parameters are also shown in Fig. 8. The stellar parameters of the
majority of the other stars were obtained from the catalogue of
magnetic hot stars published by Petit et al. (2013) and references
therein; those of β CMa and ε CMa were obtained from Fossati et al.
(2015a). ξ¹ CMa is one of the most evolved stars in the ensemble,
and is the most evolved star with a mass above about 12 M⊙. Its
position on the T<sub>eff</sub>-log g diagram indicates that it may be a more
evolved analogue of NU Ori, HD 63425 and HD 66665.

As an additional check on the stellar mass and radius, we utilized
the pulsation period and the relation P<sub>fun</sub> = Q/√ρ, where P<sub>fun</sub> is the
fundamental pulsation period, ρ is the mean density and Q is the
pulsation constant. With M* = 14.2 ± 0.4 M<sub>☉</sub>, R* = 7.9 ± 0.6 R<sub>☉</sub>

<figure><img src="image_11.png" /><figcaption>**Figure 9.** Low-resolution VLTI-AMBER *H* and *K* data and VLTI-PIONIER data, showing squared visibilities (top) and closure phases (bottom).</figcaption></figure>

and the theoretical value $Q = 0.035$ (Davey 1973; Stothers 1976),
this yields $P_{\text{fun}} = 0.18 \pm 0.02$ d, very close to the true pulsation
period $P \sim 0.21$ d. If the actual pulsation period is used to calculate
$Q$, we obtain $Q = P_{\text{fun}}\sqrt{\rho} = 0.041 \pm 0.003$.

5 BINARITY

5.1 Interferometry

Fourtune-Ravard et al. (2011) reported that ξ¹ CMa hosts weak H α emission, which is atypical for a star of ξ¹ CMa’s stellar properties. There are two possibilities to explain this: first, that it originates in a stellar magnetosphere, and secondly, that it originates in the decretion disc of a heretofore undetected Be companion star. It seems reasonable to expect that ξ¹ CMa has a binary companion, as the binary fraction is ~65 per cent for B0 stars (Chini et al. 2012). Furthermore, early claims that the H α emission of β Cep, a similar magnetic early B-type pulsator, originated in its magnetosphere (Donati et al. 2001) proved unfounded following the detection of a Be companion star by Schnerr et al. (2006).

The peak H α emission has a maximum strength of 28 per cent of the continuum (see Section 7.1). Be star emission lines can range up to ~10× of the continuum, although the emission is typically much less than this. We therefore expect that the putative companion must have a luminosity of ≥2.8 per cent of that of ξ¹ CMa, yielding log (L/L<sub>☉</sub>) ≥ 2.9 or V < 7.6. Assuming the pair to be coeval, and therefore locating the star on the same isochrone as ξ¹ CMa, the companion’s mass should be near 5 M<sub>☉</sub>, compatible with a classical Be star.

ξ¹ CMa is listed in the Washington Double Star Catalog as possessing a binary companion (Mason et al. 2001). However, the reported companion is both too dim to account for the H α emission (V = 14) and too far away, located approximately 28 arcsec from the primary, i.e. well outside the 1.6 arcsec ESPaDOnS aperture.

We acquired *H* and *K* low-resolution VLTI-AMBER data, together with VLTI-PIONIER data, in order to search for the presence of a binary star. The squared visibilities *V*² and closure phases *φ*<sub>C</sub> of these data are shown as functions of spatial frequency in Fig. 9. Both the AMBER *K*-band and PIONIER *V*² measurements are entirely flat, compatible with a single unresolved source. There is furthermore no signal in the PIONIER *φ*<sub>C</sub> measurements, which are more precise than those available from AMBER.
We analysed the data using the standard LITPRO package<fnref n="6" /> (Tallon-Bosc et al. 2008), fitting a two-point model, with one point fixed in the centre of the map and the (x, y) coordinates of the second free to move. The overall upper limit on the flux of a secondary component is 1.7 per cent of the primary star’s flux. In order to constrain the maximum flux of a binary companion at different distances from the primary, we repeated the two-point model fit in successively wider boxes (as the current version of LITPRO does not support polar coordinates, a Cartesian approximation to annuli was used). As these flux ratios are in the H and K bands, but our estimated minimum flux ratio is in the V band, we used synthetic TLUSTY SEDs (Lanz & Hubeny 2007) to convert the H- and K-band flux ratios to V-band flux ratios. In this step, we assumed that the secondary’s $T_{\text{eff}}$ is near 20 kK, appropriate to a 5–6 $M_{\odot}$ star near the main sequence. A companion of the required brightness is ruled out beyond ~40 au.

Close binaries containing magnetic, hot stars are extremely rare, with <2 per cent of close binary systems containing a magnetic companion earlier than F0 (Alecian et al. 2015). This does not mean that a close companion can be ruled out a priori. Such a companion may be detectable via RVs. As found in Section 3, the RV curve of $\xi^1$ CMa is extremely stable, with a standard deviation of the residual RVs of ~0.5 km s⁻¹. To determine if a Be companion should have been detected, we computed RVs across a grid of models with secondary masses $1 < M_{\odot} < M_S < 10 < M_{\odot}$, semi-major axes $0.05 \text{ au} < a < 500 \text{ au}$, eccentricities $0 < e < 0.9$ and inclinations of the orbital axis from the line of sight $1^{\circ} < i < 90^{\circ}$. We then phased the JDs of the observations with the orbital periods $P_{\text{orb}}$, using a single zero-point if $P_{\text{orb}}$ was less than the time-span of the observations and multiple, evenly spaced zero-points if $P_{\text{orb}}$ was greater than this span. We then calculated the expected RV of the primary at each orbital phase, and compared the standard deviation $\sigma_{\text{orb}}$ of these RVs with the observed standard deviation $\sigma_{\text{obs}}$ (RV amplitudes were not used as, for orbits with periods longer than the time-span of observations, the full RV variation would not have been sampled). A given orbit was considered detectable if $\sigma_{\text{orb}} > 3\sigma_{\text{obs}}$. Numerical experiments with synthetic RV curves including Gaussian noise with a standard deviation of 0.5 km s⁻¹ indicate that this criterion is likely conservative: with 624 RV measurements, input periods can generally be recovered even when the semi-amplitude of the RV curve is similar to the noise level.

Each orbit was assigned a probability $P(i)$ assuming a random distribution of $i$ over $4\pi$ sr (i.e. $P(i) = 1 - \cos i$), and a flat probability distribution for $M_S$, $e$ and $a$. From this, we obtained the $1\sigma$, $2\sigma$ and $3\sigma$ upper limits on a companion’s mass as a function of $a$. This mass upper limit was then transformed into a V magnitude lower limit by interpolating along the 10 Myr isochrone in Fig. 8. Within the 40 au inner boundary of the interferometric constraints, a binary companion of sufficient brightness to host the $H\alpha$ emission can be ruled out entirely at $1\sigma$ confidence, and almost entirely at $2\sigma$.

We have proceeded under the very conservative assumption of a Be star with $H\alpha$ emission of 10 times its continuum level; however, the majority of Be stars have emission of at most a few times the continuum, and many have much less than this. Thus, the upper limits on companion mass and brightness determined above rule out all but the most exceptional of classical Be stars as a possible origin for the $H\alpha$ emission. We conclude that there is unlikely to be an undetected binary classical Be companion star that is sufficiently bright to be the source of the star’s $H\alpha$ emission.

<figure><img src="image_12.png" /><figcaption>**Figure 10.** High-resolution VLTI/AMBER data showing (clockwise from top left) flux, closure phase $\psi$, differential phases $\phi_i$ and squared visibilities $V_i^2$. Black shows the raw data. Red lines show data binned to 20 per cent of the mean error bar in the raw data. Vertical dashed lines indicate the central wavelength of the $Br\gamma$ line, and solid vertical lines denote the extent of the line.</figcaption></figure>

## 5.2 Spectro-interferometry

The formation of the line emission around the $\beta$ Cep star itself is supported by the absence of any spectro-interferometric signal across the $Br\gamma$ line. Such observations were taken with AMBER, and no signature in phase was detected on the level of $2^\circ$ (Fig. 10). This means that the photometric position as a function of wavelength remained stable on the level of 60 µas (Lachaume 2003; Štefl et al. 2009). An emission component of about 30 per cent of the strength of the continuum ($Br\gamma$ emission typically being comparable in strength to that of $H\alpha$ in Be stars) must therefore have its photocentre within 60 µas from the photocentre of the nearby continuum, as otherwise it would have produced a detectable offset of the phase signal. This not only excludes a general offset from the central star, i.e. formation around a companion, but also an extended orbiting structure, in which the blue emission would be formed at an offset opposite to the red emission.

# 6 MAGNETIC FIELD

## 6.1 Least-squares deconvolution

While $\xi^1$ CMa’s sharp spectral lines and the high S/N of the ES-PaDOnS observations mean that Zeeman signatures are visible in numerous individual spectral lines, in order to maximize the S/N we employed the usual least-squares deconvolution (LSD; Donati et al. 1997) multi-line analysis procedure. In particular, we utilized the ‘improved LSD’ (iLSD) package described by Kochukhov, Makaganiuk & Piskunov (2010). iLSD enables LSD profiles to be extracted with two complementary line masks, improving the reproduction of Stokes I and V obtained from masks limited to a select number of lines.

We used a line mask obtained from the Vienna Atomic Line Database (VALD3; Piskunov et al. 1995; Ryabchikova et al. 1997; Kupka et al. 1999, 2000) for a solar metallicity star with $T_{\text{eff}} = 27$ kK. While magnetic early-type stars are often chemically peculiar, $\xi^1$ CMa is of essentially solar composition (Niemczura & Daszyńska-Daszkiewicz 2005), albeit with a mild N enhancement (Morel et al. 2008); thus, a solar metallicity mask is appropriate. The mask was cleaned and tweaked as per the usual procedure, described in detail by Grunhut (2012), such that only metallic lines unblended

<footnote n="6">⁶ LITPRO software is available at http://www.jmmc.fr/litpro</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 11.** Top: ESPaDOnS LSD profiles: Stokes *I* (bottom), *N* (middle), Stokes *V* (top). Bottom: ESPaDOnS (black triangles) and MuSiCoS (red diamonds) LSD profiles, extracted with identical line masks, from observations acquired at similar pulsation phases (~0.44). The broader and shallower MuSiCoS profile is a consequence of the much longer exposure times of the MuSiCoS observations, the instrument's lower spectral resolution and the intrinsic variability of the star. While the MuSiCoS LSD profile is obviously noisier, the Stokes *V* signature is clearly of the opposite polarity of the ESPaDOnS observation.</figcaption></figure>

with H, He, interstellar or telluric lines remained, with the strengths of the remaining lines adjusted to match as closely as possible the observed line depths. Of the initial 578 lines in the mask, 338 remained following the cleaning/tweaking procedure. The resulting LSD profiles are shown in Fig. 11 (top panel). Note that, due to intrinsic line profile variability and the longer sub-exposure times used for the 2017 data, the line profiles of the most recent data are slightly broader than the data acquired previously.

In order to directly compare MuSiCoS and ESPaDOnS results, a second line mask was employed with all lines outside of the MuSiCoS spectral range removed, leaving 139 lines. The bottom panel of Fig. 11 shows a comparison between the highest S/N MuSiCoS LSD profile and the ESPaDOnS LSD profile with the closest pulsation phase (0.4377 versus 0.4425). The Stokes *I* profiles from the two instruments agree well, considering both the lower spectral resolution and the much longer exposure times of MuSiCoS (corresponding to ~15 per cent of a pulsation period, as compared to ~2.3 per cent of a pulsation period for ESPaDOnS). The key

point of interest is in Stokes V, which is clearly negative in the MuSiCoS LSD profile, but positive in all ESPaDOnS LSD profiles. The second MuSiCoS observation also yields a negative Stokes V profile. The reliability of these MuSiCoS observations is evaluated in detail in Appendix A.

False alarm probabilities (FAPs) were calculated inside and outside of the line profile, and classified as definite detections (DDs), marginal detections (MDs) or non-detections (NDs) according to the methodology and criteria described by Donati, Semel & Rees (1992) and Donati et al. (1997). Detection flags for Stokes V and diagnostic null *N* are given in online Table B1. All MuSiCoS observations are formal NDs. All ESPaDOnS observations are DDs (with the exception of the two discarded measurements from 2017/02/11, both of which yielded NDs due to their low S/N). However, in numerous ESPaDOnS observations, *N* also yields a DD. This phenomenon is considered in greater detail in Section 9.1.

## 6.2 Longitudinal magnetic field

The longitudinal magnetic field $\langle B_z \rangle$ was measured from the LSD profiles by taking the first-order moment of the Stokes *V* profile normalized by the EW of Stokes *I* (Mathys 1989). The same measurement using *N* yields the null measurement $\langle N_z \rangle$, which should be consistent with 0 G. In order to ensure a homogeneous analysis, the LSD profiles were first shifted by their measured RVs to zero velocity, and an identical integration range of $\pm 30 \text{ km s}^{-1}$ around line centre employed. $\langle B_z \rangle$ and $\langle N_z \rangle$ measurements are reported in online Table B1. As expected from the LSD profiles, $\langle B_z \rangle < 0$ for both MuSiCoS measurements, while $\langle B_z \rangle > 0$ for all ESPaDOnS measurements. The ESPaDOnS data are of much higher quality, with a median $\sigma_B = 6$ G, as compared to 57 G for the MuSiCoS data. In contrast to the FAPs, in which many ESPaDOnS *N* profiles yield DDs, $\langle N_z \rangle / \sigma_N$ is typically very low, with a maximum of 2.4 and a median of 0.6, i.e. statistically identical to zero.

Because $\langle B_z \rangle$ is expected to be modulated by stellar rotation, it can be used to determine the rotation period. The periodograms for $\langle B_z \rangle$ and $\langle N_z \rangle$ are shown in Fig. 12, where the latter shows the variation arising from noise. There are numerous peaks at periods of a few tens to a few hundreds of days that appear in both the $\langle N_z \rangle$ and the $\langle B_z \rangle$ periodograms. Phasing $\langle B_z \rangle$ with the periods corresponding to these peaks does not produce a coherent variation (e.g. phasing the data with the highest peak in this range, at 177 d, and fitting a first-order sinusoid, yields a reduced $\chi^2$ of 1162). The strongest peak in the $\langle B_z \rangle$ periodogram is at 5100 $\pm$ 300 d. This peak does not appear in the $\langle N_z \rangle$ periodogram. This is similar to the timespan of the ESPaDOnS data set, and so the formal uncertainty is certainly underestimated, with the 5100 d period representing a lower limit. However, the S/N of this peak is 26, above the threshold for statistical significance, indicating that the long-term variation is probably real.

Fig. 13 shows the $\langle B_z \rangle$ measurements, both individual (filled symbols) and in annual bins (open circles), as a function of time. There is an obvious long-term modulation, with $\langle B_z \rangle$ steadily declining from a peak of ~330 G in 2010 (HJD 2455200) to value of ~80 G in 2017 (HJD 2457800). The annual mean $\langle B_z \rangle$ measurements are provided in Table 2.

The MuSiCoS observations, both of which are of negative polarity, indicate (in combination with the ESPaDOnS data) that the rotational period must be longer than 5100 d. The time difference between the MuSiCoS measurements and the maximum $\langle B_z \rangle$ ESPaDOnS observations is ~3600 d. If these two epochs sample the $\langle B_z \rangle$ curve at its positive and negative extrema, then the rotational
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Periodograms for ESPaDOnS ($\langle B_z \rangle$) (top) and ($\langle N_z \rangle$) (bottom). Many of the peaks in the $\langle B_z \rangle$ periodogram also appear in the $\langle N_z \rangle$ periodogram, with the exception of the maximum amplitude peak at ~5100 d (vertical dashed line). The period determined from FORS1/2 data by Hubrig et al. (2011) is indicated by the black arrow; the period determined by Fourtune-Ravard et al. (2011) from an earlier, smaller ESPaDOnS data set is indicated by the red arrow. There is very little amplitude at either period.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Top: $\langle B_z \rangle$ versus time. The solid curve shows the least-squares sinusoidal fit to the annual mean measurements using a 30 yr period, and the dashed curves indicate the $1\sigma$ uncertainty in the amplitude and mean value of the fit. The horizontal solid line shows $\langle B_z \rangle = 0$. Vertical lines indicate the times at which $\langle B_z \rangle = 0$ (dotted), $\langle B_z \rangle = \langle B_z \rangle_{\text{max}}$ (dot-dashed) and $\langle B_z \rangle = \langle B_z \rangle_{\text{min}}$ (triple-dot-dashed). Bottom: $H\alpha$ EWs. Curves show the least-squares second-order sinusoidal fit to the annual means, which is performed independently of the fit to $\langle B_z \rangle$. Note that maximum emission corresponds to $\langle B_z \rangle_{\text{max}}$, and minimum emission corresponds to $\langle B_z \rangle = 0$. Also note that the fit predicts a second, weaker emission peak corresponding to $\langle B_z \rangle_{\text{min}}$.</figcaption></figure>

<table><caption>**Table 2.** Annual weighted mean $\langle B_z \rangle$ measurements and $H\alpha$ EWs. Months/years and HJDs correspond to the mean of the cluster of observations within each bin.</caption><thead><tr><th>Month/year</th><th>HJD -<br>2450000</th><th>⟨B<sub>z</sub>⟩<br>(G)</th><th>H&alpha; EW<br>(nm)</th></tr></thead><tbody><tr><td>02/2000</td><td>2451585</td><td>-127 &plusmn; 40</td><td>-</td></tr><tr><td>10/2000</td><td>2451822</td><td>-</td><td>0.1492 &plusmn; 0.0002</td></tr><tr><td>11/2001</td><td>2452234</td><td>-</td><td>0.1462 &plusmn; 0.0003</td></tr><tr><td>11/2002</td><td>2452580</td><td>-</td><td>0.1494 &plusmn; 0.0007</td></tr><tr><td>12/2003</td><td>2452989</td><td>-</td><td>0.1440 &plusmn; 0.0003</td></tr><tr><td>05/2004</td><td>2453156</td><td>-</td><td>0.1413 &plusmn; 0.0002</td></tr><tr><td>01/2008</td><td>2454489</td><td>314 &plusmn; 6</td><td>0.130 &plusmn; 0.001</td></tr><tr><td>01/2010</td><td>2455203</td><td>328 &plusmn; 2</td><td>0.1314 &plusmn; 0.0004</td></tr><tr><td>12/2010</td><td>2455556</td><td>318 &plusmn; 6</td><td>0.136 &plusmn; 0.001</td></tr><tr><td>02/2012</td><td>2455965</td><td>305 &plusmn; 3</td><td>0.1327 &plusmn; 0.0007</td></tr><tr><td>01/2013</td><td>2456294</td><td>281 &plusmn; 3</td><td>0.1373 &plusmn; 0.0006</td></tr><tr><td>01/2014</td><td>2456674</td><td>253 &plusmn; 5</td><td>0.1392 &plusmn; 0.0008</td></tr><tr><td>02/2017</td><td>2457797</td><td>78 &plusmn; 2</td><td>0.1507 &plusmn; 0.0003</td></tr></tbody></table>

period must be at least 7200 d (~20 yr), or a half-integer multiple if there is more than one cycle between the observed extrema. The curvature of the ESPaDOnS $\langle B_z \rangle$ suggests that $\langle B_z \rangle_{\text{max}}$ occurred at ~HJD 2455200; however, as it cannot be ruled out that $\langle B_z \rangle_{\text{min}} < -127 \pm 40$ G, it is possible that $P_{\text{rot}} > 7200$ d. Indeed, a longer period seems likely: phasing $\langle B_z \rangle$ with a 20 yr period, and fitting a least-squares first-order sinusoid (as expected for a dipolar magnetic field), produces a very poor fit as compared to longer periods. In Fig. 13, the illustrative sinusoidal fit was performed using a 30 yr period, which achieves a reasonable fit to the data (the reduced $\chi^2$ is 2.6). Note that with this fit the MuSiCoS data do not define $\langle B_z \rangle_{\text{min}}$. Longer periods can also be accommodated, however at the expense of $|\langle B_z \rangle_{\text{min}}| > |\langle B_z \rangle_{\text{max}}|$. We thus adopt $P_{\text{rot}} = 30$ yr as the most conservative option allowed by the data.

## 6.3 Comparison to previous results

There are two competing claims for rotational periods in the literature. The first, based on spectropolarimetry collected with FORS1, FORS2 and SOFIN, is approximately 2.18 d (Hubrig et al. 2011). The second period, based on a preliminary analysis of an earlier, smaller ESPaDOnS data set, is ~4.27 d (Fourtune-Ravard et al. 2011). Period analysis of the ESPaDOnS measurements using Lomb–Scargle statistics rules out both of these periods. The rotational periods provided by Hubrig et al. (2011) and Fourtune-Ravard et al. (2011) are respectively indicated with black and red arrows in the top panel of Fig. 12. There is no significant power in the periodogram at either period. The peak corresponding to the Fourtune-Ravard et al. (2011) period appears in both the $\langle B_z \rangle$ and $\langle N_z \rangle$ periodograms, suggesting it to be a consequence of noise. Phasing the data with the periods given by Hubrig et al. (2011) or Fourtune-Ravard et al. (2011) does not produce a coherent variation, as is demonstrated in Fig. 14. It is impossible for a period of the order of days to account for the systematic decline of ~200 G between the earliest data and the most recent data. We also note that the ESPaDOnS data set presented here includes two epochs with superior time sampling to that of the FORS1/2 data sets: 14 observations over 10 d in 2010, and 20 observations over 38 d in 2017, as compared to 13 FORS1 observations over 1075 d and 11 FORS2 observations over 60 d. The ESPaDOnS data set thus enables a much better probe than the FORS1/2 data set of short-term as well as long-term variability.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Top: ESPaDOnS, FORS1, FORS2 and SOFIN $\langle B_z \rangle$ measurements phased with the rotational period proposed by Hubrig et al. (2011). Bottom: ESPaDOnS $\langle B_z \rangle$ measurements phased with the rotational period proposed by Fourtune-Ravard et al. (2011). Neither of the previously proposed periods provides an acceptable phasing of the data.</figcaption></figure>

All FORS1/2 and SOFIN measurements are of positive polarity (Hubrig et al. 2006, 2011; Fossati et al. 2015b), in agreement with ESPaDOnS data. Furthermore, the magnitude, ~300 G, is similar. The two SOFIN measurements also agree well with the ESPaDOnS data. The much shorter 2.18 d period determined by Hubrig et al. (2011) is due to an apparent $\langle B_z \rangle$ variation with a semi-amplitude that is significant at $1.8\sigma$ in comparison with the median uncertainties in these measurements. Bagnulo et al. (2012) showed that systematic sources of uncertainty such as instrumental flexures must be taken into account in the evaluation of $\langle B_z \rangle$ uncertainties from FORS1 data, that the uncertainties in FORS1 $\langle B_z \rangle$ measurements should thus be about 50 per cent higher, and therefore that FORS1/2 detections are only reliable at a significance of >5$\sigma$, a much higher threshold than is satisfied by the variation the 2.18 d period is based upon. The catalogue of FORS1 $\langle B_z \rangle$ measurements published by Bagnulo et al. (2015) additionally reveals differences of up to 150 G for observations of $\xi^1$ CMa between results from different pipelines. Using the published uncertainties, the FORS1/2 measurements presented by Hubrig et al. (2011) are generally within $2\sigma$ of the sinusoidal fit to the ESPaDOnS and MuSiCoS data, and only three differ by greater than $5\sigma$. Curiously, the FORS2 data published by Hubrig et al. (2011) are systematically 80 G lower than the FORS1 measurements. We conclude that the previously published low-resolution magnetic data are not in contradiction with the long-term modulation inferred from high-resolution $\langle B_z \rangle$ measurements.

## 7 EMISSION LINES

Having rejected the possibility that $\xi^1$ CMa's H$\alpha$ emission originates in the decretion disc of a classical Be companion star (Section 5), we proceed under the assumption that it is formed within the stellar magnetosphere. In this section, we explore the star's UV and H$\alpha$ emission properties. We evaluate short-term variability with the

<figure><img src="image_17.png" /><figcaption>**Figure 15.** H$\alpha$ continuum-normalized flux (top), and residual flux (bottom) after subtraction of a photospheric model, at similar pulsation phases ($\phi \sim 0.75$) in successive epochs. The thick red line shows a model line profile disc integrated with the appropriate pulsation velocity. Vertical dashed lines indicate the integration range for EW measurement.</figcaption></figure>

stellar pulsations, as well as long-term variability consistent with the slow rotation inferred from the magnetic data.

### 7.1 H$\alpha$ emission

#### 7.1.1 Long-term modulation

Fig. 15 shows H$\alpha$ in 2000, 2004, 2010, 2014 and 2017, selected so as to share the same pulsation phase ($\phi = 0.75 \pm 0.03$). A synthetic line profile calculated using the physical parameters obtained in Section 4 is overplotted as a thick line. Comparison of observed to synthetic line profiles demonstrates that emission is present at all epochs, although it is substantially weaker in the earlier CORALIE data. The agreement of H$\beta$ with the model is reasonable at all epochs, although there is also weak emission present in the line core (Fig. 4). Close analysis of the line core of H$\beta$ shows evidence of evolution, with the same trend as in H$\alpha$, although this evolution is of a much lower amplitude.

H$\alpha$ EWs were measured with an integration range of $\pm 0.4$ nm of the rest wavelength in order to include only the region with emission (Fig. 15). The spectra were first shifted to a rest velocity of 0 km s$^{-1}$ by subtracting the RVs measured in Section 3. The periodogram for these measurements is shown in the middle panel of Fig. 16. It shows a peak at 6700 d or 18 yr. The S/N of this peak is 29, indicating that it is statistically significant. This is consistent with the very long timescale of variation inferred from the ESPaDOnS $\langle B_z \rangle$ periodogram (Fig. 12), and with the minimum 7200 d period inferred from the positive and negative $\langle B_z \rangle$ extrema in the ESPaDOnS and MuSiCoS data.

The bottom panel of Fig. 13 shows the H$\alpha$ EWs as a function of time. Maximum emission and maximum $\langle B_z \rangle$ occur at approximately the same phase, as do the $\langle B_z \rangle$ minimum and the minimum emission strength. The variation in EW is seen more clearly in the annual mean EWs (open circles), which are tabulated in Table 2. The relative phasing of $\langle B_z \rangle$ and EW is consistent with a corotating magnetosphere, and supports the accuracy of the adopted period. Since $\langle B_z \rangle$ is both positive and negative, two local emission maxima are expected at the positive and negative extrema of the $\langle B_z \rangle$ curve, as at these rotational phases the magnetospheric plasma is seen
<figure><img src="image_18.png" /><figcaption>**Figure 16.** Periodogram for $H\alpha$ EWs (top), and after pre-whitening with the highest amplitude period (bottom). The highest amplitude period is indicated by the dashed line. The red arrow indicates the pulsation period determined from the RVs.</figcaption></figure>

closest to face on. The solid curve in the bottom panel of Fig. 13 shows a second-order sinusoidal fit to the annual mean EWs, which indeed yields two local maxima, the strongest corresponding to $\langle B_z \rangle = \langle B_z \rangle_{\text{max}}$, and the second maximum predicted to occur at $\langle B_z \rangle = \langle B_z \rangle_{\text{min}}$. Note that, due to the incomplete phase coverage (less than half of a rotational cycle), in the event that the $H\alpha$ variation is indeed a double wave a Lomb–Scargle periodogram should show maximum power at close to half of the rotational period. The periodogram peak at 18 yr would then indicate $P_{\text{rot}} = 36$ yr, consistent with the rotational period inferred from $\langle B_z \rangle$.

### 7.1.2 Short-term modulation

As a first step to investigating whether $H\alpha$ is affected by pulsation, residual EWs were obtained by subtracting the least-squares second-order sinusoidal fit to the annual mean EWs. The bottom panel of Fig. 16 shows the period spectrum of the residual EWs. The strongest peak is at the stellar pulsation period, with an S/N of 12. After pre-whitening the EWs with both the rotational and pulsational frequencies, the S/N of the highest peak in the periodogram is close to 4, suggesting that all significant variation is accounted for by these two frequencies.

While the semi-amplitude of the residual EW variation is similar to the median error bar, binning the residual EWs by pulsation phase does not change the semi-amplitude, but increases the significance of the variation to ~8σ with respect to the mean error bar. The phase-binned residual EWs are shown phased with the pulsation period in the top panel of Fig. 17. There is a coherent variation of the phase-binned residual EWs with the pulsation period, with a semi-amplitude of approximately 0.003 nm.

One obvious candidate mechanism for producing this pulsational modulation is the change in the EW of the underlying photospheric profile due to the changing $T_{\text{eff}}$, as explored in Section 4.5. To investigate this hypothesis, we calculated synthetic spectra via linear interpolation between the grid of TLUSTY BSTAR2006 models (Lanz

<figure><img src="image_19.png" /><figcaption>**Figure 17.** Top: residual $H\alpha$ EWs, after removal of the sinusoidal variation shown in the bottom panel of Fig. 13, phased with the pulsation period and binned with pulsation phase. The thick solid black line indicates the EW variation expected from a photospheric model calculated using the ±300 K $T_{\text{eff}}$ variation determined in Section 4.5, normalized by subtracting the mean EW of the model. The amplitude of the residual EW variation is three times larger than predicted by the photospheric model, and is about 0.5 cycles out of phase. Bottom: $H\alpha$ $V/R$ variation. Symbols and colours indicate the year of observation, the last two digits of which are given in the legend. The amplitude of the variation is much higher than predicted by the photospheric model, which is essentially flat on this scale and therefore is not shown. The amplitude furthermore changes between epochs, which is not expected for photospheric pulsations given the high stability of pulsation properties in other lines.</figcaption></figure>

& Hubeny 2007). We used the physical parameters from Section 4, including the ±300 K $T_{\text{eff}}$ variation found in Section 4.5. Deformation of the line profile due to pulsation was accounted for using the RV curve from Section 3, and modelling the (assumed radial) pulsations as described in Section 4.1. Spectra were calculated at 20 pulsation phases, and the EW was measured at each phase using the same integration range as used for the observed data, and after moving the synthetic spectra to zero RV. The black line in the top panel of Fig. 17 shows the resulting variation, normalized by subtracting the mean EW across all models. The semi-amplitude of the EW variation expected due to changes in the photospheric profile due to pulsation is about 0.001 nm, much smaller than observed. Furthermore, it is out of phase with the observed variation by about 0.5 pulsation cycles. Note that the maximum and minimum EWs of the predicted variation correspond to the minimum and maximum of the $T_{\text{eff}}$ variation (Fig. 7), as expected for H lines, which grow weaker with increasing $T_{\text{eff}}$ in this $T_{\text{eff}}$ range.

Radial pulsation introduces asymmetry into the line profile, which can be quantified by $V/R$, defined as the ratio of the EW in the blue half to the EW in the red half of the line. We measured $V/R$ from EWs calculated from -0.4 nm to line centre (defined at the laboratory rest wavelength shifted by the RV measured from metallic lines in Section 3), and line centre to +0.4 nm. The bottom panel of Fig. 17 shows the $V/R$ variation. The semi-amplitude of the variation is about 0.25. This is much higher than the semi-amplitude of 0.0007 predicted by synthetic spectra calculated using a radially pulsating photospheric model; since this is essentially flat on the scale of Fig. 17, the photospheric variation is not shown. The low
level of line asymmetry in the synthetic spectra is a consequence of the large Doppler broadening of the Hα line, approximately 30 km s⁻¹ or about twice the semi-amplitude of the RV curve. For narrower lines, e.g. the C II 656.3 nm line for which the Doppler velocity is similar to the RV semi-amplitude, the predicted and observed V/R variations are in reasonable agreement.

The semi-amplitude of the Hα V/R variation increases steadily from 0.16 in 2000–2001, to 0.20 in 2002–2004, to 0.35 in 2008–2010, and then declines to 0.24 in 2012–2014 and 0.13 in 2017. This is the same pattern as the change in total emission strength; thus, the amplitude of the V/R variation correlates to the total emission strength and, therefore, to ⟨B_z⟩.

For a more detailed view of the pulsational modulation of Hα, we calculated dynamic spectra phased with the pulsation period. These are shown in Fig. 18, using the synthetic line profile from Fig. 15 as a reference spectrum. Individual line profiles were moved to 0 km s⁻¹ by subtracting their RVs, and then binned by pulsation phase using phase bins of 0.05 cycles. The two panels of Fig. 18 show dynamic spectra in the epochs of minimum emission (2000–2001) and maximum emission (2008–2010). The emission line morphology and pattern of variability are essentially identical in other epochs. In both cases, the emission peaks near the centre of the line. In the CORALIE data, the Hα emission peak is ~20 per cent of the continuum. This rises to about 28 per cent of the continuum in the ESPaDOnS data. The emission peak anticorrelates slightly with the RV, as shown by the overplotted white lines. The strong red and blue emission variability revealed by the V/R variation in the bottom panel of Fig. 17 is due to secondary emission peaks that occur at phases 0.0 and 0.5, respectively blue- and redshifted with respect to the line centre. As with the central emission peak, these are stronger in the ESPaDOnS data, reflecting the change in V/R amplitude over time. Note that, in the data acquired at earlier epochs, these secondary emission bumps are apparently separated from the main emission peak, while in later epochs they are connected (although this depends on the choice of the reference spectrum).

Both the emission strength (as measured by the residual EWs after pre-whitening with the second-order fit in Fig. 13) and the line asymmetry (as quantified by V/R) vary coherently with the pulsation phase. Synthetic photospheric spectra calculated using a radially pulsating model are unable to reproduce the residual EW variation, which is both three times larger than predicted, and 0.5 cycles out of phase. The amplitude of the observed V/R variation is about 300 times larger than predicted by the model, which does not reproduce the prominent blue- and redshifted secondary emission bumps. The amplitude of V/R is furthermore variable with time, increasing and decreasing in strength in the same fashion as the total emission strength. As the star’s pulsation amplitude is extremely regular, a change in the amplitude of V/R with epoch cannot be explained by photospheric pulsation. These discrepancies between model and observation suggest either that the origin of the pulsational modulation of Hα is not photospheric or that it cannot be explained due to Tₐff variation alone. Exotic processes, such as temperature inversions related to shockwaves produced by the star’s supersonic pulsations, may be one explanation. Alternatively, the origin of the pulsational modulation may reside within the magnetosphere.

## 7.2 UV emission lines

Emission is present in four of the UV doublets often used to diagnose the wind properties of early-type stars: N v 1239, 1243 Å; Si iv 1394, 1403 Å; C iv 1548, 1551 Å; and Al III 1854, 1863 Å. The

<figure><img src="image_20.png" /><figcaption>**Figure 18.** Hα dynamic spectra in the epochs of minimum (2000–2001) and maximum (2008–2010) emission. The synthetic spectrum shown in Fig. 15 was used as a reference spectrum. Observed spectra were shifted to 0 km s⁻¹ by subtracting their RVs, and then co-added in 20 phase bins. The white curves indicate the RV curve due to pulsations. In all epochs, the emission profile has three components: a central peak near 0 km s⁻¹, which shifts slightly about line centre, but stays nearly constant in strength; a blueshifted peak, which is at maximum strength at phase 0 and disappears by phase 0.5; and a redshifted peak, which disappears at phase 0 and is at a maximum at phase 0.5. All three components increase in strength from 2000 to 2010. Note the greater amplitude of the variation in the secondary emission lobes in the 2008–2010 data.</figcaption></figure>
<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /><figcaption>**Figure 19.** Comparison between the mean N v (top) and Al III (bottom) emission lines of $\xi^1$ CMa (solid, black) and those of $\beta$ Cep at maximum emission (red, dashed) and maximum absorption (blue, dot-dashed). Vertical dashed lines indicate rest wavelengths. $\xi^1$ CMa's emission lines are similar to those of $\beta$ Cep at maximum emission, albeit somewhat stronger. The lower panels show temporal variance spectra (TVS, see the text), calculated from all spectra, which characterize the variability of the spectra normalized to the continuum. The horizontal line indicates the continuum variation, in comparison to which there is no statistically significant variability inside the line profile.</figcaption></figure>

emission profiles of these doublets are similar to those of the magnetic $\beta$ Cep pulsator $\beta$ Cep at maximum emission. Schnerr et al. (2008) performed this comparison for the C IV doublet. Fig. 19 compares the mean line profiles for $\xi^1$ CMa's N v and Al III lines to those of $\beta$ Cep, and demonstrates that these lines are also similar to those of $\beta$ Cep at maximum emission. Such emission is unique to magnetic stars, and as an indirect diagnostic of stellar magnetism has historically prompted the search for magnetic fields in such stars (e.g. Henrichs 1993; Henrichs et al. 2013; indeed, the detection of $\xi^1$ CMa's UV emission motivated the collection of the MuSiCoS data).

The presence of emission in the N v line is particularly interesting. This line is not seen in normal early B-type stars (Walborn, Parker & Nichols 1995), nor is it present in most normal late O-type stars (e.g. Brandt et al. 1998). In their study of IUE data for four magnetic B-type stars, including $\beta$ Cep, Smith & Groote (2001) concluded that this doublet must be formed at ~30 kK, somewhat higher than the photospheric $T_{\text{eff}}$. The presence of N v lines in the UV spectra of relatively cool stars is thought to be a consequence of Auger ionization due to the presence of X-rays (Cassinelli & Olson 1979). Since $\xi^1$ CMa has a strong magnetic field, it is expected to be

overluminous in X-rays due to magnetically confined wind shocks, and is indeed observed to be overluminous in X-rays to the degree predicted by models (Oskinova et al. 2011; Nazé et al. 2014).

$\beta$ Cep's wind lines show clear variability synchronized with its rotation period (Henrichs et al. 2013). In contrast, $\xi^1$ CMa's wind lines show only a very low level of variability, more similar to that seen in a normal (magnetically unconfined) stellar wind (Schnerr et al. 2008). The bottom panels of Fig. 19 show temporal variance spectra (TVS), which compare the variance within spectral lines to the variance in the continuum (Fullerton, Gies & Bolton 1996). To minimize variation due to pulsation, spectra were moved to 0 km s$^{-1}$ by subtracting the RV computed on the basis of the RV curve and ephemeris determined in Section 3. With the uncertainty in $P_0$ and in $\dot{P}$ (Section 3), the uncertainty in pulsation phase is about 0.028, corresponding to a maximum uncertainty in RV of 3 km s$^{-1}$, less than the ~7.8 km s$^{-1}$ velocity pixel of the IUE data. RV correction reduces the TVS pseudo-continuum level by a factor of about 2–3.

Comparison of EW measurements of these lines to EWs of synthetic spectra, using the same radially pulsating model described above in Section 7.1, yielded ambiguous results, due to the small number of high-dispersion IUE observations (13 spectra), and formal uncertainties similar to the maximum level of variability, ~0.02 nm. This low level of variability is consistent with the weak variability of the residual $H\alpha$ EWs.

The lack of variability in $\xi^1$ CMa's wind lines is consistent with a rotational pole aligned with the line of sight, an aligned dipole, a long rotation period or some combination of these. The similarity to $\beta$ Cep's emission lines at maximum emission suggests that the magnetic pole was close to being aligned with the line of sight when the UV data were acquired. The long-term modulations of both $H\alpha$ and $\langle B_z \rangle$ favour an oblique dipole with a long rotational period, in which case the UV data should have been acquired at a rotational phase corresponding to one of the extrema of the $\langle B_z \rangle$ curve. Phasing the UV data with the same 30 yr rotational period as in Fig. 13 yields a phase close to 1.0, i.e. they would indeed have been obtained close to magnetic maximum. Assuming that the data must have been acquired near an extremum of the $\langle B_z \rangle$ curve (i.e. at a phase close to 0.5 or 1.0) would require a period of 10, 15, 20, 30 or 60 yr, of which only 30 and 60 yr are not excluded by the magnetic data. This assumes that the UV data were, in fact, acquired close to an extremum. Given that the UV emission is only slightly stronger than that of $\beta$ Cep, which has a weaker magnetic field and no $H\alpha$ emission, it may be the case that the true maximum UV emission strength is significantly in excess of observations, in which case the UV data cannot be used to infer $P_{\text{rot}}$.

## 8 MAGNETIC AND MAGNETOSPHERIC PARAMETERS

In this section, the ~30 yr rotation period inferred from magnetic and spectroscopic data is used with the star's physical parameters and $\langle B_z \rangle$ measurements to establish constraints on the properties of $\xi^1$ CMa's surface magnetic field, circumstellar magnetosphere and spin-down time-scales, using the self-consistent Monte Carlo method described by Shultz (2016). These results are summarized in Table 3.

The $\langle B_z \rangle$ variation of a rotating star with a dipolar magnetic field can be reproduced with a three-parameter model: the inclination angle $i$ between the rotational axis and the line of sight, the obliquity angle $\beta$ between the magnetic and rotational axes, and the polar strength of the magnetic dipole at the stellar surface $B_d$. The two angular parameters are related via
<table><caption>**Table 3.** Magnetic, magnetospheric and rotational properties. The final column gives the section in which the parameter is derived.</caption><thead><tr><th>Parameter</th><th>Value</th><th>Ref.</th></tr></thead><tbody><tr><td>P<sub>rot</sub> (d)</td><td>&gt;30 yr</td><td>6.1</td></tr><tr><td>r</td><td>-0.78 ± 0.02</td><td>8.1</td></tr><tr><td>β (°)</td><td>83<sup>+5</sup><sub>-12</sub></td><td>8.1</td></tr><tr><td>ε</td><td>0.36 ± 0.02</td><td>4.2</td></tr><tr><td>B<sub>d</sub> (kG)</td><td>&gt;1.1</td><td>8.1</td></tr><tr><td>W</td><td>&lt;7 × 10<sup>-5</sup></td><td>8.2</td></tr><tr><td>R<sub>K</sub> (R<sub>*</sub>)</td><td>&gt;580</td><td>8.2</td></tr><tr><td>log M(M<sub>☉</sub> yr<sup>-1</sup>)</td><td>-8.0 ± 0.1</td><td>8.2, Vink et al. (2001)</td></tr><tr><td>v<sub>∞</sub> (km s<sup>-1</sup>)</td><td>2070 ± 60</td><td>8.2, Vink et al. (2001)</td></tr><tr><td>η<sub>*</sub></td><td>&gt;340</td><td>8.2</td></tr><tr><td>R<sub>A</sub> (R<sub>*</sub>)</td><td>&gt;4.6</td><td>8.2</td></tr><tr><td>τ<sub>J</sub></td><td>4 ± 1</td><td>8.3</td></tr><tr><td>t<sub>S,max</sub></td><td>42 ± 12</td><td>8.3</td></tr></tbody></table>

<figure><img src="image_23.png" /><figcaption>**Figure 20.** Oblique rotator model parameters determined via a Monte Carlo approach. Red, blue and pink show the 1σ, 2σ and 3σ contours of the grid. The bottom panels show PDFs for $B_d$, $β$ and $i$, where $i$ was drawn from a random distribution over a solid angle of $4π$.</figcaption></figure>

tan $β = (1 - r)/(1 + r)$ cot $i$, where $r$ is the ratio defined by Preston (1967) as $r = (|B_0| - B_1)/(|B_0| + B_1)$, with $B_0$ and $B_1$ the mean and semi-amplitude of the sinusoidal fit to $\langle B_z \rangle$ when phased with $P_{rot}$. We determined $r = -0.78 \pm 0.02$ from the sinusoidal fit to the annual mean $\langle B_z \rangle$ measurements shown in the top panel of Fig. 13. In consequence, these results rely on the 30 yr period.

The degeneracy between $i$ and $β$ is usually broken by determining $i$ independently, e.g. from $v\sin i$, $P_{rot}$ and $R_*$: Using this method, Hubrig et al. (2011) found $i \sim 3°$; however, this was based on their rotation period of 2.18 d, which was shown in Section 6.2 to be incorrect. Indeed, the equatorial rotational velocity implied by a 30 yr period, $v_{eq} < 0.04$ km s$^{-1}$, is much less than the upper limit on $v\sin i < 8$ km s$^{-1}$ found in Section 4.1 (and indeed, much less than the 1.8 km s$^{-1}$ velocity resolution of ESPaDOnS data). Since $i$ cannot be constrained, we applied a probabilistic prior, requiring that $i$ be drawn from a random distribution over $4π$ sr such that $P(i) = 1 - \cos i$ (e.g. Landstreet & Mathys 2000). The corresponding probability density function (PDF) is shown in the bottom-right panel of Fig. 20. The bottom-middle panel shows the resulting PDF for $β$, which peaks at $83°$. Despite the inverse relationship of $i$ and

$β$ obtained for $r = -0.78$ (middle panel of Fig. 20), and the bias towards large $i$ in the prior, large $β$ are favoured overall. We note that we obtain a similar $β$ to that given by Hubrig et al. (2011); however, this is simply a coincidence as their value was obtained with an incorrect rotation period and a very different value of the Preston $r$ parameter (0.56, as compared to -0.78).

For each $(i, β)$ pair, $B_d$ was determined using equation 1 from Preston (1967). This also requires knowledge of the limb darkening coefficient $ε$, which we obtained from the tables calculated by Díaz-Cordovés et al. (1995) as $ε = 0.36 \pm 0.02$ (Section 4.1). The resulting PDF peaks sharply at the minimum value permitted by Preston's equation, $B_d = 1.1$ kG. As demonstrated in the middle and upper-left panels of Fig. 20, $B_d < 2$ kG over the range $20° < i < 85°$, with a low-probability tail extending out to several kG for very large and very small $i$ and $β$. This is much lower than the value of ~5.3 kG given by Hubrig et al. (2011), which arose from their very small $i$.

If instead the ESPaDOnS and MuSiCoS measurements define the extrema of $\langle B_z \rangle$, we obtain $r = -0.33 \pm 0.02$. This does not affect $B_d$, although a slightly smaller $β \sim 75°$ is favoured.

Magnetic wind confinement is governed by the balance of kinetic energy density in the radiative wind to the magnetic energy density. This is expressed by the dimensionless wind magnetic confinement parameter $η_*$ (ud-Doula & Owocki 2002). If $η_* > 1$, the star possesses a magnetosphere.

Oskinova et al. (2011) analysed IUE observations of $\xi^1$ CMa using the PoWR code in order to determine the wind parameters, obtaining $log(M/M_{☉} yr^{-1}) = -10$ and $v_{∞} = 700$ km s$^{-1}$. However, magnetic confinement reduces the net mass-loss rate and, more seriously, strongly affects line diagnostics due to the departure from spherical symmetry in the circumstellar environment. Indeed, Oskinova et al. were unable to achieve a simultaneous fit to the Si IV and C IV doublets. Comparison of magnetohydrodynamic (MHD) simulations and spherically symmetric models to the magnetospheric emission of Of?p stars has demonstrated that MHD simulations yield a superior fit, and require mass-loss rates comparable to those obtained from the Vink, de Koter & Lamers (2001) recipe, whereas spherically symmetric models in general require much lower mass-loss rates to achieve a relatively poor match to the observations (Grunhut et al. 2012b; Marcolino et al. 2012, 2013). In addition to this, magnetic wind confinement is not expected to modify the surface mass flux, which is the relevant quantity in calculating the strength of magnetic confinement (ud-Doula & Owocki 2002). Thus, $η_*$ should be determined using the wind parameters as they would be in the absence of a magnetic field, rather than those measured via spectral modelling (e.g. Petit et al. 2013). Using the mass-loss recipe of Vink et al. (2001) with $T_{eff}$, log L and $M_*$ from Table 1, assuming the metallicity $Z/Z_{☉} = 1$, and determining the wind terminal velocity $v_{∞} = 2070 \pm 60$ km s$^{-1}$ via scaling the star's escape velocity by 2.6 as suggested by Vink et al., yields a mass-loss rate $log(M/M_{☉} yr^{-1}) = -8.0 \pm 0.1$. From equation 7 of ud-Doula & Owocki (2002), $η_* > 340$, so the wind is magnetically confined.

The physical extent of the magnetosphere is given by the Alfvén radius $R_A$, defined as the maximum extent of closed magnetic field lines in the circumstellar environment. $R_A$ can be calculated heuristically from $η_*$ (equation 7 in ud-Doula, Owocki & Townsend 2008): we find $R_A > 4.6R_*$: Unless $i$ is particularly large or small, in which case $B_d$ is significantly in excess of 2 kG, the magnetic and magnetospheric parameters will be close to the derived lower limits. $R_A$ is almost certainly below $20 R_*$.

$\xi^1$ CMa has the highest X-ray luminosity and the hardest X-ray spectrum of any of the magnetic $β$ Cep stars (Cassinelli et al. 1994; Oskinova et al. 2011). The X-ray luminosity, corrected
for interstellar absorption corresponding to $E(B-V) = 0.015$, is $2.4 \times 10^{31} \text{ erg s}^{-1}$ (Nazé et al. 2014). Using 2D MHD simulations, ud-Doula et al. (2014) calculated X-ray emission from magnetically confined wind shocks and developed an X-ray Analytic Dynamical Magnetosphere (XADM) scaling for log $L_X$ with $B_d$, $R_*$, $\dot{M}$ and $v_\infty$. Comparison to available X-ray data for magnetic early-type stars has indicated that log $L_X$ predicted by XADM should be scaled by ~5–20 per cent to match the observed log $L_X$ (Nazé et al. 2014). This efficiency factor accounts for dynamical infall of the plasma. XADM successfully predicts the X-ray luminosity of magnetic, hot stars across three decades in log L, two decades in $B_d$ and five decades in log $L_X$ (Nazé et al. 2014), i.e. the range of the model’s successful application is much larger than the uncertainty introduced by the efficiency factor. Assuming an efficiency of 10 per cent, $B_d = 1.1$ kG, and taking into account the uncertainties in the stellar parameters in Table 1, XADM predicts log $L_X = 31.5 \pm 0.1$, in excellent agreement with the observed X-ray luminosity log $L_X = 31.47$. Adopting the higher value of $B_d = 5.3$ kG suggested by Hubrig et al. (2011) yields log $L_X = 31.7$, slightly higher than observed (although this can be reconciled by lowering the efficiency factor to 5 per cent). If the lower $\dot{M}$ and $v_\infty$ determined by Oskinova et al. (2011) are used instead, XADM predicts log $L_X = 28.9$ with an efficiency factor of 100 per cent, 2.6 dex lower than observed. Since the efficiency factor can only lower log $L_X$, it is impossible for the XADM model to match $\xi^1$ CMa’s observed X-ray luminosity with a mass-loss rate significantly lower than the Vink et al. (2001) prediction.

Petit et al. (2013) divided magnetic, massive stars into two classes, those with dynamical magnetospheres (DMs) only and those also possessing centrifugal magnetospheres (CMs). CMs appear when $R_K < R_A$, where $R_K$ is the Kepler radius, defined as the radius at which the centrifugal force due to corotation compensates for the gravitational force (Townsend & Owocki 2005; ud-Doula et al. 2008). Solving for $R_K$ using equation 12 from Townsend & Owocki (2005) with $M_*$ from Table 1 and $P_{\text{rot}} > 30$ yr yields $R_K > 580R_*$. As $R_K \gg R_A$, the magnetosphere does not include a CM. The Kepler radius is related to the dimensionless rotation parameter $W \equiv v_{\text{eq}}/v_{\text{orb}} = R_K^{-3/2}$, where $v_{\text{orb}}$ is the velocity required to maintain a Keplerian orbit at the stellar surface (ud-Doula et al. 2008). Critical rotation corresponds to $W = 1$, and no rotation to $W = 0$. For $\xi^1$ CMa, $W < 7 \times 10^{-5}$.

Magnetic wind confinement leads to rapid spin-down due to angular momentum loss via the extended moment arm of the magnetized wind (Weber & Davis 1967; ud-Doula, Owocki & Townsend 2009). The rotation period will decrease exponentially, with a characteristic angular momentum loss time-scale $\tau_J$ of (ud-Doula et al. 2009)

$$ \tau_J = \frac{3}{2} f \tau_M \left( \frac{R_*}{R_A} \right)^2, \qquad (2) $$

where $\tau_M \equiv M_*/\dot{M}$ is the mass-loss time-scale and f is the moment of inertia factor, which can be evaluated from the star’s radius of gyration $r_{\text{gyr}}$ as $f = r_{\text{gyr}}^2$. Consulting the internal structure models calculated by Claret (2004), $f \simeq 0.06$ for a star of $\xi^1$ CMa’s mass and age. Solving equation (2) then yields $\tau_J = 4 \pm 1$ Myr. The rotation parameter at a time t after the birth of the star is

$$ W(t) = W_0 e^{-t/\tau_J}, \qquad (3) $$

where $W_0$ is the initial rotation parameter. Assuming $W_0 = 1$ yields the maximum spin-down age $t_{S,\text{max}} = 42 \pm 12$ Myr (Petit et al. 2013). This is about four times longer than the age inferred from evolutionary tracks. Solving equation (3) for $W_0$ yields

$W_0 < 0.02$. Thus, either magnetic braking must have been much more rapid than predicted or almost all of the star’s angular momentum loss must have occurred before it began its main-sequence evolution, i.e. the star was already a slow rotator at the zero-age main sequence (ZAMS).

Given the important role played by the mass-loss rate in determining $R_A$ and $\tau_J$, it is of interest to explore the sensitivity of these results to different mass-loss prescriptions. Muijres et al. (2012) found that they were unable to reproduce the observed mass-loss rates of stars with log L < 5.2, suggesting that $\dot{M}$ may be lower than predicted by the Vink et al. (2001) recipe. We first note that, since $\tau_J$ increases with decreasing $\dot{M}$, a lower mass-loss rate cannot resolve the discrepancy between $t_{S,\text{max}}$ and the age inferred from the HRD. Secondly, we note that satisfying the condition for a CM ($R_A > R_K$) requires log ($\dot{M}/\text{M}_\odot \text{yr}^{-1}$) $\le -16.5$. If the lower $\dot{M}$ and $v_\infty$ found by Oskinova et al. (2011) are used to calculate the magnetospheric and spin-down parameters, none of the above conclusions are fundamentally changed ($R_A \sim 14R_* \ll R_K$, and $t_{S,\text{max}} = 250_{-100}^{+80}$ Myr). Lucy (2010) computed mass-loss rates for late O-type stars in the weak-wind domain using the theory of reversing layers (Lucy 2007), and for a star with $\xi^1$ CMa’s $T_{\text{eff}}$ and log g predict log ($\dot{M}/\text{M}_\odot \text{yr}^{-1}$) $\sim -8.8$. Krtička (2014) provided an alternate calculation for $\dot{M}$, intended specifically for B-type stars, and for $\xi^1$ CMa predicted log ($\dot{M}/\text{M}_\odot \text{yr}^{-1}$) $= -8.9$ (see table 5 in Krtička 2014), similar to the Lucy (2010) prediction. Thus, while there are systematic differences between theoretical mass-loss rates, these are much smaller than would be required to change the basic conclusions that the star lacks a CM and has a spin-down age much less than its evolutionary age.

# 9 DISCUSSION

## 9.1 Impact of pulsation on the magnetometry

While $\langle B_z \rangle$ exhibits a clear long-term modulation, it also shows evidence for short-term variability. In particular in 2010 (HJD~2455200), the most densely time-sampled epoch, there is substantial apparent scatter in the measurements: the mean error bar is 5 G, but the standard deviation of $\langle B_z \rangle$ in this epoch is 15 G.

We calculated residual $\langle B_z \rangle$ measurements by subtracting the mean $\langle B_z \rangle$ at each annual epoch in order to remove the long-term trend. These are shown phased with the pulsation period in Fig. 21. The reduced $\chi^2$ of a least-squares sinusoidal fit to all residual $\langle B_z \rangle$ measurements is 1.2, as compared to 2.8 for the null hypothesis of no variation. In 2010, there is an apparently coherent variation

<figure><img src="image_24.png" /><figcaption>Figure 21. Residual ESPaDOnS $\langle B_z \rangle$ measurements, after subtraction of the mean $\langle B_z \rangle$ at each epoch, phased with the pulsation period. The legend indicates the year in which each data set was acquired. Only epochs with at least four observations are shown.</figcaption></figure>
with a semi-amplitude of 16 ± 2 G, consistent with the standard deviation in $\langle B_z \rangle$ of 15 G, and three times larger than the median error bar of 5 G. The reduced $\chi^2$ of the least-squares sinusoidal fit to the 2010 residual $\langle B_z \rangle$ measurements is 0.9, as compared to 6.1 for the null hypothesis. In 2017, the epoch with the largest number of observations and the most systematic sampling of the phase curve, the same fit yields a semi-amplitude of 7 ± 2 G, which is identical to both the standard deviation of, and the median uncertainty in, $\langle B_z \rangle$ during this epoch. The reduced $\chi^2$ of a sinusoidal fit and the null hypothesis to the 2017 data is close to 1 in both cases, indicating that in 2017 there is no evidence for a coherent variation with pulsation phase. In 2013, the only other epoch with sufficient observations to constrain a sinusoidal fit, the semi-amplitude is 10 ± 5 G, the reduced $\chi^2$ of the fit is 1.3 and the reduced $\chi^2$ of the null hypothesis is 2.3, i.e. 2013 yields results intermediate between 2010 and 2017. This suggests that the observed variation of $\langle B_z \rangle$ may be explained as the superposition of a short-term variation with the stellar pulsation period on top of the long-term variation due to rotation, with the short-term variation declining in amplitude with $\langle B_z \rangle$.

We next turn our attention to the LSD profiles. It is clear from the numerous 'definite' detections in LSD *N* profiles that pulsation has affected our data to some degree. Signatures in *N* profiles are sometimes seen in stars in which the stellar lines vary during acquisition of the polarized sub-exposures used to construct the final Stokes V spectrum. Following the analysis performed by Neiner et al. (2012), we modelled the *N* profiles by creating disc-integrated Stokes I and V profiles at the RVs corresponding to the pulsation phases of each sub-exposure, and then combining the individual model profiles in such a way as to simulate the double-ratio method's combination of I ± V beams to yield Stokes I, Stokes V and *N* (Donati et al. 1997). Each model profile was created with $v\sin i = 0$ km s$^{-1}$, a macroturbulent velocity of 8 km s$^{-1}$, the measured RVs, a projection factor of 1.45, and a constant dipolar magnetic field with $i = \beta = 90^\circ$ and $B_d = 1.1$ kG, with the positive magnetic pole at the centre of the stellar disc. The model profiles were normalized to the observed Stokes I EW.

Illustrative results are shown for Stokes I, *N* and Stokes V in Fig. 22. The model reproduces the variation in *N* reasonably well. We conclude from this that the *N* profile signatures are principally a consequence of line profile variability between sub-exposures.

There is no change in $\langle B_z \rangle$ measured from model LSD profiles created with or without introducing RV shifts between sub-exposures. Furthermore, $\langle N_z \rangle$ remains null, as expected. This result is in agreement with the observational results of Neiner et al. (2012), who found that while shifting sub-exposures by their pulsation velocities greatly reduced the *N* profile signatures, $\langle B_z \rangle$ remained unchanged within error bars.

If the modulation of $\langle B_z \rangle$ with pulsation phase is not a consequence of an instrumental effect, it might be due to a real change in the intrinsic surface magnetic field strength of the star. An obvious mechanism that might produce such an effect is conservation of magnetic flux throughout radial pulsation cycles. Magnetic flux is conserved as $1/R_*^2$, and the semi-amplitude of the change in stellar radius is $\Delta r/R_* \sim 1-1.5$ per cent; thus, $\langle B_z \rangle$ could vary by about ±2–3 per cent or ±7–10 G at $\langle B_z \rangle_{max} = 328$ G, which is close to the observed semi-amplitude in 2010 of 16 ± 2 G. In this scenario, maximum $\langle B_z \rangle$ should occur at phase 0.25, when the star is at its most contracted, and minimum $\langle B_z \rangle$ at phase 0.75, when the star is at its most extended; this is indeed what is observed. Thus, a modulation of $\langle B_z \rangle$ arising due to an actual change in the surface magnetic field strength with the pulsation period is consistent with both the magnitude and the phasing of the observed effect. The

<figure><img src="image_25.png" /><figcaption>**Figure 22.** Stokes I (left), null *N* (middle) and Stokes V (right) LSD profiles. Observed profiles are shown by black circles, and synthetic profiles by red lines.</figcaption></figure>

apparent proportionality of the semi-amplitude of the pulsational modulation of $\langle B_z \rangle$ to the magnitude of $\langle B_z \rangle$ is also consistent with this hypothesis, as a change in $\langle B_z \rangle$ due to a change in $B_d$ should be largest when there is the least amount of magnetic flux cancellation across the stellar disc (i.e. at magnetic maximum).

While our model does not predict a change in $\langle B_z \rangle$ due to RV shifts between sub-exposures, we cannot entirely rule out the possibility that the modulation shown in Fig. 21 is due to line profile variability. Whatever its origin, the influence of pulsation on $\langle B_z \rangle$ is much smaller than the long-term modulation of $\langle B_z \rangle$; thus, the fundamental conclusion that $P_{\text{rot}} > 30$ yr should not be affected by pulsation.

## 9.2 Does extremely slow rotation in an evolved magnetic B-type star make sense?

Due to magnetic braking, evolved magnetic hot stars are expected to rotate more slowly than either similar, younger magnetic stars or non-magnetic stars of comparable age and mass. However, as was shown in Section 8, the spin-down age inferred for $\xi^1$ CMa is difficult to reconcile with the age inferred from its position on the HRD. The Ekström et al. (2012) evolutionary models imply an age of 11 ± 0.7 Myr, while its spin-down time-scale is $\tau_J = 4 \pm 1$ Myr. Thus, at most two or three e-foldings of the initial rotational period can have occurred since the ZAMS. Assuming initially critical rotation, $\xi^1$ CMa should now have a rotational period of about 6 d; conversely, calculating from its actual rotation period, its initial rotation fraction $W_0$ must have already been very close to 0, implying a rotational period on the ZAMS already of the order of years.

Petit et al. (2013) calculated spin-down ages for all stars in Fig. 8 except $\beta$ CMa and $\epsilon$ CMa, finding $40 < t_{S,max} < 130$ Myr for these stars, i.e. their maximum spin-down ages are 4 to 10 times greater than their main-sequence ages, which for stars in this mass range are typically 8–15 Myr. We note that the rotational periods of $\tau$ Sco
and $\beta$ Cep are both known to high precision (Donati et al. 2006; Henrichs et al. 2013). For $\tau$ Sco, $t_{S,max} = 126_{-32}^{+16}$ Myr, while its position on the HRD indicates $t = 3.3_{-1.1}^{+0.5}$ Myr, a discrepancy of 1.6 dex; similarly for $\beta$ Cep, $t_{S,max} = 204_{-19}^{+9}$ Myr, differing by 1.2 dex with its age of $t = 13.2_{-0.7}^{+0.4}$ Myr (Shultz 2016).

These discrepancies between the stellar age and the maximum spin-down age cannot be reconciled by utilizing a different set of evolutionary models: while the absolute ages inferred from different prescriptions vary depending on assumptions about rotation, mixing or mass-loss, these uncertainties in the models are typically of the order of 1–2 Myr for a 15 $M_{\odot}$ star, much less than the ~1 dex difference between gyrochronological and evolutionary ages. It should be noted that the Ekström et al. (2012) models do not include magnetic fields. Fully self-consistent evolutionary models of magnetic stars are not yet available; however, it seems unlikely that such models will resolve the discrepancy via modification of the age inferred from a star’s position on the HRD.

Fossati et al. (2016) also found that spin-down time-scales are typically much longer than stellar ages for magnetic, massive stars, which they interpreted as evidence for rapid magnetic flux decay. At first order, the surface magnetic field should weaken as $1/R_*^2$ due to conservation of fossil magnetic flux. As an evolved star, $\xi^1$ CMa’s radius has increased by a factor of $1.8\pm0.1$ since the ZAMS; thus, its minimum surface magnetic field strength on the ZAMS should have been $B_{d,ZAMS} > 3.2$ kG. Recalculating $\tau_J$ on the ZAMS using Vink mass-loss rates suggests that $\tau_J$ should not have been very different from its present value, as the lower mass-loss rate on the ZAMS compensates for the stronger magnetic field and greater Alfvén radius. Thus, magnetic flux conservation alone cannot account for $\xi^1$ CMa’s slow rotation, and some other mechanism or combination of mechanisms – magnetic flux decay, rapid angular momentum loss on the pre-main sequence or modification of the internal structure of the star due to the magnetic field – is necessary. While it is of course possible that the theoretical mass-loss rates are simply incorrect, resolving the discrepancy by changing $\dot{M}$ alone would require an increase of ~1 dex, well outside the range of uncertainty of current models that typically differ by a factor of 2–3, and are in addition usually lower than the Vink et al. (2001) mass-loss rates (e.g. Sundqvist, Puls & Owocki 2014; Keszthelyi, Puls & Wade 2017).

The irreconcilability of gyrochronological and evolutionary time-scales amongst magnetic early B-type stars is worthy of future investigation. Here, it is pointed out primarily to emphasize that the difficulty in explaining $\xi^1$ CMa’s extremely slow rotation should not count as an argument against a long rotational period, as this problem is shared by all similar stars.

Fig. 23 shows $P_{\text{rot}}$ as a function of fractional main-sequence age $\tau_{\text{MS}}$ for those stars from Fig. 8 for which both $P_{\text{rot}}$ and $B_d$ are known. The majority of the dipolar field strengths and rotational periods were obtained from the catalogue published by Petit et al. (2013), and references therein. $\beta$ CMa’s magnetic, rotational and stellar properties were obtained from Fossati et al. (2015a). The rotational periods of NU Ori and HD 63425 were determined using ESPaDOnS data (Shultz 2016). Masses and fractional main-sequence ages were determined via linear interpolation between the Ekström et al. (2012) evolutionary models on the basis of the stars’ positions in Fig. 8. For HD 63425 and HD 66665, we used the $T_{\text{eff}}$–log g diagram to infer stellar parameters, as the uncertainties in log L are very large for these stars.

In addition to having by far the longest rotational period of any magnetic B-type star, $\xi^1$ CMa has a stronger magnetic field than any star in Fig. 23 but NU Ori. It is also amongst the most evolved

<figure><img src="image_26.png" /><figcaption>**Figure 23.** $P_{\text{rot}}$ as a function of fractional main-sequence age $\tau_{\text{MS}}$ for magnetic early B-type stars with masses similar to $\xi^1$ CMa, with colour corresponding to stellar mass and symbol size to $B_d$. Uncertainties in $P_{\text{rot}}$ are negligible on the scale of the diagram.</figcaption></figure>

stars in the sample: the only stars of a comparable fractional age are $\beta$ Cep and $\beta$ CMa, both of which are less massive (~12 $M_{\odot}$), and have much weaker magnetic fields (~0.1 kG). Since mass-loss rates increase with $M_*$, more massive stars should spin down more rapidly than less massive stars. With the exception of NU Ori, stars more massive than $\xi^1$ CMa all have rotational periods of 10–100 d, although they are less evolved. NU Ori, which has a stronger magnetic field, a much shorter rotational period and a slightly higher stellar mass than $\xi^1$ CMa, is also the youngest star on the diagram, likely accounting for its rapid rotation. This comparison shows that, viewed in the context of stars with similar magnetic and stellar parameters, $\xi^1$ CMa’s very slow rotation is not anomalous given its age.

## 9.3 A B-type star with an optically detectable DM?

In DMs, plasma deposited in the magnetosphere by the stellar wind falls back due to gravity, while in CMs centrifugal support due to the magnetically enforced corotation of the plasma with the star prevents infall above $R_K$. Petit et al. (2013) noted that stars without CMs only display emission if their mass-loss rates are high enough to replenish the magnetosphere on dynamical time-scales, thus enabling the circumstellar environment to remain optically thick despite continuous depletion as material returns to the photosphere. As a consequence, the only stars with optical emission lines originating in their DMs are the magnetic O-type stars, whereas the only magnetic B-type stars listed by Petit et al. as hosting $H\alpha$ emission possess CMs. While they were aware of $\xi^1$ CMa’s $H\alpha$ emission, Petit et al. assumed the 4.26 d rotation period found by Fourtune-Ravard et al. (2011), and thus that the emission originated in the star’s CM. However, as $R_K \gg R_A$, $\xi^1$ CMa does not have a CM.

The $H\alpha$ emission profile shown in Fig. 18 is furthermore entirely distinct from the double-horned profile characteristic of a CM, which generally produces two emission peaks located at approximately $\pm R_K$ (e.g. Leone et al. 2010; Grunhut et al. 2012a; Oksala et al. 2012; Rivinius et al. 2013). The strong central peak is much more similar to the typical profile of a magnetic O-type star with a DM (e.g. Sundqvist et al. 2012; ud-Doula et al. 2013; Wade et al. 2015), although the emission peaks of such stars tend
<figure><img src="image_27.png" /><figcaption>**Figure 24.** $R_A$ as a function of log L for magnetic early B-type stars (open black circles) and the magnetic O-type stars (filled blue squares). $\xi^1$ CMa is shown as a large, filled red circle. With the exception of $\xi^1$ CMa, only the O-type stars display emission in $H\alpha$. The ad hoc dashed line suggests a possible relationship between the onset of emission in DMs, and $R_A$ and log L.</figcaption></figure>

to be slightly redshifted, whereas $\xi^1$ CMa's is closer to line centre. $\xi^1$ CMa's ultraviolet emission line morphology is also consistent with an origin in a DM, at least in so far as it is very similar to the emission lines of $\beta$ Cep, which does not have a CM.

This leads to the question of how a B-type star can have an optically detectable DM. The quantity of plasma in a DM is a function of two parameters: the size of the magnetosphere and the mass-loss rate that feeds it. As is clear from Fig. 8, $\xi^1$ CMa is one of the most luminous magnetic B-type stars known: while $\tau$ Sco and NU Ori have similar log L, no stars are known to have a higher luminosity. In Fig. 24, $R_A$ is plotted as a function of log L for the stars in Fig. 8. Not only does $\xi^1$ CMa have a higher luminosity than most of the B-type stars, it also has one of the largest Alfvén radii.

Also plotted in Fig. 24 are the magnetic O-type stars from the Petit et al. (2013) catalogue. All of these stars have emission, despite the very small $R_A$ of the stars with higher luminosities. The dashed diagonal line indicates a possible division between stars with and without emission line DMs. Under the assumption that a star with a smaller $R_A$ requires a higher mass-loss rate in order to show emission, stars above the line should show emission, while those below should not. Alternatively, the presence or absence of emission could be related purely to the mass-loss rate, as indicated by the vertical dotted line: in this case, stars with $\log L/L_{\odot} \gtrsim 4.5$ should show emission. $\xi^1$ CMa's position in either case is consistent with the presence of $H\alpha$ emission. The presence or absence of magnetospheric emission in magnetic OB stars with $4.5 \le \log L/L_{\odot} \le 5.0$, but $R_A \ge 3R_*$, would help to distinguish between these scenarios.

Only one other magnetic B-type star, NU Ori, is in the same part of the diagram as $\xi^1$ CMa. As noted above, NU Ori is a rapid rotator, and possesses a CM. This star shows neither $H\alpha$ nor UV emission (Petit et al. 2013). The absence of emission, given its high luminosity and comparable magnetic confinement strength to $\xi^1$ CMa, is a distinct puzzle deserving of future investigation.

## 9.4 A pulsating magnetosphere?

Due to the limited number and quality of IUE observations, the evidence for pulsational modulation of UV resonance doublets is inconclusive. Evidence for pulsational modulation of the star's $H\alpha$ emission is much stronger. The residual $H\alpha$ EWs, after subtraction of a sinusoidal fit to the annual mean EWs phased with the presumed rotation period, yield a clear peak in the periodogram at the pulsation period. When phased with the pulsation period, the amplitude of the residual EW variation is three times higher than the amplitude expected from a photospheric modulation of the absorption line's EW, and is furthermore 0.5 pulsation cycles out of phase with the predicted photospheric variation. Finally, $H\alpha$ shows a much stronger $V/R$ variation than is predicted by the photospheric model.

Stochastic changes in the EW of the $H\alpha$ line have been observed for some magnetic O-type stars with DMs, e.g. the well-studied object $\theta^1$ Ori C (Stahl et al. 2008). 3D MHD simulations are able to reproduce this short-term variability as a consequence of turbulent plasma infall (ud-Doula et al. 2013). As $\dot{M}$ is constant for such stars, their short-term variability is wholly a consequence of the complex behaviour of MHD plasmas. Since $\xi^1$ CMa's mass-loss rate should change with the pulsation phase, it is conceivable that this could force a periodicity in the outflow and infall of the magnetically confined plasma. The observation of prominent blue- and redshifted secondary peaks in the $H\alpha$ emission structure occurring at distinct pulsation phases, together with the stronger modulation in $H\alpha V/R$ as compared to residual $H\alpha$ EWs, further suggests that, throughout the course of a pulsation cycle, there are larger local changes in the distribution of plasma within the magnetosphere, as compared to the total mass of magnetically confined plasma.

The majority of the $H\alpha$ emission is expected to be produced in the innermost regions (ud-Doula et al. 2013). The wind-flow time-scale, calculated using a standard $\beta = 1$ velocity law,<fnref n="7" /> ranges from ~0.13 d at a distance of $2 R_*$ to ~0.19 d at the minimum Alfvén radius of $4.6 R_*$. The free-fall time-scale is somewhat longer, ranging from 0.2 to 1.2 d over the same distance. These time-scales are similar to the pulsation period, suggesting that periodic density waves arising due to the varying mass-loss rate might lead to a periodicity in the plasma infall. The strongest blueshifted secondary emission peak occurs at about phase 1.0, and the strongest secondary redshifted emission peak at phase 0.5. If the blue- and redshifted emission peaks are respectively consequences of outflowing and infalling plasma, and assuming that the majority of the emission is formed relatively close to the star, a lag of approximately 0.7 cycles between the peak mass-loss rate (at phase 0.3) and the strongest blueshifted emission (at phase 1.0) would make sense.

While no other magnetic, massive star's $H\alpha$ emission has ever been found to be modulated by pulsation as well as rotation, Oskinova et al. (2014) reported that $\xi^1$ CMa's X-ray light curve exhibits a coherent variation with the pulsation period, with both hard and soft X-ray emission peaking ~0.1 cycles after the maximum of the visible light curve.

Since $\xi^1$ CMa is a pulsating star, $\dot{M}$ is not constant. The changes in $T_{\text{eff}}$ and $\log L_*$ should lead to a variation in $\dot{M}$ with a semi-amplitude of 0.05 dex. Changes in $v_{\infty}$ should be negligible. Calculation of $\eta_*$ and $R_A$ as a function of pulsation phase suggests the overall changes in these parameters to be small (~20 per cent in $\eta_*$ and

<footnote n="7">⁷ While confinement in a dipolar magnetic field does modify the wind-flow time-scale, in practice this is a second-order effect (ud-Doula et al. 2014).</footnote>
4 per cent in $R_A$). Similar X-ray pulsations have also been reported for the $\beta$ Cep star $\beta$ Cru (Cohen et al. 2008). As $\beta$ Cru does not have a detected magnetic field (Grunhut, private communication), this argues against a purely magnetic explanation for the X-ray variability.

Oskinova et al. (2014) considered a variation in $\dot{M}$ as a possible mechanism behind the star’s X-ray pulsations, and rejected it based upon the absence of X-ray pulsations in other $\beta$ Cep stars, either magnetic or not. However, they did not consider their results in the context of the XADM model. The XADM model predicts a modulation in $\log L_X$ of approximately $\pm 10$ per cent to arise from a variation in $\dot{M}$ of $\pm 0.05$ dex, with the maximum and minimum of the predicted X-ray light curve corresponding to the maximum and minimum of the visible light curve. This is the same amplitude as that reported by Oskinova et al. (2014).

As a sanity check, we performed similar calculations for $\beta$ Cep and $\beta$ Cen, the only other magnetic $\beta$ Cep stars with X-ray light curves. $\beta$ Cep is slightly cooler ($T_{\text{eff}} \sim 26$ K), less strongly magnetized ($B_d \sim 260$ G; Henrichs et al. 2013) and has somewhat lower amplitude pulsations leading to a smaller change in luminosity ($\pm 0.015$ mag and $\pm 0.006$ dex, respectively, as evaluated by phasing the star’s *Hipparcos* observations with its primary pulsation period of 0.19 d; Shibahashi & Aerts 2000) and, hence, a smaller change in the mass-loss rate ($\pm 0.03$ dex). The XADM model predicts that $\beta$ Cep should have X-ray pulsations with a semi-amplitude of $\le 6$ per cent, close to the upper limit determined by Favata et al. (2009). For $\beta$ Cen, it is not clear which of the many pulsation frequencies identified by Pigulski et al. (2016) belongs to the magnetic secondary. However, even the highest amplitude frequency ($\pm 0.0015$ mag) leads to a luminosity modulation of just 0.006 dex, hence a change in $\dot{M}$ of less than $\pm 0.01$ dex, and, with $B_d = 200$ G (Shultz 2016), a variation in the X-ray luminosity of $< 3$ per cent, below the upper limit established by Raassen et al. (2005).

It is probably not the case that simple calculation of global parameters at each pulsation phase provides an accurate picture. Strong local density variations due to the changing $\dot{M}$ may lead to correspondingly strong local changes in the magnetic confinement strength. MHD simulations will be necessary to properly explore the impact of stellar pulsation on magnetospheric dynamics.

# 10 CONCLUSIONS AND FUTURE WORK

The principal result of this paper is that both magnetic and spectroscopic data indicate that the rotation period of $\xi^1$ CMa is very long, of the order of decades. Rotational periods of the order of days, as previously suggested in the literature, are conclusively ruled out by the ESPaDOnS $\langle B_z \rangle$ measurements alone, which require a period of at least 5100 d. Inclusion of the MuSiCoS results requires that $P_{\text{rot}}$ be about 30 yr, a conclusion supported by the modulation of the $H\alpha$ emission strength seen in the CORALIE and ESPaDOnS data sets. A period of approximately 30 yr is consistent with the timing of the star’s strong UV emission, which should correspond to an extremum of the $\langle B_z \rangle$ curve. Given that phase coverage of the rotational cycle is incomplete to an unknown degree, the constraints that can be placed upon the geometry and strength of the magnetic field are necessarily somewhat loose. Further spectropolarimetric observation over decades is essential to characterization of the magnetic field. If the 30 yr period is correct, magnetic crossover ($\langle B_z \rangle = 0$) will occur in 2018, and magnetic minimum in 2025.

We have conducted the first detailed examination of the star’s weak $H\alpha$ emission, which we detect in all available spectroscopic data. There is no evidence in the interferometric data of a binary

companion, and combined constraints from interferometry and RV residuals rule out a binary companion of sufficient luminosity to host the weak emission. Furthermore, there is no indication of a Kepplerian disc in the high-resolution AMBER spectro-interferometry. Thus, the interferometric data do not support the hypothesis that the emission originates in the Keplerian disc of a classical Be companion star, indicating that this emission is produced in $\xi^1$ CMa’s circumstellar environment. This makes $\xi^1$ CMa the first magnetic $\beta$ Cep star with a magnetosphere detectable in visible light. Given that its extremely slow rotation implies a DM rather than a CM, $\xi^1$ CMa is also the coolest star with a DM detectable in $H\alpha$.

Finally, in addition to the already reported modulation of X-ray flux with the pulsation period, the $H\alpha$ emission also appears to correlate with the pulsations in a fashion that seems likely to be intrinsic to the circumstellar plasma. The best evidence for this is in the strong $H\alpha V/R$ variations, which are much greater than expected from a radially pulsating photospheric model. The total $H\alpha$ EW in any given epoch also varies in antiphase with the expected EW variation due to photospheric $T_{\text{eff}}$ variation, and furthermore has a larger amplitude than predicted by a photospheric model. However, the EW is only weakly variable compared to $V/R$, suggesting that the distribution of plasma within the magnetosphere, rather than the total quantity of magnetically confined material, is the primary origin of the variability. Given the complex velocity fields within DMs, MHD simulations will be necessary to understand the interplay between pulsation, mass-loss and magnetic confinement. As the magnetic wind confinement parameter $\eta_*$ is of the order of $10^2$, modelling $\xi^1$ CMa’s magnetosphere should be within the capabilities of the current generation of 2D MHD codes.

We also find a weak modulation of $\langle B_z \rangle$ with the pulsation period. Further modelling is needed to understand the origin of this phenomenon. Signatures in the diagnostic null profile are accurately reproduced as a consequence of RV variations between sub-exposures; however, this model does not predict any modulation of $\langle B_z \rangle$. Polarized radiative transfer in a moving atmosphere will be essential to answering the question of whether these modulations are an artefact of the measurement methodology, or intrinsic to the stellar magnetic field.

# ACKNOWLEDGEMENTS

This work is based on observations obtained using the MuSiCoS spectropolarimeter at Pic du Midi Observatory; ESPaDOnS at the Canada-France-Hawaii Telescope (CFHT), which is operated by the National Research Council of Canada, the Institut National des Sciences de l'Université de la Centre Nationale de la Recherche Scientifique of France and the University of Hawaii; on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programme 292.D-5028(A) at the Paranal Observatory, ESO Chile, with the NIR interferometers AMBER and PIONIER; on observations obtained at the La Silla Observatory, ESO Chile, with the echelle spectrograph CORALIE at the 1.2 m Euler Swiss telescope; and on NEWSIPS data from the *IUE* satellite. We acknowledge the assistance of Conny Aerts, who provided the CORALIE data set. This work has made use of the VALD data base, operated at UPSala University, the Institute of Astronomy RAS in Moscow and the University of Vienna. MS acknowledges the financial support provided by the European Southern Observatory studentship programme in Santiago, Chile. MS and GAW acknowledge support from the Natural Science and Engineering Research Council of Canada (NSERC). CN and the MiMeS Collaboration acknowledge
financial support from the Programme National de Physique Stel-
laire (PNPS) of INSU/CNRS. We acknowledge the Canadian As-
tronomy Data Centre (CADC). The *IUE* data presented in this paper
were obtained from the Multimission Archive at the Space Tele-
scope Science Institute (MAST). STScI is operated by the Associ-
ation of Universities for Research in Astronomy, Inc., under NASA
contract NAS5-26555. Support for MAST for non-*HST* data is pro-
vided by the NASA Office of Space Science via grant NAG5-7584
and by other grants and contracts. This project used the facilities of
SIMBAD and *Hipparcos*. This research has made use of the Jean-
Marie Mariotti *LITPRO* service co-developed by CRAL, IPAG and
LAGRANGE.

REFERENCES

Aerts C., De Cat P., 2003, Space Sci. Rev., 105, 453

Aerts C., Christensen-Dalsgaard J., Kurtz D. W., 2010, Asteroseismology.
Springer, Berlin

Aerts C., Simón-Díaz S., Groot P. J., Degroote P., 2014, A&A, 569, A118

Aleilian E. et al., 2011, A&A, 536, L6

Aleilian E. et al., 2015, in Meynet G., Georgy C., Groh J., Stee P., eds, Proc. IAU Symp. 307, New Windows on Massive Stars: Asteroseismology, Interferometry, and Spectropolarimetry. Cambridge Univ. Press, Cambridge, p. 330

Bagnulo S., Landstreet J. D., Fossati L., Kochukhov O., 2012, A&A, 538, A129

Bagnulo S., Fossati L., Landstreet J. D., Izzo C., 2015, A&A, 583, A115

Baranne A. et al., 1996, A&AS, 119, 373

Baudrand J., Bohm T., 1992, A&A, 259, 711

Bessell M. S., Castelli F., Plez B., 1998, A&A, 333, 231

Bloomfield P., 1976, Wiley Series in Probability and Mathematical Statistics.
Wiley, New York

Brandt J. C. et al., 1998, AJ, 116, 941

Breger M. et al., 1993, A&A, 271, 482

Cassinelli J. P., Olson G. L., 1979, ApJ, 229, 304

Cassinelli J. P., Cohen D. H., Macfarlane J. J., Sanders W. T., Welsh B. Y., 1994, ApJ, 421, 705

Chini R., Hoffmeister V. H., Nasseri A., Stahl O., Zinnecker H., 2012, MNRAS, 424, 1925

Claret A., 2004, A&A, 424, 919

Clayton G. C., Wolff M. J., Sofia U. J., Gordon K. D., Misselt K. A., 2003, ApJ, 588, 871

Cohen D. H., Cassinelli J. P., MacFarlane J. J., 1997, ApJ, 487, 867

Cohen D. H., Kuhn M. A., Gagné M., Jensen E. L. N., Miller N. A., 2008, MNRAS, 386, 1855

Davey W. R., 1973, ApJ, 179, 235

Díaz-Cordovés J., Claret A., Giménez A., 1995, A&AS, 110, 329

Donati J.-F., Semel M., Rees D. E., 1992, A&A, 265, 669

Donati J.-F., Semel M., Carter B. D., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658

Donati J.-F., Catala C., Wade G. A., Gallou G., Delaigue G., Rabou P., 1999, A&AS, 134, 149

Donati J.-F., Wade G. A., Babel J., Henrichs H. F., de Jong J. A., Harries T. J., 2001, MNRAS, 326, 1265

Donati J.-F. et al., 2006, MNRAS, 370, 629

Ducati J. R., 2002, VizieR Online Data Catalog, 2237, 0

Ekström S. et al., 2012, A&A, 537, A146

Favata F., Neiner C., Testa P., Hussain G., Sanz-Forcada J., 2009, A&A, 495, 217

Fitzpatrick E. L., Massa D., 1990, ApJS, 72, 163

Fossati L. et al., 2015a, A&A, 574, A20

Fossati L. et al., 2015b, A&A, 582, A45

Fossati L. et al., 2016, A&A, 592, A84

Fourtune-Ravard C., Wade G. A., Marcolino W. L. F., Shultz M., Grunhut J. H., Henrichs H. F., Henrichs, 2011, in Neiner C., Wade G., Meynet G., Peters G., eds, Proc. IAU Symp. 272, Active OB Stars: Structure, Evolu-

tion, Mass Loss, and Critical Limits. Cambridge Univ. Press, Cambridge,
p. 180

Frost E. B., 1907, ApJ, 25, 59

Fullerton A. W., Gies D. R., Bolton C. T., 1996, ApJS, 103, 475

Gray R. O., 1998, AJ, 116, 482

Grunhut J. H., 2012, PhD thesis, Queen's University, Canada

Grunhut J. H. et al., 2012a, MNRAS, 419, 1610

Grunhut J. H. et al., 2012b, MNRAS, 426, 2208

Hauck B., Mermilliod M., 1998, A&AS, 129, 431

Henrichs H. F., 1993, in Nemec J. M., Matthews J. M., eds, IAU Colloq. 139: New Perspectives on Stellar Pulsation and Pulsating Variable Stars. Cambridge Univ. Press, Cambridge, p. 186

Henrichs H. F. et al., 2013, A&A, 555, A46

Heynderickx D., Waelkens C., Smeyers P., 1994, A&AS, 105, 447

Hubrig S., Briquet M., Schöller M., De Cat P., Mathys G., Aerts C., 2006, MNRAS, 369, L61

Hubrig S., Ilyin I., Schöller M., Briquet M., Morel T., De Cat P., 2011, ApJ, 726, L5

Jerzykiewicz M., 1999, New Astron. Rev., 43, 455

Keszthelyi Z., Puls J., Wade G. A., 2017, A&A, 598, A4

Kochukhov O., Bagnulo S., Wade G. A., Sangalli L., Piskunov N., Landstreet J. D., Petit P., Sigut T. A. A., 2004, A&A, 414, 613

Kochukhov O., Makaganiuk V., Piskunov N., 2010, A&A, 524, A5

Krtička J., 2014, A&A, 564, A70

Kupka F. G., Piskunov N., Ryabchikova T. A., Stempels H. C., Weiss W. W., 1999, A&AS, 138, 119

Kupka F. G., Ryabchikova T. A., Piskunov N. E., Stempels H. C., Weiss W. W., 2000, Balt. Astron., 9, 590

Kuschnig R., Weiss W. W., Gruber R., Bely P. Y., Jenkner H., 1997, A&A, 328, 544

Lachaume R., 2003, A&A, 400, 795

Landstreet J. D., Mathys G., 2000, A&A, 359, 213

Lanz T., Hubeny I., 2007, ApJS, 169, 83

Le Bouquin J.-B. et al., 2011, A&A, 535, A67

Lefever K., Puls J., Morel T., Aerts C., Decin L., Briquet M., 2010, A&A, 515, A74

Lenz P., Breger M., 2005, Commun. Asteroseismol., 146, 53

Leone F., Bohlender D. A., Bolton C. T., Buemi C., Catanzaro G., Hill G. M., Stift M. J., 2010, MNRAS, 401, 2739

Lucy L. B., 2007, A&A, 468, 649

Lucy L. B., 2010, A&A, 512, A33

McNamara D. H., 1955, ApJ, 122, 95

Marcolino W. L. F. et al., 2012, MNRAS, 422, 2314

Marcolino W. L. F., Bouret J.-C., Sundqvist J. O., Walborn N. R., Fullerton A. W., Howarth I. D., Wade G. A., ud-Doula A., 2013, MNRAS, 431, 2253

Mason B. D., Wycoff G. L., Hartkopf W. I., Douglass G. G., Worley C. E., 2001, AJ, 122, 3466

Mathys G., 1989, Fundam. Cosm. Phys., 13, 143

Moon T., 1985, Ap&SS, 117, 261

Morel T., Hubrig S., Briquet M., 2008, A&A, 481, 453

Muijres L. E., Vink J. S., de Koter A., Müller P. E., Langer N., 2012, A&A, 537, A37

Nardetto N., Mathias P., Fokin A., Chapellier E., Pietrzynski G., Gieren W., Graczyk D., Mourard D., 2013, A&A, 553, A112

Nazé Y., Petit V., Rinbrand M., Cohen D., Owocki S., ud-Doula A., Wade G. A., 2014, ApJS, 215, 10

Neilson H. R., Ignace R., 2015, A&A, 584, A58

MiMeS CollaborationNeiner C., Landstreet J. D., Alecian E., Owocki S., Kochukhov O., Bohlender D., 2012, A&A, 546, A44

Niemczura E., Daszyńska-Daszkiewicz J., 2005, A&A, 433, 659

Oksala M. E., Wade G. A., Townsend R. H. D., Owocki S. P., Kochukhov O., Neiner C., Alecian E., Grunhut J., 2012, MNRAS, 419, 959

Oskinova L. M., Todt H., Ignace R., Brown J. C., Cassinelli J. P., Hamann W.-R., 2011, MNRAS, 416, 1456

Oskinova L. M., Nazé Y., Todt H., Huenemoerder D. P., Ignace R., Hubrig S., Hamann W.-R., 2014, Nat. Commun., 5, 4024
Oskinova L. M., Todt H., Huenemoerder D. P., Hubrig S., Ignace R., Hamann W.-R., Balona L., 2015, A&A, 577, A32

Petit V., Wade G. A., 2012, MNRAS, 420, 773

Petit V. et al., 2013, MNRAS, 429, 398

Petrov R. G. et al., 2007, A&A, 464, 1

Pigulski A. et al., 2016, A&A, 588, A55

Piskunov N. E., Kupka F., Ryabchikova T. A., Weiss W. W., Jeffery C. S., 1995, A&AS, 112, 525

Preston G. W., 1967, ApJ, 150, 547

Queloz D. et al., 2000, A&A, 354, 99

Queloz D. et al., 2001, The Messenger, 105, 1

Raassen A. J. J., Cassinelli J. P., Miller N. A., Mewe R., Tepedelenlioğlu E., 2005, A&A, 437, 599

Rivinius T., Townsend R. H. D., Kochukhov O., Štefli S., Baade D., Barrera L., Szeifert T., 2013, MNRAS, 429, 177

Ryabchikova T. A., Piskunov N. E., Kupka F., Weiss W. W., 1997, Balt. Astron., 6, 244

Saesen S., Briquet M., Aerts C., 2006, Commun. Asteroseismol., 147, 109

Schnerr R. S., Henrichs H. F., Oudmaijer R. D., Telting J. H., 2006, A&A, 459, L21

Schnerr R. S. et al., 2008, A&A, 483, 857

Shibahashi H., Aerts C., 2000, ApJ, 531, L143

Shorlin S. L. S., Wade G. A., Donati J.-F., Landstreet J. D., Petit P., Sigut T. A. A., Strasser S., 2002, A&A, 392, 637

Shultz M. E., 2016, PhD thesis, Department of Physics, Engineering Physics, and Astrophysics, Queen's University, Canada

Shultz M., Wade G., Rivinius T., Marcolino W., Henrichs H., Grunhut J., 2015, in Meynet G., Georgy C., Groh J., Stee P., eds, Proc. IAU Symp. 307, New Windows on Massive Stars: Asteroseismology, Interferometry, and Spectropolarimetry. Cambridge Univ. Press, Cambridge, p. 399

Silvester J. et al., 2009, MNRAS, 398, 1505

Smith M. A., Groote D., 2001, A&A, 372, 208

Stahl O., Wade G., Petit V., Stober B., Schanne L., 2008, A&A, 487, 323

Štefli S. et al., 2009, A&A, 504, 929

Stothers R., 1976, ApJ, 210, 434

Sundqvist J. O., ud-Doula A., Owocki S. P., Townsend R. H. D., Howarth I. D., Wade G. A., 2012, MNRAS, 423, L21

Sundqvist J. O., Simón-Díaz S., Puls J., Markova N., 2013, A&A, 559, L10

Sundqvist J. O., Puls J., Owocki S. P., 2014, A&A, 568, A59

Tallon-Bosc I. et al., 2008, in Schöller M., Danchi W. C., Delplancke F., eds, Proc. SPIE, Vol. 7013, Optical and Infrared Interferometry. SPIE, Bellingham, p. 70131J

Townsend R. H. D., Owocki S. P., 2005, MNRAS, 357, 251

ud-Doula A., Owocki S. P., 2002, ApJ, 576, 413

ud-Doula A., Owocki S. P., Townsend R. H. D., 2008, MNRAS, 385, 97

ud-Doula A., Owocki S. P., Townsend R. H. D., 2009, MNRAS, 392, 1022

ud-Doula A., Sundqvist J. O., Owocki S. P., Petit V., Townsend R. H. D., 2013, MNRAS, 428, 2723

ud-Doula A., Owocki S., Townsend R., Petit V., Cohen D., 2014, MNRAS, 441, 3600

van Leeuwen F., 2007, A&A, 474, 653

Vink J. S., de Koter A., Lamers H. J. G. L. M., 2001, A&A, 369, 574

Wade G. A. et al., 2006, A&A, 458, 569

Wade G. A. et al., 2015, MNRAS, 447, 2551

Wade G. A. et al., 2016, MNRAS, 456, 2

Walborn N. R., Parker J. W., Nichols J. S., 1995, VizieR Online Data Catalog, 3188, 0

Weber E. J., Davis L., Jr, 1967, ApJ, 148, 217

Worthey G., Lee H.-c., 2011, ApJS, 193, 1

<caption>**Table C2.** Log of CORALIE observations, RV measurements and H$\alpha$ EW measurements. Peak signal-to-noise ratios (S/N) are per spectral pixel.</caption>

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the article.

## APPENDIX A: RELIABILITY OF MUSICOS DATA

Fossati et al. (2015b) objected to a long rotation period, first suggested by Shultz et al. (2015), on the grounds that the period is largely constrained by the MuSiCoS measurements, and further suggesting that as these data were collected at high airmass the measurements were unreliable. As we have shown, it is not the case that a long rotation period is implied primarily by the MuSiCoS data: to the contrary, a gradual decline in $\langle B_z \rangle$ is clearly evident in the ESPaDOnS data alone, and the conclusion that the star is slowly rotating is also supported by the long-term modulation of the H$\alpha$ emission strength. However, the MuSiCoS measurements are an important constraint, and in this appendix we address the question of their reliability. While we are aware of no reason why the sign of the Stokes V signature should be flipped due to the star being observed at high airmass, the fibres were plugged in differently each observing run, which could result in an apparent polarity change. To check this possibility, we acquired the MuSiCoS data set of the Bp star 36 Lyn reported by Wade et al. (2006), of which one measurement was acquired on 02/10/2000, during the same observing run as the $\xi^1$ CMa data. We compare these data with Narval observations of 36 Lyn. Narval, a clone of ESPaDOnS, obtains essentially identical results to its sibling instrument (Wade et al. 2016). LSD profiles were extracted from the two data sets using the same method described in Section 6.1, with a line mask appropriate to 36 Lyn's $T_{eff}$.

Fig. A1 shows the resulting $\langle B_z \rangle$ measurements, phased with the rotational ephemeris determined by Wade et al. (2006). The two data sets agree well at all phases. While the measurement obtained on 02/10/2000 is in its expected position, $\langle B_z \rangle = -40 \pm 170$ G is close to 0; thus, a sign flip would not change its position on the $\langle B_z \rangle$ curve. As an additional check, Fig. A2 compares the LSD profile of this MuSiCoS observation to the Narval LSD profile with the closest rotational phase. The Stokes V profiles of the two observations are a good match to one another. Moreover, a change in polarity of the MuSiCoS signature would render it in conflict with the Narval

<figure><img src="image_28.png" /><figcaption>**Figure A1.** $\langle B_z \rangle$ measurements of 36 Lyn obtained with Narval and MuSiCoS, phased with the ephemeris determined from the MuSiCoS data by Wade et al. (2006). The MuSiCoS measurement obtained during the same observing run as the $\xi^1$ CMa observations is indicated by the filled red circle.</figcaption></figure>

## SUPPORTING INFORMATION

Supplementary data are available at MNRAS online.

**Table B1.** Log of spectropolarimetric observations.

**Table C1.** RVs and H$\alpha$ EWs from ESPaDOnS sub-exposures.
<figure><img src="image_29.png" /><figcaption>**Figure A2.** Comparison between the MuSiCoS LSD profile of the magnetic Bp star 36 Lyn obtained during the same observing run as the $\xi^1$ CMa data and a Narval LSD profile of 36 Lyn at a similar rotational phase.</figcaption></figure>

profiles. We conclude that there is no compelling reason to doubt the reliability of MuSiCoS, which obtained results fully consistent with those of modern high-resolution spectropolarimeters, and thus that the negative polarity of $\xi^1$ CMa's Stokes V profile during this epoch is real.

APPENDIX B: LOG OF MAGNETIC MEASUREMENTS

APPENDIX C: RV AND H$\alpha$ EW MEASUREMENTS

APPENDIX D: LOG OF INTERFEROMETRIC OBSERVATIONS

<table><caption>**Table D1.** Log of VLTI observations. Type is photometry (H or K band) or high-resolution (HR).</caption><thead><tr><th>Cal. date</th><th>Type</th><th>VLTI config.</th></tr></thead><tbody><tr><td></td><td>AMBER</td><td></td></tr><tr><td>2014-03-08</td><td>H</td><td>A1-K0-G1</td></tr><tr><td>2014-03-08</td><td>K</td><td>A1-K0-G1</td></tr><tr><td>2014-03-08</td><td>HR</td><td>A1-K0-G1</td></tr><tr><td></td><td>PIONIER</td><td></td></tr><tr><td>2014-03-23</td><td>H</td><td>A1-G1-J3-K0</td></tr></tbody></table>

This paper has been typeset from a TeX/Latex file prepared by the author.