The spectral variability and magnetic field characteristics of the Of?p star HD 148937

We report magnetic and spectroscopic observations and modeling of the Of?p star HD 148937 within the context of the MiMeS LP at the CFHT. Thirty-two high signal-to-noise ratio circularly polarised (Stokes V) spectra and 13 unpolarised (Stokes I) spectra of HD 148937 were acquired in 2009 and 2010. A definite detection of a Stokes V Zeeman signature is obtained in the grand mean of all observations (in both LSD mean profiles and individual spectral lines). The longitudinal magnetic field inferred from the Stokes V LSD profiles is consistently negative, in contrast to the essentially zero field strength measured from the diagnostic null profiles. A period search of equivalent width measurements confirms the previously-reported 7.03 d variability period. The variation of equivalent widths is not strictly periodic: we present evidence for evolution of the amount or distribution of circumstellar plasma. Interpreting the 7.03 d period as the stellar rotational period within the context of the ORM, we have phased the equivalent widths and longitudinal field measurements. The longitudinal field measurements show a weak sinusoidal variation of constant sign, with extrema out of phase with the H{\alpha} variation by about 0.25 cycles. The inferred magnetic configuration confirms the suggestion of Naz\'e et al (2010), who proposed that the weaker variability of HD 148937 as compared to other members of this class is a consequence of the stellar geometry. Based on the derived magnetic properties and published wind characteristics, we find a wind magnetic confinement parameter \eta\ast \simeq 20 and rotation parameter W = 0.12, supporting a picture in which the Halpha emission and other line variability have their origin in an oblique, rigidly rotating magnetospheric structure resulting from a magnetically channeled wind. (Abridged.)


INTRODUCTION
The enigmatic Of?p stars are identified by a number of peculiar and outstanding observational properties. The classification was first introduced by Walborn (1972) according to the presence of C iii λ4650 emission with a strength comparable to the neighbouring N iii lines. Well-studied Of?p stars are now known to exhibit recurrent, and apparently periodic, spectral variations (in Balmer, He i, C iii and Si iii lines), narrow P Cygni or emission components in the Balmer lines and He i lines, and UV wind lines weaker than those of typical Of supergiants (see Nazé et al. 2010 and references therein).
Only 5 Galactic Of?p stars are known (Walborn et al. 2010): HD 108, HD 148937, HD 191612, NGC 1624-2 and CPD−28 • 2561. Three of these stars -HD 108, HD 148937 and HD 191612 -have been studied in detail. In recent years, HD 191612 and HD 108 have been carefully examined for the presence of magnetic fields (Donati et al. 2006;Martins et al. 2010;Wade et al 2011), and both have been clearly detected. In retrospect, the similarity of some of the observational properties of the Of?p stars to the O7 V magnetic oblique rotator θ 1 Ori C appears to have been a strong indicator that they are also oblique rotators with strong, stable, organised magnetic fields. As the brightest known Of?p star, and in light of the report of the marginal detection of a longitudinal magnetic field by Hubrig et al. (2008Hubrig et al. ( , 2011, the time is ripe for a more intensive investigation of HD 148937. HD 148937 has been extensively observed; its observational properties are discussed in the literature by Nazé et al. (2008aNazé et al. ( , 2010. Similar to other well-studied Of?p stars, HD 148937 is revealed to be a high mass (50-60 M ⊙ ) main sequence O-type star with an effective temperature of about 40 kK. The physical parameters of this star (as first derived by Nazé et al. 2008a and revised in this paper) are reported in Table 1. The HR diagram positions of the 3 well-studied Of?p stars are illustrated in Fig. 1. While it is clear that the Of?p stars are high-mass O-type stars, their HR diagram positions are relatively uncertain due to their poorly constrained distances. This translates into important absolute and relative uncertainties on their masses and ages, as is reflected in the relative large error bars in Fig. 1.
Like HD 108 and HD 191612, HD 148937 is a spectroscopic variable star. However, it is distinguished from those stars by its lower-amplitude variability, particularly in the primary diagnostic Hα line. Whereas HD 108 and HD 191612 exhibit Hα equivalent width (EW) variability of more than a factor of 2, the Hα line of HD 148937 changes in EW by only about 20% (e.g. Nazé et al. 2008a). Despite this much weaker variability, Nazé et al. (2008a) were able to tentatively identify a period of 7.031 ± 0.003 d in EW measurements of HD 148937 spanning over 3 years. This period was subsequently confirmed using higher quality data by Nazé et al. (2010). Photometry of HD 148937 reveals no significant broadband flux variation (Balona 1992, van Genderen et al. 1989, 2001, Nazé et al. 2008a. Nazé et al. (2008a) also found that this star exhibits an enhanced nitrogen abundance (by about a factor of four relative to the sun). At X-ray wavelengths, HD 148937 resembles HD 108 and HD 191612 in having a thermal spectrum dominated by a relatively cool component (kT = 0.2 keV), broad lines (> 1700 km s −1 ; although see detection of narrower features in high-Z lines by Nazé et al. 2011), and an order-of-magnitude overluminosity compared to normal O stars (log[L unabs X /L BOL ] ∼ −6). In this paper we perform a first detailed investigation of the combined magnetic and variability properties of HD 148937 using According to its position on this diagram, HD 148937 is the youngest, and most massive and luminous of the well-studied Of?p stars (although it should be noted that the detailed positions are subject to significant uncertainties in the distances to these stars, as reflected in the error bars). Evolution tracks and isochrones include rotational mixing, and are from Meynet & Maeder (2005). Adapted from Martins et al. (2010). an extensive high-resolution spectropolarimetric and spectroscopic dataset. In Sect. 2 we discuss the data acquired and the methods of analysis used. In Sect. 3 we re-examine the physical properties of the star, as well as its projected rotational velocity. In Sect. 4 we examine the spectral variability, confirming the periodic variability of the Hα and other emission and absorption lines, and the period derived by Nazé et al. (2008aNazé et al. ( , 2010. In Sect. 5 we analyse in detail the magnetic data acquired at the Canada-France-Hawaii Telescope (CFHT). In Sect. 6 we employ the CFHT magnetic data, in combination with the stellar physical and rotational parameters, to constrain the magnetic field strength and geometry. Finally, in Sect. 7 we explore the implications of our study, particularly regarding the variability and other properties of HD 148937, the confinement and structure of its stellar wind, and of the properties of the general class of Of?p stars.

ESPaDOnS spectropolarimetry
Spectropolarimetric observations of HD 148937 were obtained using the ESPaDOnS spectropolarimeter at the CFHT in 2009 and 2010 within the context of the Magnetism in Massive Stars (MiMeS) Large Program and Project. Altogether, 32 Stokes V sequences were obtained -two in May 2009, two in September 2009, 14 during 7 consecutive nights in June 2010, and 14 during 7 consecutive nights in July 2010. In particular, the final observing runs were planned so as to sample two complete rotations of the star, hypothesising that the 7.03 d period identified by Nazé et al. was in fact the stellar rotational period.
Each polarimetric sequence consisted of four individual Table 1. Summary of stellar physical, wind and magnetic properties of HD 148937, derived by Nazé et al. (2008a) and in this paper. The wind magnetic confinement parameter η * , the rotation parameter W and the characteristic spin down time τ spin are defined and described in Sect. 7.

Spectral type
This paper log(L X /L Bol ) -6.1 Nazé et al. (2008a) subexposures taken in different polarimeter configurations. From each set of four subexposures we derive a mean Stokes V spectrum following the procedure of Donati et al. (1997), ensuring in particular that all spurious signatures are removed at first order. Diagnostic null polarization spectra (labeled N) are calculated by combining the four subexposures in such a way that polarization cancels out, allowing us to check that no spurious signals are present in the data (see Donati et al. 1997 for more details on how N is defined). All frames were processed using the Upena pipeline feeding Libre ES-pRIT (Donati et al. 1997), a fully automatic reduction package installed at CFHT for optimal extraction of ESPaDOnS spectra. The peak signal-to-noise ratios per 1.8 km s −1 velocity bin in the reduced spectra range from 60 to nearly 1000, with a median of 690, depending on the exposure time and on weather conditions. Due to the low declination of HD 148937, the star was generally observed at high airmass (2.8 − 3.5) from CFHT. The log of CFHT observations is presented in Table 2.

FEROS spectroscopy
In addition to the 32 circular polarisation spectra, 13 highresolution unpolarised spectra were acquired on 10 nights in 2009 March using the FEROS spectrograph on the MPG/ESO 2.2-m telescope at La Silla. These spectra cover a useful wavelength range of ∼360-920 nm at a resolving power of R ≃ 48, 000, with a S/N of ∼125 at Hα. Theéchelle data were extracted and merged using the pipeline FEROS Data Reduction System (see Kaufer et al. 1999 for details). The log is given in Table 3.

Coralie spectroscopy
Finally, we have re-analysed the 20échelle spectra of HD 148937 obtained in 2008 with the Coralie spectrograph on the 1.2m Euler Swiss telescope at La Silla (Chile), described by Nazé et al. (2010).

PHYSICAL AND ROTATIONAL PROPERTIES
We have re-determined the main stellar properties of HD 148937 using the atmosphere code CMFGEN (Hillier & Miller 1998). The results are basically unchanged compared to the analysis of Nazé et al. (2008a). The only minor difference is a slightly higher luminosity (log L/L ⊙ = 5.8 vs. 5.75). We consider this to be a better determination, based on the fit of the UV-optical-IR SED rather than just on the V-band magnitude. We also derive the colour excess E(B − V) = 0.67 from this fitting process, assuming a classical Galactic extension law (Seaton 1979, Howarth 1983. Fits illustrating our results are presented in Fig. 2. Based on our spectral analysis and comparison with the line profiles of HD 108 and HD 191612, we conclude that macroturbulence is likely the dominant line broadening mechanism. Tests performed using He i λ4713 and C iv λλ5801/5812 reveal that convolution of our synthetic spectra with a Gaussian profile (mimicking isotropic macroturbulence) gives better results compared to pure rotational broadening. In practice, a macroturbulence of about 50 km s −1 gives satisfactory fits. However, many profiles -in particular those of He i -have a strong (wind-induced?) asymmetry, and those that have only small asymmetries have broadening that is dominated by 'turbulence'. The projected rotational velocity is an important parameter needed to constrain the stellar and magnetic geometry. To place a stronger constraint on the v sin i of HD 148937, we focused mainly on the C iv λ5801 line as it is reasonably strong, and reasonably symmetric.
We used a 40 kK, log = 4.0 TLUSTY (Lanz & Hubeny 2003) model for the intrinsic line profiles, slightly scaled to match the observed line strength. The model profiles were rotationally broadened with a simple convolution and a linear limb-darkening coefficient of 0.4 (numerical integration of specific intensities gives results that are indistinguishable); isotropic Gaussian turbulence; and a Gaussian instrumental response, before remapping to the sampling rate of the observations, and adding matching noise. Comparison of models & observations was performed for a range of v sin i and turbulence, in both wavelength and Fourier space.
We find that the broadening is roughly described by the constraint that Table 2. Log of CFHT MiMeS observations of HD 148937, including results of the magnetic analysis described in Sect. 3. S/N is peak signal-to-noise ratio per 1.8 km s −1 pixel in the reduced spectrum (occurring around Hα). Phase is computed using the ephemeris expressed by Eq. (2). z = B z /σ is the detection significance of the longitudinal magnetic field. where ζ is the FWHM (2.354σ) of the Gaussian turbulence. Interestingly, this is essentially identical to the result for HD 191612 (Howarth et al. 2007). This is remarkable, given the large difference in their inferred rotational periods. As suggested by , this supports the view that rotational broadening is a negligible contributor to the observed line widths. The upper limit to v sin i from C iv λ5801 is 60 km s −1 , and v sin i = 0 provides an acceptable fit (with ζ = 106 km s −1 ). This is illustrated in Fig. 3. We also explored a variety of weaker lines, but even with such high S/N spectra we are limited by S/N issues. Nevertheless, in the end we concluded that the N iv λ5226 line gives a slightly better upper limit on v sin i of about 45 km s −1 .
These numbers are compatible with those reported by Nazé et al. (2008a). The properties of HD 148937 as derived in this section are summarised in Table 1.

SPECTRAL VARIABILITY
As described by e.g. Nazé et al. (2008a), the complex spectrum of HD 148937 exhibits many similarities to HD 108 and HD 191612, as well as some important differences. Like HD 108 and HD 191612, the C iii λ4650 lines are in emission, with strength    The left-hand panel shows that a rotationally broadened model matching the line width fails to match the line shape; an additional broadening mechanism is required. The remaining panels show that the projected rotational velocity is poorly constrained, because of this significant turbulent broadening.
comparable to the nearby N lines. However, according to Nazé et al., unlike those stars, the C iii lines are not observed to vary. In fact, due to the substantially higher S/N of the new CFHT data, we detect weak variability in a large number of spectral lines, including the C iii lines. Following Nazé et al. (2008aNazé et al. ( , 2010 we characterise line variability using the equivalent width (EW). Before measuring the EW, each spectral line was locally re-normalized, and a telluric correction algorithm was applied to regions redward of 5790 Å. The EWs were then obtained by numerically integrating over the line profile.
The 1σ uncertainties were calculated by propagating the individual pixel uncertainties in quadrature. Neither the Coralie nor the FEROS spectra had pre-computed uncertainties; therefore, a single uncertainty value was estimated from the RMS scatter in the continuum regions around the line profile and assigned to each pixel.
The line showing the strongest variability is Hα, with a variability amplitude of 20-30%. The variations of other photospherically-important lines showing distinct emission contributions -Hγ, He ii λ4686, Hβ, He i λ5876 -are also statistically significant. He i λ6678 is an exception -this line is in full emission, Figure 5. Phased EW variations of lines in the ESPaDOnS spectra of HD 148937: Balmer lines Hα, Hβ and Hγ; He ii λ4686, He i λ5876, He i λ6678; He i λ4471, He i λ4713; C iv λ5801, C iv λ5812; C iii λ4650 lines; and the DIB located at 5797 Å. Symbols/colours distinguish different cycles. Dashed curves represent least-squares sinusoidal fits to the combined ESPaDOnS data. but varies only weakly (the variation is detected only in the higher-S/N ESPaDOnS measurements, and then only marginally). These results are consistent with the report of Nazé et al. (2010).
Other lines of He i (e.g. λ4471, λ4713), while asymmetric, generally show no distinct emission. These lines show no significant variability. The C iv λ5801, λ5812 lines are slightly variable. This appears to be the first detection of EW variability of the optical C iv lines in an Of?p star.
We performed period searches using a Lomb-Scargle algorithm. As a first step, the procedure was applied to each of the timeseries of EW measurements of the strongly variable emission lines. We then used the method described by Folsom et al. (2008), treating each pixel within the line profile as an independent timeseries of measurements. Each line profile was interpolated onto the wavelength grid corresponding to one of the ESPaDOnS spectra. Periodograms were constructed for individual pixels in the profile. These periodograms were then averaged together, weighted by the amplitude of the variability in their respective pixels. The resulting weighted mean periodogram for the profile was then examined for a best fit period.
All periodograms of the emission line EWs show significant power near 7d. These periods, which are reproduced in the pixel-by-pixel periodograms with somewhat larger uncertainties, are all consistent with the best-fit Hα period derived by Nazé et al. (2008a). To derive a single period characterising the variability of HD 148937, the periodograms resulting from the analysis of the EWs and line profile variations of the individual lines were all averaged. The final resulting periodogram, shown in Fig. 4, exhibits a single acceptable period at 7.032 ± 0.003 d. This period is consistent with that reported by Nazé et al. (2010). Based on these results, for the remainder of the current study we adopt the following ephemeris for HD 148937: where we have adopted a zero point corresponding to maximum Hα emission as derived from a sinusoidal fit to the entire set of phased measurements (i.e. ESPaDOnS, FEROS and Coralie).
The phased EW measurements are illustrated in Fig. 5 (for illustrative clarity we show only the higher-quality ESPaDOnS spectra). The variations of all variable emission lines are in approximate phase with the Hα line. Variability with a similar phasing is observed in the C iii 4650 Å lines (which we note are predicted to be blended with an O ii line). Interestingly, the EW variations of the C iv λ5801 line appears to be shifted by approximately 0.25 cycles relative to the variations of the emission lines, although this shift is only marginally significant. We also include EWs measured from the DIB at 5797 Å to illustrate the lack of variability of a reference line that is not formed in the environment of the star.
We find that in many cases the cycle-to-cycle scatter of the EWs is larger than the uncertainties attributable to noise, normalisation uncertainties and telluric line contamination. This appears to be a consequence of evolution of the shapes of the EW variations of some lines. This is especially evident in the Hα and He ii λ4686 lines which exhibit clear changes in amplitude and shape of the EW variation between the two well-sampled cycles (illustrated as green diamonds and blue circles in Fig. 5). The cycle-to-cycle differences in the EW curves of these two lines track each other almost perfectly, and correspond to clear systematic differences in the profile shapes ( Fig. 6) that are expressed in essentially all lines to some degree. As these lines are located throughout the spectrum, these differences cannot be a consequence of systematic errors in local normalisation. Nor, given the location of some lines (e.g. λ4686) well into the blue part of the spectrum, can it be attributed to telluric lines. We therefore conclude that this variability is intrinsic to HD 148937.
A potential explanation of these differences would be an error in the ephemeris expressed by Eq. (2). To check if this could be a plausible cause, we varied the period to see if we could improve the agreement between the EW variations from different epochs. We still find that the 7.032 d period yields the most consistent EW variations amongst different epochs and that changing the period within the uncertainty had a negligibly small effect on the relative phasing of the different epochs. We also compared the profiles from different epochs corresponding to similar EW values but with small relative phase offsets to see if the observed profiles are a better match than those with different EWs but essentially identical phases (e.g. as in Fig. 6). These comparisons show that we find larger discrepancies between the profiles with similar EW with small relative phase offsets than those shown in Fig. 6, confirming that adjusting the ephemeris is not a solution to this phenomenon.
Based on these results, we conclude that the differences in the line profiles are best explained by intrinsic changes in the amount or distribution of emitting material with time, i.e. evolution of the magnetosphere of HD 148937. It is unclear if this evolution is secular or stochastic. This will be discussed further in Sect. 7.

Diagnosis from nightly CFHT spectra
Least-Squares Deconvolution (LSD, Donati et al. 1997) was applied to all CFHT observations. In their detection of the magnetic field of HD 191612 in 2006, Donati et al. developed and applied an LSD line mask containing 12 lines. Given the similarity between the spectra of HD 148937 and HD 191612, we began by using this line list to extract, for all collected spectra, mean circular polarization (LSD Stokes V), mean polarization check (LSD N) and mean unpolarized (LSD Stokes I) profiles. All LSD profiles were produced on a spectral grid with a velocity bin of 36 km s −1 , providing reasonable sampling of the observed mean profile and enhancing the S/N per LSD pixel. The LSD Stokes I profile shows a strong extension to high (∼ 350 km s −1 ) velocities in the blue wing. This asymmetry is similar to that observed in the He i lines of HD 148937, and likely is reflective of the inclusion of such lines in the mask. An additional dip and contribution to the depression in the blue wing is attributed to the DIB located at ∼ 579.8 nm, which was not included in the line mask but which is blended with C iv λ5801 (which is itself included).
Using the χ 2 signal detection criteria described by Donati et al. (1997), we evaluated the significance of the signal in both the Stokes V and N LSD profiles. We obtain a marginal detection (i.e. a false alarm probability (fap) below 10 −3 ) in one spectrum (#1206203). In no case is any signal detected in the N profiles extracted from the individual spectra. If we coadd the LSD profiles of spectra acquired during a single night (typically 2 spectra per night), we obtain one definite detection of signal (false alarm probability fap < 10 −5 ) in the V profiles (for the coadded profile corresponding to spectra IDs 1217692 and 1217696), and one marginal detection (1218069 and 1218073). An examination of this LSD profile shows a weak asymmetry (positive in the red wing, negative in the blue wing) that is not present in N.
To quantitatively characterise the magnetic field implied by our data, from each set of LSD profiles we measured the mean longitudinal magnetic field in both V and N using the first moment method (Rees & Semel 1979) as expressed by Eq. (1) of Wade et al. (2000), integrating from -280 to +180 km s −1 (as did Wade et al. 2011 for HD 191612). The longitudinal field measured from Stokes V is detected significantly (i.e. |z| = |B ℓ |/σ 3) in 3 of our observations. In no case is the longitudinal field significantly detected in N. The results of this analysis are summarised in Table 1.
Due to the relatively large width of the profile, and the limited S/N of our observations, our longitudinal field error bars are relatively large: from about 100 G, with a median of 138 G. As expected, the longitudinal field measured from the N profiles is consistent with the absence of any signal in the N spectrum. The inferred field values are distributed randomly around zero in a manner that is statistically consistent with a Gaussian distribution (as characterised by the detection significance z). The longitudinal field measured from the N profiles exhibits a mean value of 62 ± 26 G (2.4σ). On the other hand, the longitudinal field measured from the V profiles is generally negative, with a (weighted) mean value of −143 ± 26 G (5.5σ). This again suggests the presence of a signal in V that is not detected in N.
As a further test, we have investigated if the distribution of longitudinal field detection significance from Stokes V is in fact significantly different from that derived from N. We employed a two-sided Kolmogorov-Smirnov test. For the cumulative distributions of z V and z N , we obtain D = 0.5625, indicating that the null hypothesis (that the distributions are the same) can be rejected at over 99.9% confidence.
We therefore conclude that the data provide strong evidence of the presence of a magnetic field in the photospheric layers of HD 148937.

Diagnosis from coadded CFHT spectra
Because of the magnitude of the mean longitudinal field (∼ 150 G) and that it appears to remain consistently negative, it is reasonable to assume we are seeing an organised magnetic field, viewed predominantly from a single magnetic hemisphere (i.e. near the nega- Figure 6. Cycle-to-cycle changes of the of Hα, Hβ, He i λ5876, He ii λ4686, He i λ4471, C iv λ5801 lines of HD 148937. Shown are profiles obtained at nearly identical phases according to Eq. (2). Profiles acquired close in time to each other (#1206203 -phase 0.8080 -and #1206207 -phase 0.8123) match perfectly, while a profile obtained 408 days earlier (approximately 58 cycles earlier) at almost identical phase (spectrum #1076968 -phase 0.8008) shows systematic differences unattributable to local normalisation. tive (south) magnetic pole). In this case, even if the star is rotating according to the ∼ 7.03 d period, the Stokes V Zeeman signature's shape would vary relatively little from one observation to the next. We therefore proceed to average all LSD profiles to obtain a single grand average profile, with the hope that we are able to detect a significant Stokes V signature.
The grand average LSD profiles are shown in Fig. 7. The N and V profiles are characterised by an LSD noise level (per 36 km s −1 pixel) of 2.4 × 10 −5 . A clear variation is visible in Stokes V, whereas N shows no variation whatsoever. The detection probability of a signal within the line in Stokes V is effectively 100% (i.e. fap of 1.1 × 10 −8 ), corresponding to a definite detection. No detection is obtained in N (fap of 1.2 × 10 −1 ). The longitudinal field inferred from the grand average V profile is −193 ± 26 G (7.4σ), while no similar field is detected in N (37 ± 26 G; 1.4σ). The detection of a clear signature in the grand average profile provides clear confirmation that HD 148937 hosts a magnetic field, with a strong disc-averaged and phase-averaged longitudinal component.
We have also searched the grand average spectrum for Zeeman signatures in individual line profiles. As illustrated in Fig. 8, marginal signatures are visible in Hβ, He i λ5411 and C iv λ5801, while a rather strong signature is visible in He i λ5876. The shapes of these signatures are compatible with that of the Stokes V LSD profile. Investigation of individual and phase-binned (see below) spectra provided no detection of signatures in individual lines -they are only detectable in the grand average spectrum.
To test the sensitivity of these results to the detailed mask composition, we follow the procedure employed by Wade et al. (2011) for HD 191612, extracting LSD profiles for two additional line masks. We began using a generic line mask based on a 40000 K extract stellar request from the VALD database, using a line depth threshold of 0.1. We sought to construct a mask that yielded a mean LSD Stokes I line profile that exhibited as little variability as possible, and that was as symmetric as possible, while still maximizing Stokes V signal. At each step of the mask development, we visu-alised the agreement of the LSD model (i.e. the convolutions of the Stokes I and V LSD profiles with the line mask) with the reduced spectrum, and evaluated the symmetry and variability of the LSD profile. After removal of the most significantly affected lines, we were left with a mask containing 26 lines, for which we adjusted the mask depths to best match the mean observed line depths. This mask is dominated by lines of neutral He, but also contains lines of ionised He, C, N O and Si. LSD profiles extracted using this mask -which contains more than twice the number of lines in the Donati et al. mask for HD 191612 -yielded LSD profiles nearly identical to those obtained from the Donati et al. mask. This is likely a consequence of the inclusion of many of the same strong He lines in both masks. Measurements of the longitudinal fields extracted from these LSD profiles are in good agreement with those obtained above.

Further tests
As a second test, we extracted LSD profiles using a more restricted mask containing only a half-dozen metallic lines (lines of C, O and Si). This mask yields profiles whose shapes differ substantially from those obtained from the first two masks, as is expected due to the exclusion of the much broader He lines. Longitudinal fields computed using this mask (measured using an integration range of -175 to +115 km s −1 ) are in formal agreement with those obtained using the masks containing He lines, with a mean of −160 ± 23 G. The grand average metallic line LSD profile yields a marginal detection in Stokes V (and no detection in N), and corresponds to a longitudinal field of −121 ± 23 G (5.3σ).
As a third test, we compared the nightly-averaged profiles produced by the metallic line mask to a grid of synthetic Stokes V profiles using the method of Petit & Wade (2011). We restrict ourself to the metallic line mask for this test because in order to treat the LSD profiles as a single line profile for modelling, it is best to use lines with similar shapes. The synthetic Stokes V profiles were produced by a simple algorithm performing a disk integral and assuming the weak field approximation to compute the local Stokes V profiles. The parameters of the synthetic line were chosen to fit the intensity LSD profile. The magnetic model is a simple centered dipolar field, parametrized by the dipole field strength B d , the rotation axis inclination i with respect to the line of sight, the positive magnetic axis obliquity β and the rotational phase ϕ. Assuming that only ϕ may change between different observations of the star, the goodness-of-fit of a given rotation-independent (B d , i, β) magnetic configuration can be computed to determine configurations that provide good posterior probabilities for all the observed Stokes V profiles, in a Bayesian statistic framework. Notice that in order to stay general, we do not constrain the phases nor the inclination. The Bayesian prior for the inclination is described by a random ori- Figure 9. The probability density functions marginalised for B d have been normalized by their peak values in order to facilitate graphical representation. The parameter evaluation for the dipole model treat any difference with the model as additional Gaussian noise that is marginalized, leading to the most conservative estimate of the parameters. entation (p(i) ∝ sin(i) di), the prior for the dipolar field strength has a modified Jeffreys shape with a cut at B d = 40 G corresponding to twice the grid step in order to avoid a singularity at B d = 0 G. The obliquity and the phases have simple flat priors. In order to assess the presence of a dipole-like signal in our observations, we compute the odds ratio of the dipole model with the null model (no magnetic field, Stokes V = 0). We also perform the same analysis on the null profiles. As expected, no single observation is strongly in favour of the magnetic model, although most of the observations that do favour the magnetic model are clustered around the same phases. Taking into account all the observations simultaneously, the odds ratio is 22:1 in favour of the magnetic model. For N, the combined odds ratio is 2:1 in favour of the null model. Note that as the case B d = 0 G is included in the magnetic model, the difference between the two models for noise is expected to be dominated by the ratio of priors, i.e the Occam factor that penalizes the magnetic model for its extra complexity. The derived probability distribution functions (PDF) marginalised for the dipole strength are shown in Fig. 9. The most probable dipole strength according to the metallic line Stokes V profiles is about 350 G -approximately 3 times greater than the longitudinal field inferred from those profiles.
Based on this analysis, we conclude that the shapes of the LSD profiles depend on the composition of the line mask. However, all analyses, independent of the line mask used, imply the existence of a magnetic field with a magnitude of the longitudinal component of 100-200 G. Given these results, and for consistency with the analysis of HD 191612 by Wade et al. (2011), we have based our following analysis on the LSD profiles extracted using the 12 line mask of Donati et al. (2006).

Longitudinal magnetic field variation
To study the magnetic field geometry, we model the longitudinal field variation of HD 148937 as a function of rotational phase. To this end, we proceed as follows. First, we adopt the oblique rotator model as the framework in which we interpret the magnetic observations. This implies that the 7.03 d period (and in particular the ephemeris as expressed by Eq. 2) corresponds to the stellar rotational ephemeris. We have binned the LSD profiles in 7 phase bins (phases determined according to Eq. 2) typically spanning Table 4. Longitudinal magnetic field of HD 148937, binned in phase according to Eq. (2). P V and P N are the detection probabilities derived from the LSD Stokes V and N profiles, respectively. z = B z /σ is the detection significance of the longitudinal magnetic field.

Phase
Stokes 0.03-0.04 cycles (but as large as 0.06 cycles), and re-measured the longitudinal field from these averaged LSD profiles. The binned V longitudinal field measurements are consistently negative, and a majority of them correspond to detections at 3σ confidence or greater. In contrast, the associated measurements obtained from the N profiles scatter randomly about zero. (While one of the N longitudinal field measurements is formally significant at over 3σ, no corresponding coherent signal is apparent in the N profile). The resultant uncertainties are reasonably uniform, about 60-75 G. These binned results are reported in Table 4. The phase variation B ℓ (φ) of the longitudinal field is illustrated in Fig. 10. A fit by Least-Squares of a cosine curve of the form B ℓ (φ) = B 0 + B 1 cos(2π(φ − φ 0 )) to the data yields a reduced χ 2 of 0.3, with parameters B 0 = −200 ± 9 G, B 1 = 90 ± 12 G and φ 0 = 0.24 ± 0.15. Therefore, according to Least-Squares, the longitudinal field is significantly different from zero (at over 20σ), and the variation of the field is detected at over 7σ. On the other hand, the phases of the extrema are relatively weakly constrained. The reduced χ 2 of the data relative to the straight line B ℓ (φ) = 0 (the hypothesis of a null field) is 10.2, indicating that the null hypothesis does not provide an acceptable representation of the data. While the reduced χ 2 of the cosine fit to the data is relatively low, it is not unreasonable given the small number of degrees of freedom.
For illustration, we include the 4 measurements reported by Hubrig et al. (2008) and Hubrig et al. (2011) in Fig. 10, for both their full spectrum and hydrogen line results. Those measurements are qualitatively and statistically consistent with the variation we derive.

Surface magnetic field geometry
With the inferred rotational period and derived radius (from Table 1), HD 148937 should have an equatorial rotational velocity v e = 108 km s −1 . Based on the upper limit on v sin i and its uncertainty obtained in Sect. 3 (see Table 1), we obtain i 30 • . We therefore conclude that HD 148937 is viewed relatively close to a rotational pole, a conclusion consistent with the proposal of Nazé et al. (2010).
Adopting i = 30 • (which ultimately will yield a lower limit on the inferred magnetic field strength), we have fit the phase-binned longitudinal field measurements with synthetic longitudinal field variations corresponding to a grid of oblique dipole surface magnetic fields characterised by the obliquity angle β and the polar field strength B d . The best-fit model (according to the χ 2 statistic)  is characterised by B d = −1020 G and β = 38 • . The 1σ uncertainty contours (shown in Fig. 11) permit models with B d in the range −700 G to −1400 G and β from 20 to 65 • . At 3σ the lower limit on the inferred dipole strength is still 450 G, i.e. the presence of a magnetic field remains confident. We have also performed the same procedure using the unbinned longitudinal field measurements, and we obtain results that are statistically identical.
The sum of the inclination and obliquity angles for models able to reproduce the observed variation approximately obeys the general rule i + β ≃ 70 • .
The extrema of the best-fit dipole model occur at phases 0.24 and 0.74, with approximate 1σ uncertainties of ±0.15 cycles. On the other hand, the emission line variations (in particular those of Hα show extrema at phases 0.0 and 0.5. This will be discussed further in Sect. 7.

STELLAR GEOMETRY ACCORDING TO THE Hα MODULATION AND HIPPARCOS LIGHTCURVE
We have also attempted to constrain the stellar geometry using the equivalent width variation of the Hα line. We have used the "toy" model of Howarth et al. (2007) which assumes an oblique, rigidly rotating, geometrically thin, optically thick Hα disc that introduces modulation of the Hα emission due to its varying projected area as the star rotates. The model includes limb darkening and inner/outer disc radii that are specified input parameters.
We have computed a number of grids of models, varying the inner and outer disk radii (inner:outer=1.05:1.10, 2.0:2.5, 1.05:2.5 in units of R * ) and linear limb darkening (from 0-1), for rotational axis inclination i from 0-90 • and disc obliquity angle α (which we assume to be equivalent to the magnetic obliquity β in the magnetic model) from 0 − 90 • , with steps of 1 • .
For each model we evaluated the χ 2 relative to the observed variation, identifying the 100 best models, which define the locus of best-fit α as a function of i. Figure 12 illustrates the distribution of these best-fit models in the i/α plane.
We find that the results are very insensitive to disc size, and that the modest sensitivity to limb darkening is readily compensated by selecting different pairs of (α/i) values. The models are also highly degenerate with respect to the angles i and α. This is in strong contrast to the results obtained for HD 191612 (Howarth et al. 2007. Why do these models constrain the geometry so poorly compared to 191612? Because the flat-bottomed section of 191612's Hα EW variation gives a strong constraint (the disk must be seen nearly edge-on for a significant range of phase)a factor which isn't available for 148937.
Ultimately, we conclude that for i 30 • (as implied by the upper limit on v sin i and the 7.03 d period) the disc obliquity α 10 • , and for i = 30 • , 10 • α 25 • . These results are compatible with the constraint derived from the magnetic data. Although this modeling is very approximate, it reinforces the applicability of the ORM and supports the interpretation of the emission line variations as due to the varying projection angle of a plasma disc confined to the magnetic equatorial plane.
A second potential constraint on the geometry is potentially derived from the photometric variability. Nazé et al. (2008a) reviewed published analyses of photometric observations of HD 148937, and performed their own analysis of Hipparcos and Tycho photometry of this star. They concluded that no significant change or periodicity is present. Here we compare the observed Hipparcos H p photometric variation, converted to Johnson V-band magnitudes according to Harmanec (1998), with the predictions of Townsend's Monte-Carlo radiative transfer code for simulating light scattering in circumstellar envelopes described by Wade et al. (2011). The models were developed for HD 191612, computed for wind confinement parameter η * (defined in the following section) that differs somewhat from that of HD 148937 (η * ≃ 50 versus 20), but with lightcurves computed for the constraint derived from the magnetic field variation of HD 148937, i + β ≃ 70 • . Unfortunately, the available photometry is not sufficiently precise to allow us to obtain constraints on the geometry. Nevertheless, the predictions demonstrate that high-precision photometry and broadband polarimetry would enable the determination of the individual values of the angles i and β.

DISCUSSION AND CONCLUSIONS
We have presented new spectroscopic and magnetic measurements of the Of?p star HD 148937 based on extensive spectropolarimetric monitoring with the ESPaDOnS spectropolarimeter at the Canada-France-Hawaii Telescope, and spectroscopic monitoring using FEROS at the ESO 2.2m telescope. The observations and their analysis were undertaken within the context of the Magnetism in Massive Stars (MiMeS) Project.
Using the new spectroscopic observations, in tandem with archival UV and IR fluxes, we perform a careful re-determination of the stellar physical properties, and place a slightly improved upper limit on the projected rotational velocity.
Analysis of the new and previously published optical spectroscopic observations confirms the variability of a variety of spectral lines, and detects variability of others for the first time. Most, and probably all, variable lines vary in phase with Hα. (The C iv λ5801 line may appear to vary out of phase by ∼ 0.25 cycles; however, the formal uncertainly of the phases of the extrema is sufficient large that this cannot be robustly established based on our measurements.) The equivalent widths measured from these lines confirm the variation period reported by Nazé et al. (2008aNazé et al. ( , 2010. We find that the variability of the line profiles and equivalent widths is not strictly periodic, and that the absorption and emission line variability exhibits cycle-to-cycle changes. This is tentatively attributed to evolution of the quantity and distribution of emitting material in the magnetosphere of HD 148937. Such a phenomenon has also been reported for the other well-studied Of?p stars (Howarth et al. 2007, Nazé et al. 2008b. Our new magnetic observations of HD 148937 confirm the existence of an organised magnetic field in the photosphere of this Of?p star. The measurements are consistent with a longitudinal magnetic field of constant sign and weak variability. If we interpret the magnetic field and spectral variability within the context of the Oblique Rotator Model (ORM), then the variability period (corresponding to the stellar rotational period) implies an equatorial rotational velocity of 108 km s −1 . To explain the low inferred v sin i, the rotational axis of HD 148937 must be viewed at relatively low inclination, i 30 • . Rotational modulation at low inclination is able to qualitatively and quantitatively reproduce many of the observed properties of HD 148937, in addition to the low v sin i: a longitudinal magnetic field of constant sign and approximately constant strength, the single-wave nature, sinusoidal shape and low amplitude of the Hα EW variation, and the lack of observed photometric variability.
We modelled the longitudinal field variation of HD 148937 according to the dipole ORM, assuming a rotation axis inclination to the line of sight i = 30 • (providing a lower limit on the strength of the (assumed dipolar) surface magnetic field). We derive the dipole polar intensity B d = 1020 −380 +310 G and obliquity angle β = 38 +17 −28 • . For other values of i, the angles follow the approximate relation i + β ≃ 70 • . Using the stellar and wind parameters reported in Table 1, we compute the wind magnetic confinement parameter η * = B 2 eq R 2 /Ṁv ∞ ≃ 20 (neglecting clumping) and rotation parameter W = v eq /v crit = 0.12. This places the Alfven radius at about R Alf = η 1/4 * R * = 2.0 R * . The Kepler (or corotation) radius is located further from the star, at about R Kep = W −2/3 R * = 4.0 R * . As pointed out by ud Doula et al. (2008), for any material trapped on magnetic loops inside the Kepler radius, the outward centrifugal support is less than the inward pull of gravity; since much of this material is compressed into clumps that are too dense to be significantly line-driven, it eventually falls back to the star following complex   Howarth et al. (2007). Different curves correspond to differing disc inner-outer radii and /limb-darkening coefficients. patterns along the closed field loops. Hence, for HD 148937 (as for HD 191612; Wade et al. 2011), all magnetically-confined plasma (i.e. all wind plasma located inside the Alfven radius) is unstable to this phenomenon. HD 148937 therefore appears to lack a region in its circumstellar environment where centrifugal forces are able to support magnetically-channelled wind plasma against infall back toward the stellar surface (although we recall that the uncertainties on η * are quite large, and therefore the location of R Alf insecure). In any case, while the wind confinement parameter of HD 148937 is comparable to those of the other known magnetic Of?p stars HD 191612 and HD 108 and θ 1 Ori C (∼ 50) , its relatively rapid rotation (P rot ∼ 7 d versus ∼ 538 d (HD 191612), ∼ 55 y (HD 108) and 15.5 d (θ 1 Ori C) may lead to a greater influence of rotational centrifugal effects on its wind confinement.
The above discussion implicitly assumes that the variable line emission arises in a reasonably static (in the rotating stellar frame) distribution of plasma located principally in the magnetic equatorial plane. We point out that the He ii λ 4686 and H α λ6563 line shapes could also be interpreted as P Cyg profiles produced in an outflow. If we assume that these profiles arise from a spherically-symmetric wind (admittedly a rather rough assumption given the observed line variability, and the sensitivity of mass-loss rates determined from recombination lines to the (rather uncertain) volume filling factor), we can use the CMFGEN code to model these lines using the physical parameters in Table 1. We find that the profiles of these lines can be reproduced with an average "density"Ṁ/v ∞ of about 1.8 × 10 −9 . This is a of factor five larger than that is needed to fit the UV lines (Ṁ/v ∞ = 3.8 × 10 −10 , neglecting clumping). The width of the emission line profiles however suggests a substantially lower flow velocity in the line formation region (∼ 900 km s −1 ). This results in an only moderately higher mass loss rate (1.6 ± 0.4 × 10 −6 M ⊙ /yr) than inferred by Nazé et al. (1 × 10 −6 M ⊙ /yr). As Hα is usually formed within one stellar radius above the stellar surface (i.e. below the Alfven radius in the present case), the line emission may thus originate from a slower, and thus denser outflow in regions of closed magnetic loops.
One peculiarity of the derived magnetic geometry is its phasing relative to the emission line variations. In the case of the Of?p star HD 191612 , extrema of the Hα EW variation occur simultaneously with phases of magnetic extrema. In the case of HD 148937, however, the magnetic extrema appear to be offset from the Hα extrema by about 0.25 cycles. In the context of a scattering disc located in the magnetic equatorial plane, this phenomenon is difficult to explain. However, the offset is only a 1.5-2σ effect. Given the low amplitude of the longitudinal field variation, we speculate that this is likely an artifact. New, higher-precision measurements are required to clarify this issue.
HD 148937 is distinguished from the other well-studied Of?p stars (HD 108 and HD 191612) by its relatively short rotational period. If we assume that angular momentum is lost through the moment arm generated by the stellar wind and the magnetic field, the characteristic "magnetic braking" spintown time τ spin (Eq. (25) of ud Doula et al. 2009) is about 1.6 Myr (assuming the characteristics of the UV wind, neglecting clumping), and adopting k = 0.17 (ud Doula et al. 2009). The spindown time we compute for HD 148937 is somewhat longer than those of HD 191612 and HD 108 (respectively 0.33 Myr and 0.83 Myr using recent physical and magnetic data for those stars). Given the ages of the Of?p stars (2-4 Myr;ud Doula et al. 2009), these differences alone are insufficient to explain the large range of rotational periods of these stars (assuming they all began their lives with similar rotation rates) assuming similar ages. The short rotational period of HD 148937 could therefore imply that this star is significantly younger than HD 108 or HD 191612.
In closing, we note that HD 148937 is surrounded by the bipolar nebula NGC6164-5. It is known to be enriched (e.g. Dufour et al. 1988) and was therefore proposed to be the result of an giant eruption of the star. A larger "bubble" surrounds the system, and is attributed to wind expansion in earlier phases of the star's life (Fairall et al. 1985). The morphology of NGC6164-5 is axisymmetric, with two bright "blobs" to the north and south of the star. While the origin of the nebula is not well understood, such morphologies are often explained by ejection in the same general direction as the stellar rotation axis,. However, the magnetic results presented here rather suggest that this axis is nearly perpendicular to the plane of the sky. A detailed modelling of the nebular ejections should thus be undertaken, taking into account the geometry suggested by magnetic measurements.
The main conclusions of this study are that an organised magnetic field with a dipole surface intensity of approximately 1 kG is detected in HD 148937. Within the context of the ORM, the Hα period of 7.03 d is interpreted to be the stellar rotational period. This is, by more than a factor of two, the shortest rotational period found so far for a magnetic O star. The stellar radius, rotation period, v sin i upper limit, weak longitudinal magnetic field variation, Hα equivalent width variation and lack of photometric variability can all be reconciled with the ORM if the star is viewed close to a rotational pole (within 30 • ). While the detection of a magnetic field in HD 148937 is now relatively robust, the parameters of the oblique rotator model remain rather uncertain -this is a challenging object to observe, and its longitudinal magnetic field intensity and variability are both weak.
Improved geometrical constraints should be obtained by modeling the Hα equivalent width variation using more sophisticated models, and measuring and modeling the continuum flux and linear polarisation variation with high precision. In fact, first measurements of the Hα line linear polarisation of this star were obtained by Vink et al. (2009), but with no detection of any stellar contribution to the measured polarised flux. (We note however that there was no detection either for θ 1 Ori C, suggesting that higher precision measurements would be desirable). This is the 3rd Of?p star in which a magnetic field has been firmly detected, confirming that the underlying cause of their spectral peculiarities and variations is almost certainly the magnetic field. We therefore propose that the spectroscopically-identified Of?p stars are themselves a class of magnetic stars -the first among the O-type stars -whose distinctive spectral peculiarities and variability identify them unambiguously as magnetic objects. The test of this proposal will be to search for an unambiguously detect magnetic fields in the two newest members of this class, CPD -28 2561 and NGC 1624-2 (Walborn et al. 2010).