Stellar magnetic parameters from Bayesian analysis

# Stellar magnetic field parameters from a Bayesian analysis of high-resolution spectropolarimetric observations††thanks: Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) and at the Télescope Bernard Lyot (TBL). CFHT is operated by the National Research Concil of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. TBL is operated by CNRS/INSU.

Dept. of Geology & Astronomy, West Chester University, West Chester, PA 19383
Dept. of Physics, Royal Military College of Canada, PO Box 17000, Stn Forces, Kingston, Canada, K7K 4B4
E-mail: VPetit@wcupa.edu
Accepted 2011 October 26. Received 2011 October 22; in original form 2011 August 17
###### Abstract

In this paper we describe a Bayesian statistical method designed to infer the magnetic properties of stars observed using high-resolution circular spectropolarimetry in the context of large surveys. This approach is well suited for analysing stars for which the stellar rotation period is not known, and therefore the rotational phases of the observations are ambiguous. The model assumes that the magnetic observations correspond to a dipole oblique rotator, a situation commonly encountered in intermediate and high-mass stars. Using reasonable assumptions regarding the model parameter prior probability density distributions, the Bayesian algorithm determines the posterior probability densities corresponding to the surface magnetic field geometry and strength by performing a comparison between the observed and computed Stokes profiles.

Based on the results of numerical simulations, we conclude that this method yields a useful estimate of the surface dipole field strength based on a small number (i.e. 1 or 2) of observations. On the other hand, the method provides only weak constraints on the dipole geometry. The odds ratio, a parameter computed by the algorithm that quantifies the relative appropriateness of the magnetic dipole model versus the non-magnetic model, provides a more sensitive diagnostic of the presence of weak magnetic signals embedded in noise than traditional techniques.

To illustrate the application of the technique to real data, we analyse seven ESPaDOnS and Narval observations of the early B-type magnetic star LP Ori. Insufficient information is available to determine the rotational period of the star and therefore the phase of the data; hence traditional modelling techniques fail to infer the dipole strength. In contrast, the Bayesian method allows a robust determination of the dipole polar strength,  G.

###### keywords:
stars: magnetic fields– stars: early-type – techniques: polarimetric – methods: statistical.
pagerange: Stellar magnetic field parameters from a Bayesian analysis of high-resolution spectropolarimetric observationsthanks: Based on observations obtained at the Canada-France-Hawaii Telescope (CFHT) and at the Télescope Bernard Lyot (TBL). CFHT is operated by the National Research Concil of Canada, the Institut National des Sciences de l’Univers of the Centre National de la Recherche Scientifique of France, and the University of Hawaii. TBL is operated by CNRS/INSU.11pubyear: 2010

## 1 Introduction

The evolution of a massive star is strongly determined by its rotation, as well as the mass lost through its stellar wind, both of which can be strongly influenced by the presence of a magnetic field (e.g. Townsend et al., 2010; Wade et al., 2011). A magnetic field can furthermore couple different layers of a star’s interior, thereby modifying internal differential rotation and circulation currents (Maeder & Meynet, 2005). If a field has a large-scale component that extends outside the stellar surface, it can also channel the outflowing stellar wind, creating a structured wind - a magnetosphere - which will modify the rate and geometry of mass loss along with its observational properties (Landstreet & Borra, 1978; Shore & Brown, 1990; Babel & Montmerle, 1997a, b). Furthermore, if the field couples the rotating surface of the star with the outflowing stellar wind, both effects will result in angular momentum loss (via the outflowing stellar wind), which differs strongly from that of a non-magnetic star (ud-Doula et al., 2009). As angular momentum and mass loss are cornerstone inputs in stellar evolution calculations, it is crucial that the effect of magnetic fields in massive stars be understood properly. For example, evolutionary tracks and isochrones can be used to interpret large datasets associated with OB associations like the VLT-FLAMES surveys of massive stars (Hunter et al., 2008).

Over the last decade, our knowledge of the properties of massive star magnetic fields has significantly improved, in large part due to a new generation of powerful high-resolution spectropolarimetric instrumentation. Traditionally, stellar magnetic fields were modelled using measurements of the longitudinal magnetic field of the star, yielding a single quantity: the strength of the disc-integrated line-of-sight component of the magnetic field (e.g. Bagnulo et al., 2006). In the absence of a field detection, the interpretation of such data was entirely dependent on the (unknown) stellar and magnetic geometry, and therefore highly ambiguous. However, modern high-resolution spectropolarimetry yields Doppler-broadened, velocity-resolved Stokes profiles measured across spectral lines, which encode additional information about the field geometry. These data therefore allow a more robust inference of the field characteristics.

Remarkably detailed information about the strength and local/global structure of stellar magnetic fields can be extracted from high-resolution, phase-resolved spectropolarimetric datasets acquired in two or four Stokes parameters (i.e. Stokes or Stokes ) using techniques such as Magnetic Doppler Imaging (e.g. Kochukhov & Wade, 2010). However, detailed modelling of this sort relies on extensive high-quality datasets that demand significant telescope time, and which can only be obtained for a small number of stars. The lower-quality data that can be obtained for stars with less suitable observational properties can still be approximately investigated using parametric models; nevertheless, even such a simple approach often requires approximately a dozen observations per star.

Large observing programs, such as the Magnetism in Massive Stars (MiMeS) project, have dedicated significant resources to survey large samples of massive stars in search of magnetic fields (Wade et al., 2011). Such surveys typically acquire a small number of Stokes observations (typically 1-3) per star, for a large number of stars. An outstanding problem concerns how to extract useful information about the magnetic field of a star from such a limited data set.

In this paper we describe an approach to constrain the magnetic field strength of a star from small numbers of high-resolution Stokes observations, assuming the dipole oblique rotator paradigm, and expressed using the formalism of Bayesian statistics. In Sect. 2 we describe briefly the model used to synthesise the emergent local circular polarisation produced by the Zeeman effect under the weak-field approximation. We present the disk-integrated emergent Stokes profiles obtained for a dipolar magnetic topology of a rotating star. In Sect. 3 we introduce the Bayesian algorithm used to compare a set of high-resolution Stokes observations with a grid of synthetic line profiles. Sect. 4 presents the application of the Bayesian algorithm to simulated data. We demonstrate that the Bayesian odds ratio can help detect a weak magnetic signal, by quantifying which target should be re-observed. It can also discriminate, with a few additional observations, whether an observation has a noise pattern that by chance looks like a magnetic signal versus a real magnetic signal. We also show that the Bayesian algorithm provides a meaningful upper limit on the magnetic strength in the case of a non-detection, and is able to reliably estimate the dipole field strength in the case of a detection. Some limited constraints can also be obtained for the geometry of the dipole. We compare the Bayesian diagnostics with the traditional global longitudinal field and signal detection probability diagnostics. Finally, in Sect. 5, we present the Bayesian results for the magnetic B-type star LP Ori (Petit et al., 2008), obtained with observations of various signal-to-noise ratios.

## 2 Line synthesis

### 2.1 Local profiles

Splitting of spectral lines due to the Zeeman effect corresponds to about 1 km s per kG of surface field modulus. This implies that even relatively strong fields are difficult to detect reliably from Zeeman splitting in the spectrum of most stars. The situation is particularly challenging in the case of hot stars, where thermal broadening (of order a few km s), turbulent broadening (of order a few tens of km s) and rotational broadening (potentially a few hundreds of km s) combine to fully obscure any modification of the line profile due to the Zeeman effect.

The Zeeman effect provides us with a second useful tool, in the form of the polarisation of the Zeeman components. Zeeman components exist in two different types: components, spread symmetrically about the zero-field wavelength of the spectral line and components, with symmetric wavelength offsets to the red and blue of the zero-field wavelength. If the external magnetic field is oriented parallel to the observer’s line of sight (a longitudinal field), the components vanish, while the two groups of components have opposite circular polarisations. In a field aligned perpendicular to the line of sight (a transverse field), the components are linearly polarised perpendicular to the field direction, while the components are linearly polarised parallel to the field direction. Therefore, spectropolarimetric observations are an extremely useful tool for detection and characterisation of both the strength and orientation of stellar magnetic fields.

In order to correctly predict the polarised spectrum that will emerge from a stellar atmosphere, one must perform radiative transfer in an anisotropic medium, where the propagation properties depend on the propagation directions. In the weak-field regime, we can treat the anisotropic absorption and dispersion profiles as perturbations of the isotropic case (Landi degl’Innocenti & Landolfi, 2004). In that case, at first order, only the longitudinal component of the magnetic field contributes to the polarisation, and therefore only circular polarisation is predicted. The contribution of the transverse field, and the occurrence of linear polarisation, is a second-order effect. Therefore, circular polarisation (Stokes spectra) signatures are generally stronger, and used for large surveys as such observations provide the best detection threshold.

Although there are some elaborate codes that are able to accurately synthesise the detailed circular polarisation profile of an arbitrary spectral transition (for a review see Wade et al., 2001), we will use the weak-field approximation, in order to draw a simple picture of the Stokes spectrum emerging from a star. The 1st order solution for the circular polarisation to the polarised radiative transfer equations is given by:

 V(v)=−e4πmc2cgeffλ0B||dI(v)dv. (1)

The Stokes profile emerging from a point at the surface of a star (referred to as a local profile) has the same shape as the derivative of the local intensity profile , scaled by the longitudinal component of the magnetic field at that point, by the wavelength of the transition, and by the effective Landé factor of the transition. The effective Landé factor represents the separation of a triplet Zeeman pattern that approximates a more complex Zeeman pattern, and can be calculated under LS coupling, or taken from atomic experiments. Although we will be using the weak field approximation throughout this paper, the Bayesian algorithm can use synthetic profiles computed with any spectrum synthesis code.

A multi-line approach will usually be applied – in order to increase the signal to noise ratio (s/n) – to the data for which a Bayesian approach will be useful. Although these multi-line techniques have been developed and refined for more than two decades now (e.g. Borra et al., 1981; Semel & Li, 1996), the most widely-used approach is the Least Squares Deconvolution (LSD) method introduced by Donati et al. (1997). Under the weak field approximation, the LSD method assumes that all the spectral lines have the same shape, and that this shape is scaled by the line depth for the intensity line profiles (Stokes ) and by the product of the line depth, the wavelength and the effective Landé factor for the circularly polarised line profile (Stokes ). The LSD method therefore already approximates the Zeeman pattern of a line by a triplet whose separation is given by an effective Landé factor. For a complete discussion on the circumstances under which a LSD profile can be treated as a single spectral line, see Kochukhov et al. (2010).

### 2.2 Disk integration

In order to predict the Stokes spectrum of a given stellar spectral line, we must take into account the contribution from every point on the visible hemisphere. In this paper, we consider a single spectral line of arbitrary depth and a width corresponding to the sum of the spectral resolution and an ad hoc local broadening width. The contribution of each point on the visible hemisphere is dictated by the effective surface area and a wavelength-independent limb-darkening of the form:

 IcIc,0=1−ϵ+ϵcosθ, (2)

where is the intensity at the stellar disk centre, is the angle between the normal to the surface and the line of sight, and is the limb-darkening factor (Gray, 1992). We adopt for the rest of this paper. The surface velocity field of the stellar disk is defined by the projected equatorial velocity assuming rigid rotation.

We consider a simple magnetic topology with a dipolar field of strength , as the vast majority of intermediate mass and high mass stars’ magnetic fields appear to be predominantly dipolar (e.g. Borra & Landstreet, 1980; Aurière et al., 2007). The orientation of the dipole with respect to the observer’s reference frame is described by the inclination of the rotational axis to the line of sight, the rotational phase and the obliquity of the magnetic axis with respect to the rotational axis.

Figure 1 shows the radial field (top) and the longitudinal field (bottom) for a dipole field seen positive pole-on (left) and magnetic equator-on (right). The radial field is defined by its orientation with respect to the stellar surface normal. For example, when we are looking at the positive pole (left), the radial field is at its maximum (red) in the middle of the stellar disk. On the edge of the disk, the field is parallel to the surface, and the radial field is null (green).

The longitudinal field is the important quantity for the polarised radiative transfer. The longitudinal field represents the projection of the magnetic field vectors with respect to the observer’s line of sight (perpendicular to the paper plane). If we look at the magnetic field seen pole-on (bottom left), all of the magnetic vectors that are oriented toward the observer are colour-coded in red shades, green corresponds to a null longitudinal field component, and magnetic vectors oriented away from the observer are colour-coded in blue shades. As seen in Eq. (1), two longitudinal magnetic fields of same magnitude but opposite signs will produce inverse local Stokes profiles. However, a good fraction of the field components oriented away from the observer are located on the hidden hemisphere of the star. There is therefore a net positive longitudinal component on the visible hemisphere (i.e. the visible hemisphere is more red than blue). In terms of the total emergent Stokes , positive local profiles (local profiles for which the longitudinal field is positive) will dominate the total line profile. This effect is even more pronounced when limb darkening is taken into account, as the edges of the disk contribute less to the total brightness.

We now turn our attention to the case where we are looking at the magnetic equator (bottom right). In principle, there are as many magnetic vectors pointing toward the observer as away. The global longitudinal component is therefore null, as every magnetic vector on the left side of the stellar disk has its opposite on the right side of the disk. In terms of the emerging Stokes profile, each positive local Stokes profile will be cancelled out by the negative local Stokes profile of same amplitude coming from the opposite side of the visible stellar disk, even when considering limb-darkening. The resulting Stokes profile should therefore be null when the global longitudinal component is null. This is effectively the case, if the star is not rotating.

If the star is rotating around a rotation axis oriented toward the top of the paper (as represented by the black grid), the line of sight component of that rotation will introduce some Doppler shifts to different points on the stellar surface. For example, if the star rotates from left to right, the point situated to the extreme left of the stellar disk will be shifted by . The point situated at the extreme right of the stellar disk, whose local Stokes profile would in principle cancel out the one produced by the leftmost point, will now be shifted by . Therefore, the rotation is able to separate in the velocity space the local Stokes line profiles that would otherwise cancel out, and a net circular polarisation profile can be seen even if the global longitudinal field is null. It has been shown that a as low as a few km s is enough for instruments able to resolve this effect in the line profile (Wade et al., 2000b).

Figure 2 shows an example of the shapes and relative amplitudes of Stokes profiles emerging from a star rotating at 50 km s (top). The radial and longitudinal fields are shown (middle), as well as the global longitudinal field curve (bottom). The phases indicated by red dots on the curve indicate the rotational phases corresponding to each shown profile. We can see that for this dipole configuration (, ), the emerging Stokes profile has an amplitude as strong when the global longitudinal field is null (2nd and 4th phases) than when it is at its maximum (1st, 3rd and last phases). We note the Stokes profiles at  G are symmetric around the centre of the line. If we were to observe such profiles with an instrument that does not resolve the line profiles, we would indeed obtain a null global longitudinal field, even if the profiles are as clearly detectable with high-resolution observations than when the longitudinal field is at its maximum.

Figure 3 shows a second dipole configuration (, ), where the rotation is not able to separate components that will cancel out. When the longitudinal field is null (middle profile, we are looking at the magnetic equator), we can see that as the star rotates from left to right, all the points situated on the left side of the stellar disk will be shifted to the blue. However, each half of the visible hemisphere contains as many field vectors oriented toward the observer as oriented away from the observer. Therefore, all the red-shaded points situated on the left side of the stellar disk will have a corresponding inverse profile that will be Doppler shifted by the same amount, and the cancellation will occur. This particular symmetry occurs when , at a phase where we are looking at the magnetic equator.

Due to the intrinsic symmetry of a dipole field, as well as the fact that circular polarisation is only sensitive to the longitudinal field components, some parameter degeneracy is encountered for the resulting Stokes profiles. For example, if a dipolar field of parameter , and produces a profile , the following symmetries (or anti-symmetries) occur:

 ⎛⎜ ⎜ ⎜⎝90∘−i90∘−βφf⎞⎟ ⎟ ⎟⎠,⎛⎜ ⎜ ⎜⎝90∘−iβ180∘+φ−f⎞⎟ ⎟ ⎟⎠,⎛⎜ ⎜ ⎜⎝i90∘−β180∘+φ−f⎞⎟ ⎟ ⎟⎠. (3)

## 3 The bayesian algorithm

Bayesian statistics have proven to be an useful tool to find plausible models to explain astronomical data (e.g. Tuomi et al., 2011; Arregui & Asensio Ramos, 2011). A probability density distribution provides a quantitative encoding of the plausibility of a certain hypothesis (), given our current state of knowledge (). The main pathway to Bayesian statistics (Jaynes, 2003) is the Bayes theorem:

 p(Hi|D,I)=p(Hi|I)p(D|Hi,I)p(D|I). (4)

This theorem states that the posterior probability density of an hypothesis, given a new dataset and the current state of knowledge, is equal to the product of prior knowledge of the plausibility of the hypothesis and the likelihood of obtaining the new dataset if the hypothesis is true. The global likelihood, , acts as a normalisation constant. The act of summing (or integrating in the case of a continuous hypothesis set) is called marginalisation.

One difficulty of modelling the Stokes spectra of stars observed as part of a survey like MiMeS is that most of the time the rotational phase of the star at the moment of the observation is not known, because the rotational period is generally unknown. Our approach is therefore to use Bayesian probability nomenclature to find which set of phase-independent configurations can reproduce the observations. In other words, while the geometry of the field must stay the same for each observation of a given star, the phase may in principle have any value.

The hypothesis to be tested is the presence of an oblique dipolar magnetic field in a particular star. The predictions of this hypothesis are represented by the model , parametrized with the parameters (=[, , ]) and (=[..]), the latter representing a set of phases associated with a set of Stokes observations (=[..]). The Bayes theorem tells us that the joint posterior probability density for the parameters, assuming the veracity of the model , is

 p(B,Φ|D,M1)=p(B,Φ|M1)p(D|B,Φ,M1)p(D|M1). (5)

Our prior knowledge of the probability density for each model parameter, which can be as simple as its expected range, can be encoded in the prior term , whereas the information brought forward by the new observations is reflected in the likelihood term . The global likelihood is the normalisation term , which is equal to the total probability:

 p(D|M1)=∫∫p(B,Φ|M1)p(D|B,Φ,M1)dΦdB. (6)

We then treat the set of phases as nuisance parameters by marginalising the joint posterior probability:

 p(B|D,M1)=∫p(B,Φ|D,M1)dΦ. (7)

We therefore ensure that an adequate possesses Stokes profiles at some phases that match all the observations. From this posterior probability density, we can then explore the plausible region of the parameter space of an oblique dipole field given the data, assuming that the model is true.

In the case of a limited number of spectropolarimetric observations, we will be generally interested in the field strength value, as only weak constraints can be placed on the geometrical parameters. This is particularly relevant in the case of a non-detection, where we wish to put an upper limit on the strength of an undetected dipole field. To obtain the posterior probability density for a given parameter, we need to marginalise the joint probability over the other parameters. For example, the posterior probability density marginalised for the field strength is:

 p(Bp|D,M1)=∫p(B|D,M1)didβ. (8)

Another powerful application of Bayesian statistics is the ability to quantitatively test the plausibility of one hypothesis versus another. We can therefore compare the plausibility of the oblique dipole model , by computing the so-called odds ratio of this model compared to the model representing the absence of a magnetic field. To do so, we need to compute the posterior probability of a given model which is, according to the Bayes theorem:

 p(Mi|D,I)=p(Mi|I)p(D|Mi,I)p(D|I). (9)

The prior term encodes the plausibility of the model, given our current knowledge, and the global likelihood encodes the plausibility of the model, given the new observations.

Therefore, the odds ratio of our two competing models, and , can be written as:

 odds (M0/M1)=p(M0|D,I)p(M1|D,I)=p(M0|I)p(M1|I)p(D|M0,I)p(D|M1,I). (10)

Typically, as no model is preferred prior to the acquisition of Stokes observations, the ratio of priors will be set to unity. The global likelihood of the dipole model is given by Eq. (6). As the model representing the absence of magnetic field has no parameters, its global likelihood is simply the product of the likelihoods that for each observation. According to Jeffreys (1998), the evidence in favour of a model is considered weak when the odds are (3:1), moderate when (10:1), strong when (30:1) and very strong when (100:1).

### 3.1 Practical implementation

Given the complexity of disk integration for a rotating magnetic star, it is not practical to solve the Bayesian problem analytically (for an example applied to the local Stokes profiles of the Sun, see Asensio Ramos, 2011). A Markov-chain Monte-Carlo method works by choosing series of parameter sets, the parameter regions with a high likelihood being more likely to be picked. The sampling itself therefore reflects the posterior probability density. This method is not well suited to this problem because a large number of calculations are required to assess the probability of each configuration, due to the marginalisation over all possible phases. We therefore chose a numerical, brute-force approach in order to explore the parameter space by sampling it with a regular grid.

When evaluating Eqs. (5) and (7), we can make use of the fact that the prior probability densities for each parameter are independent such that . Furthermore, the prior probability for the phases are the same such that . Finally, the likelihood can be factored into given that the likelihood of one observation is only dependent on the phase . Therefore, Eq. 7 can be written as:

 p(B|D,M1)=p(B|M1)∏Nn=1∫p(φ|M1)p(Dn|B,φ,M1)dφp(D|M1). (11)

For each point in a [, , , ] parameter space, we can compute the likelihood between the synthetic profile produced by these parameters and each observation. In fact, given the Stokes profile symmetry properties described in Eq. (3), only the synthetic profiles for a quarter of the parameter space need to be computed. Furthermore, synthetic profile interpolation is possible in the direction, as the local Stokes profile amplitude increases linearly although the shape remains the same. The maginalisation integrals are performed numerically using a five point Newton-Cotes algorithm. An adequate sampling of the parameter space was chosen to obtain accurate values of .

### 3.2 Probability densities

For each parameter of the oblique dipole model, a prior probability density is needed to evaluate Eq. (5). The dipole field strength is a parameter that can vary over several decades (from a few dozen G to a few kG). We therefore used a Jeffreys prior, which sets an equal probability per decade and is therefore scale invariant. We used a modified form (Gregory, 2005a) in order to eliminate the singularity at :

 (12)

This prior function behaves as a flat prior when , and as a Jeffreys prior when . The maximum dipolar field strength is adjusted according to the strength of the Stokes signal and by the quality of the data (for example, if the s/n is lower, a field of a higher strength can hide in the noise in the case of a non-detection). For this paper, the grid has 250 values. A finer sampling might be required when using a large number of observations, in order to avoid under sampling the joint posterior probability density. We tested the effect of the choice of the parameter, and we settled using a value about twice the grid step.

The inclination of the rotation axis to the line of sight is a “position” parameter (invariant under a shift of zero position), for which a flat prior would be well suited. However, we know that if the orientations of the rotation axes of stars are generally randomly distributed, the probability that the inclination angle is between and is:

 p(i|M1)=12sini. (13)

The inclination could in principle vary from to . However, if a non-null is measured, it is possible to put a lower limit on the inclination by estimating the break-up velocity – typically around 700 km s for a massive OB star. We used a grid of 37 values, giving a sampling at least every .

The same reasoning could apply to the obliquity angle of the magnetic axis to the rotation axis. However, we ignore at this point if the magnetic axes are generally randomly distributed, or if a general relation exists for the position of the magnetic axis relative to the rotation axis. For example, Landstreet & Mathys (2000) have shown that the magnetic axis of long-period ApBp stars is generally aligned with their rotation axis. Therefore, instead of using a prior that assumes a randomness of the magnetic axis position, we kept a flat prior for the obliquity angle:

 p(β|M1)=1βmax−βmin. (14)

The obliquity varies from to , sampled every .

A flat prior probability is also used for the rotational phase:

 p(φ|M1)=1φmax−φmin, (15)

with and , sampled every .

The likelihood function of a given data set corrupted with Gaussian noise is given by:

 p(Dn|B,φ,M1)=(2π)−M2[M∏k=1σ−1k] exp{−12M∑k=1(dk−fk)2σ2k}, (16)

where represents one of the datapoints with its variance , and represents the model prediction.

When performing a parameter estimation, assuming the veracity of a model, it is a good practice to treat the variance of our data as a model parameter itself, using our estimate of the variance (i.e. the error bars) as a starting point. Proceeding in this way, the final results will be less sensitive to our estimation of the variance, as noise propagation can be problematic in heavily processed data, as can be the case for LSD profiles. Furthermore, this approach treats any features that cannot be explained by this particular model (for example higher order polar field components or distortions of the profile due to abundance spots) as additional noise. According to the maximum entropy principle, a Gaussian distribution will be the most noncommittal about information we do not have, leading to the most conservative estimates of the parameters. If our estimate of the variance of each point is denoted by , we introduce a “noise scaling” parameter (Gregory, 2005b), in order to estimate the total noise variance :

 σ2k=s2kb, (17)

where varies around unity. The resulting expression for the likelihood is:

 exp{−b2M∑k=1(dk−fk)2s2k}. (18)

The resulting posterior probability density will therefore be marginalised over . This procedure is only used for the parameter estimation of model . The noise scaling parameter is a scale parameter, and a Jeffreys prior will be used:

 p(b|M1)=1bln(bmaxbmin). (19)

The noise scaling parameter for the parameter estimation varies from to , corresponding to an underestimation of the variance by a factor of 10 and an overestimation of the variance by a factor of 2, respectively. This conservative range can be verified a posteriori, by the probability density marginalised for .

## 4 Numerical tests

In order to demonstrate our Bayesian method, we have simulated some realistic synthetic observations, and analysed them using our procedure.

### 4.1 Ideal non detection

We start with a simple test representing an ideal non-detection. The local intensity line profiles are represented by Gaussians, with the width of the point spread function of an instrument with spectral resolution of 4.2 km s (R=65000), typical of current high-resolution spectropolarimeters like ESPaDOnS, and an extra turbulent broadening of 5 km s, typical of the broadening encountered in hot stars. We used a shallow line with a depth of 0.17 . The rotational velocity field is set to  km s. The associated circular polarisation is set to zero (i.e. no magnetic field) and the error bars are set to represent a certain s/n. We used two simulated observations corresponding to a s/n of 12 000 and 24 000 respectively, values that can typically be achieved for the LSD profiles of hot stars. The effective Landé factor and the wavelength of the spectral line are set to 1.2 and 5 000 Å respectively. We performed the Bayesian analysis on a grid extending up to 1 kG, with a maximum rotation velocity of 700 km s, leading to a minimum inclination of .

Figure 4 shows the probability densities resulting from our Bayesian analysis, for the s/n of 12 000. The three bottom panels show the posterior probability densities marginalised for each parameter – , and , from left to right respectively. The densities have been normalised by their maximum, in order to facilitate the display. The three top panels show the 2D posterior probability densities, marginalised for the -, - and - planes. The contours encircle the regions that contain 68.3%, 95.4%, 99.0% and 99.7% of the probability.

As it can be seen from the probability density distribution, the bulk of the probability is concentrated at low values. The amplitude of the Stokes profiles produced for low values will remain under the noise level, for all possible dipole orientations. As increases, the amplitude of the Stokes profiles will progressively grow over the noise level, and only a restricted number of possible orientations will produce Stokes profiles with amplitudes below the noise level, as explained in section 2.2. The probability density marginalised for therefore has a decreasing exponential-like shape, with a tail extending far away from the mode of the distribution.

The shape of the 2D - and - planes gives us insight about the geometries that are more likely to produce a non-detection for a large field strength value. When we look at a star rotational equator-on (), the obliquity required to obtain a low amplitude Stokes profile is around or . Thus, as the field is nearly aligned with the rotational axis, the amplitude of the Stokes profiles will stay small for the whole phase range. For lower inclination angles, the Stokes signal will vanish only when the dipole field is seen nearly magnetic equator-on, therefore such an null observation is only possible on a restricted interval of rotation phases (as it was shown in Figure 3). Consequently, if we assume that a strong field is present but is not producing any Stokes signal, it is more likely that we are observing a dipole configuration for which the Stokes profiles always have a small amplitude, rather than a dipole configuration for which the Stokes signal vanishes for only a small fraction of the rotational period. This explains why the probability density is more extended toward higher values for and .

The 2D - plane shows that we have no constraint on the inclination and obliquity angles. The shape of the probability density marginalised for is in part due to the extra probability at higher for , but is mainly driven by the prior probability that favours large inclinations (as given in Eq. (13)). Following the same reasoning, the aligned dipoles ( or ) are slightly favoured.

Table 1 compiles the base ten logarithm of the odds ratio computed for each observation taken individually, and for the two observations combined together. In order to illustrate the effect of the priors, we also computed the odds ratios while using a flat prior for all the parameters. For a single low s/n observation, the model for the absence of a magnetic field is preferred by a factor of 2. When doubling the s/n, the odds in favour of only increase to a factor 2.8. When the two observations are combined together, the odds ratio goes up to a factor of 3.

If we had used a flat prior for all parameters, the odds ratios would have been in favour of the model by a factor of roughly 10 and 20, for the low and high s/n respectively. Jeffreys priors are used when we ignore the exact scale of a parameter. This serves to balance the high scales and the lower scales. In a non-detection case, the bulk of the probability is situated at low values. With a Jeffreys prior, the small scales have as much weight as the larger scales, hence the configurations at high values that produce bad fits dominate less the global likelihood. The flat prior effectively gives more weight to the larger scales, hence rejecting more strongly the model. We note here that the model is in fact a subset of the model, given that the field strength range extends down to 0 G. This means that both models can reproduce this simulated dataset perfectly. However, the model is penalised by its complexity, which is mathematically encoded by the priors. Therefore, for a non-detection, the odds ratios will never be very strongly in favour of the model by a large number, as the discrimination comes mainly from “Occam’s razor”. Furthermore, the posterior probability density is sensitive to the exact choice of the prior in this case. Given the decreasing exponential behaviour of the probability density, we think that the more conservative Jeffreys prior is more suitable, as it is less eager to reject the model.

Assuming that the dipole model is true, we can perform a parameter estimation to determine which dipole strengths would be admissible by our observations. This can be achieved using the marginalised probability density for . Integrating this probability density between two values gives the probability that the true value of the field strength lies between these two values. The probability density itself can be used in many applications, such as building a statistical field strength distribution for a stellar population.

In other applications, such as the study of magnetically confined wind shocks, it would be valuable to get an upper field strength limit. When a probability density has a Gaussian-like shape, it is customary to express our confidence intervals by regions enclosing a certain percentage of the total posterior probability, namely 68.3%, 95.4% and 99.7%. With a Gaussian distribution, the 95.4% region will be twice as extended as the 63.8% region, and the 99.7% will be 3 times as extended as the 68.3% region, analogous to the 1, 2, 3 contours in frequentist statistics. However, as our probability distribution does not have a Gaussian shape, but is rather shaped like a decreasing exponential, these credible regions will not be regular (i.e. the 99.7% credible region will reach much farther than 3 times the 68.3% credible region, taking into account the extended tail of the distribution). We therefore added a credible region enclosing 99.0% of the probability.

Table 2 gives the upper limits of the credible regions (the lower limit being 0 G in each case) for each percentage threshold of the probability. For example, the probability that the field strength is lower than 258 G is 99.0% for the low s/n case. For the higher s/n case, all the credible regions are narrower. When we combine the two observations together, the 68.3% and 95.4% regions are similar to the high s/n case, but the probability density is more peaked (i.e. it has a less extended tail toward the high values), as can be seen from the narrower 99.0% and 99.7% credible regions (110 G and 196 G respectively). When combining two non-detection observations, the probability density of high values narrows around and . The likelihood of the , configurations stays the same, as none of the phases produce any Stokes signal. However, the likelihood of the other degenerate configurations decreases, as it becomes less likely that we have observed the star both times during the same narrow range of phases that do not produce a Stokes signal, if the observations were taken at random times. Therefore, the high- tail becomes less important as we combine multiple observations.

Table 2 also gives the credible regions for the combined observations, while considering a flat prior for each parameter. In that case, the shape of the probability density is more weighted toward the high scales, and the credible regions are more extended toward the high values.

We also give the credible regions when considering a fixed variance of our data. However, in this particular case, as the observations can be reproduced perfectly by both models (as we set Stokes strictly to zero), the inclusion of the noise scaling parameter does not change the probability density much, except for a slight tightening of the credible regions as we allowed for the possibility of variance overestimation.

### 4.2 Ideal detections

In order to show the behaviour of the algorithm in the presence of a magnetic signal, we now add non-zero Stokes profiles to the simulated observation corresponding to a s/n of 12 000. Once again, no random noise was added, to illustrate the ideal case. To draw a connection between the credible regions found in the preceding section, we chose field values corresponding roughly to the upper limits of these credible regions: 30 G, 100 G, 250 G and 450 G. The magnetic configuration was set to and , at the rotational phase when we are looking straight at the magnetic pole (). This corresponds to a profile with a shape like the leftmost profile of Figure 2.

Table 3 gives the odds ratios for the used priors, as well as the ones computed with flat priors for comparison. We also give the credible regions containing 68.3%, 95.4%, 99.0% and 99.7% of the total probability. The lower limit is 0 G, unless indicated otherwise.

The posterior probability density for the 30 G observation is indistinguishable from the observation from the previous section. The odds ratio is still marginally in favour of the model, and the credible regions for are similar.

At 100 G, the odds ratio is still in favour of , although both models are now nearly as likely (). The odds ratio computed with a flat prior still rejects the model by a factor of three. The shape of the probability density, as shown in Figure 5, is also different than the non-detection case, with a secondary peak in the probability density marginalised for , although the probability density still peaks at 0 G. The credible regions for are therefore extended toward higher values, and the 68.3% credible region upper limit (123 G) encompasses the real field value. As it was the case for the ideal non-detection, an aligned magnetic configuration is preferred, as it maximises the chances of observing this particular profile shape. The - 2D plane is different from that of the ideal non-detection. The shape of the magnetic signal adds an extra constraint, by rejecting the configurations for which the positive pole is never located on the visible hemisphere.

At 250 G, the odds ratio has switched in favour of the dipole model, by three orders of magnitude (). However, our chosen prior is again more conservative than the flat prior, which favours more strongly (). As a flat prior overweights the larger scales compared to a Jeffreys prior, our used prior needs a more significant likelihood at large values in order to favour the magnetic model.

The posterior probability density is now typical of a detectable magnetic field (Figure 6). A sharp cut in the probability density marginalised for and the 2D planes shows a tight constraint on the lower field limit. The distribution peaks around 275 G and decreases slowly toward higher values. Therefore, the credible regions for are not symmetric with respect to the mode of the distribution.

The lowest field strength values correspond to configurations for which the positive magnetic pole is located at the stellar disk center. Larger field values are possible when we are looking at the positive pole from an angle. For example, if the field is aligned (), the amplitude of the Stokes profiles will decrease as the inclination decreases, until we reach where the magnetic signal vanishes.

Because there is only one observation, the aligned configurations are preferred, as shown by the probability density marginalised for . Obviously, when the inclination is , an aligned configuration would not provide any Stokes signal and cannot reproduce the observed profile. Some obliquity would be necessary in order to put the positive pole on the visible hemisphere, which explains the dip at in the probability density marginalised for and marginalised for the - plane.

When and the field is aligned, we are always looking directly at the positive pole, which makes this configuration quite likely. However, our prior knowledge tells us that a low inclination is less likely than a high inclination. For this reason, the posterior probability density marginalised for decreases toward and .

At the present configuration ( and ) the particular profile shape (anti-symmetrical) only occurs during a certain phase range (see Figure 2), which makes such a configuration less likely than the aligned configurations. Therefore, with only one observation, no meaningful constraints can be put on the angles. However, the field strength can be better constrained, and the probability density marginalised for can provide statistical insight into the more likely field strength values.

At 450 G, the odds ratio favours by more than 13 orders of magnitude. When the model is so strongly rejected, the difference in odds ratios between our chosen prior and the flat prior is less important than when the signal detection was more marginal. This is because the bulk of the likelihood is located at higher values, where the two priors are similar. The posterior probability density is quite similar in shape to the one with a 250 G field, with the distributions shifted to higher values, as it can be seen from the credible regions in Table 3.

### 4.3 Realistic case

So far, we have explored the behaviour of the algorithm in cases where the models can reproduce the data perfectly. Obviously, real observations will have some deviations from the predicted values, introduced by the noise. We will now explore how the noise corruption affects our detection capacity.

The top spectrum of Figure 7 shows a simulated observation, consisting of only random noise generated from a normal distribution corresponding to a s/n of 12 000. Figure 8 presents the posterior probability density for that specific observation. Notice how the distribution looks like the probability density for the 100 G observation of the previous section (Figure 6) 111The density probability is different from the 100 G ideal detection because the simulated noise pattern has a shape similar to a magnetic field whose pole is located at a non-null rotational radial velocity. Such a location requires a dipole that is not aligned.. In fact, the odds ratio is even in favour of the dipole model (), even though the signal is pure noise. This is because the observation is better fitted by a non-null magnetic field. The absolute best fit to the data is given by the maximum of the likelihood. We illustrate this fit by the red profile overplotted on the data (Figure 7). Therefore, the detection of a real signal embedded in the noise is ambiguous, as the noise could have – by chance – the shape of a magnetic profile. An observer would likely not feel confident to report a magnetic detection based on only one observation like this one.

In order to verify if this ambiguity can be lifted by re-observing the star, we generated 4 additional simulated observations, again with pure normal noise (rest of Figure 7). Each of these observations leads to an odds ratio in favour of by about a factor of two (Table 4, column 2).

The best fit achievable for each observation taken individually is again overplotted in red. However, nothing restricts the magnetic configuration to be the same for each observation. The best fit produced by a single configuration is given by the maximum of the joint likelihood, illustrated in green. Although the data can be reproduced by a non null magnetic field, the odds ratio of the combined observations is in favour of the model (). The slight improvement of the fit produced by the oblique dipole model (see the reduced indicated in Figure 7), is not enough to justify the use of a more complex model.

The posterior probability density (Figure 9) is similar in shape to the perfect non-detection that was shown in Figure 4. The high- tail is less extended that in the case of a single observation. This can by seen more easily by comparing the 2D probability densities for the - and - planes of Figures 9 and 4. As explained in the previous sections, combining multiple observations has this effect because, in the case of non-detections, the high- inclined dipole configurations become less likely, as this would mean that all the observations were taken during a narrow range of phase. Obviously, the possibility that the observations were taken at the same rotational phase can generally be assessed by the time span of the observations and the possible range of rotational periods. In the case where it is known that the rotational period is quite long compared to the observation time span, more information is available than is assumed in this algorithm (that the phase of each observation can assume any value). In that case, it would be wiser to combine all the observations together and treat them as a single observation in order to get a more meaningful probability density.

Table 5 compiles upper limits to the 68.3%, 95.4%, 99.0% and 99.7% credible regions, extracted from the posterior probability density marginalised for . Although we combined five observations, the 68.3% credible region for extends to 38 G, which is higher than the value obtained for one perfect observation of low s/n (30 G), because of the deviations introduced by the noise. However, as mentioned above, the distribution of the combined observations does not extend as far than as that from a single perfect observation, as illustrated by the 99.7% credible region that extends only to 370 G, compared to 456 G for the perfect observation. We also compiled the credible regions we would obtain without the use of the noise scaling parameter. As expected, the deviations from the model are consistent with our estimation of the variance (the error bars) and the credible regions are nearly the same in both cases.

We have therefore demonstrated that a suspicious signal that has a shape similar to that of a magnetic signal can be verified by the acquisition of a small number of additional observations. We now demonstrate that the same is true for a real signal that is sufficiently embedded in the noise to render the odds ratio of a single observation ambiguous. Given that in the ideal detection case we were able to detect a field with a strength near the 95.4% upper limit of the ideal non-detection, we chose a field strength close to the upper limit of the 95.4% credible region of the previous example. The configuration is given by  G, and . Five phases were randomly generated (), and the corresponding Stokes profiles were added to the previous simulated dataset of pure noise. The resulting simulated observations are shown in Figure 10. The underlying real Stokes profile is shown as the black dashed curve for each observation.

Table 4 gives the odd ratio for each observation (column 4), all of which favour the dipole model , although observations 2, 3 and 5 less strongly than observations 1 and 4. In fact, the odds ratios of the formers are similar to the odds ratio of the first observation of the pure noise test and on their own, each of these observations would result in an ambiguous signal detection. However, combining all the observations together, the odds ratio is now strongly in favour of the dipole model by 9 orders of magnitude ().

Figure 11 shows the posterior probability densities for the combined observations. There is a sharp lower limit on the possible magnetic strengths, illustrated by the probability density marginalised for , and for the 2D - and - planes. Given that most of the likelihood is situated at non-null field strengths, the improvement of the fit to the data justifies the more complex model, as shown by the reduced on Figure 10 for the model fits (black lines) and the fit produced by maximum of the joint likelihood for (green curves).

Different choices are possible in order to express the derived values for the model parameters. One can choose for example the maximum of the joint posterior probability density (MAP), or the mode of the marginalised probability density for each parameter. Note however that if the probability distribution is complex, the parameters given by the MAP do not necessarily correspond to the mode of the marginalised probability densities. Usually, the MAP will produce the best fit to the data given that the a-priori information does not exclude any interesting parts of the parameter space, but does not necessarily represent the bulk of the probability. The mode of each parameter represents well the bulk of the probability, but does not necessarily give an excellent fit to the data. Using the median of the marginalised probability densities is usually a good compromise (Gregory, 2005b).

In Table 6, we list the MAP, the mode and the median for each parameter, which in this case are quite similar. The MAP, the modes and the medians all yield similarly good fits to the data (the MAP and the modes fits are shown in Figure 10 in blue and magenta respectively).

The range of the credible regions of the probability density marginalised for each parameter are compiled in Table 6. Using the median with the 68.3% credible region, we would infer a dipole strength  G, which is slightly higher than the input dipole field strength of 125 G. We verified that this difference was due to the particular noise pattern used in the data, as additional noise simulations allowed us to recover the input value within the 68.3% region. The real values of the and angles are recovered by the 68.3% credible regions, although the constraints are poor, as expected.

### 4.4 Comparison with traditional diagnostics

With low-resolution instruments (such as FORS1 and FORS2 at the Very Large Telescope), one is generally only sensitive to the global longitudinal component of the magnetic field, as the rotationally-broadened spectral lines are not resolved. This global longitudinal field value is extracted from the spectrum using the “slope” method as described in Bagnulo et al. (2002, 2006). The presence of a magnetic field in a single observation is therefore diagnosed by the significance of the global longitudinal field measurement compared to its error bar.

In the case of high-resolution spectropolarimetry (e.g. ESPaDOnS at Canada-France-Hawaii Telescope, Narval at Télescope Bernard-Lyot or HARPSPol at ESO La Silla 3.6 m Telescope), the rotationally-broadened spectral lines are resolved and the field is generally diagnosed by the deviation of the circular polarisation profile with respect to . The detection can be quantified by the probability that such a deviation from is produced by random noise (the false alarm probability or FAP). The detection probability can then be expressed as . Following Donati et al. (1997), a field is generally considered detected if the FAP is less than (), and marginally detected when FAP is less than (). With high-resolution spectropolarimetry, it is also possible to integrate the signal over the line profile in order to recover a value equivalent to the global longitudinal field obtained from low-resolution instruments (Donati et al., 1997; Wade et al., 2000a). Although it is possible to detect a magnetic signal even when the global longitudinal field is null, this is a useful quantity for producing longitudinal field curves and for comparing with low-resolution data.

We applied these traditional diagnostics to our simulated datasets and we report the results in Table 7. In the case of the ideal detections, both the longitudinal field (columns 2 and 3) and the detection probability (column 4) yield a detection when the field strength reaches 450 G. However, the odds ratios are already in favour of the magnetic model at 250 G. The detection probability only looks at the total deviation, whereas the odds ratio looks for a shape similar to that of the theoretical model.

For the realistic case consisting only of noise, all five observations have longitudinal field measurements below significance (see column 3). Moreover, no signal is detected in the Stokes profiles, as shown by the detection probabilities. The same is also true for the realistic case consisting of noise plus a signal for  G, although the detection probability nearly reaches a marginal detection for Obs. 1 and 4. For these observations, the odds ratios were in favour of the magnetic model by 2 and 4 orders of magnitude respectively. For the remaining observations, the odds ratios were also in favour of the magnetic model, although at a lower significance. Therefore, if we had observed any of the simulated observations with , the traditional diagnostics would not have diagnosed the presence of a field, but the odds ratio would have indicated the possible magnetic signal. A few additional observations are enough to distinguish between real noise and a buried signal consistent with an oblique dipole field, as demonstrated by the example here. Therefore, Bayesian odds ratios provide a quantitative indication of stars worth re-observing in magnetic surveys.

## 5 Application to real data

In order to present the application of this method to real data, we will use high-resolution spectropolarimetric observations of the magnetic B-type star LP Ori (=Par 1772, HD 36982).

LP Ori is often considered a chemically peculiar star because it was first classified as a B1.5p star by Sharpless (1952). LP Ori’s status as a He-strong or He-weak star is still uncertain. An inspection of our seven spectra (described below) does not reveal any significant variation in the He-line strength that would indicate a He-strong or He-weak star. Furthermore, comparing our spectra with the BSTAR grid of synthetic spectra from non-LTE tlusty models (Lanz & Hubeny, 2007), we find a reasonable agreement with a temperature of 20 kK, a of 4.0 and a of 80 km s. LP Ori is also a candidate Herbig Ae/Be star, because of its far-infrared excess, although no emission is present in the visible spectrum. Manoj et al. (2002) suggested that LP Ori is an object transiting from the pre-main sequence to the main sequence. LP Ori seems to be a single star. No radial velocity variation was found by Abt et al. (1991) nor was any speckle companion with a K-band ratio of less than 0.04 for a separation of 150 mas or more (Preibisch et al., 1999).

The magnetic field of LP Ori was first reported by Petit et al. (2008). In that paper, they used three Stokes observations obtained with the ESPaDOnS and Narval spectropolarimeters to detect the magnetic field.

ESPaDOnS and Narval are twin high-resolution spectropolarimetric instruments located at Canada-France-Hawaii Telescope and Télescope Bernard-Lyot respectively. A polarisation measurement consists of a set of 4 sub-exposures taken with different polarimeter configurations. From this measurement set, the circular polarisation Stokes spectrum is extracted, as well as a diagnostic null polarisation spectrum (labeled ) by combining the sub-exposures in such a way that the astronomical object’s polarisation should cancel out. ESPaDOnS frames were processed using the Upena pipeline provided by CFHT. The Narval frames were processed by the TBL archive pipeline. Both pipelines use the reduction package Libre-ESPRIT (Donati et al., 1997). The spectral range of both instruments covers the 370 nm to 1050 nm wavelength band, with a resolution .

Petit et al. (2008) used the Bayesian method described here to estimate the dipole strength of LP Ori, obtaining  G. For the present analysis, we have obtained one additional ESPaDOnS observation (within the context of the MiMeS CFHT Large Program) and an additional 3 Narval observations. The observation log, which includes the new observations as well as those analysed by Petit et al. (2008), is given in Table 8222During the analysis of the complete set of observations of LP Ori for this paper, it was discovered that the sign of the Stokes profiles corresponding to two spectra employed by Petit et al. (2008) was inverted (the Jan. 2006 ESPaDOnS spectrum and the Narval spectrum). During the present analysis we have corrected the sign of these spectra and verified the sign of all others. We have verified that the results of Petit et al. (2008) are not substantially modified by this change.. The sample of observations has a large range of s/n, and is therefore suitable for testing the behaviour of the Bayesian method on real data.

As is customary, we applied the LSD procedure to our observations. We used the iLSD code described by Kochukhov et al. (2010). We chose spectral lines from a Vienna Atomic Line Database (VALD; Kupka et al., 2000) list corresponding to an atmosphere model with  kK and . From that list, we chose metallic lines and weak He lines that were unblended with strong Balmer lines and uncontaminated by telluric lines or strong nebular emission. The line depths have been slightly adjusted to match the spectra. The final line list is shown in Table 11. The intensity () and polarisation () weights of each line were normalised by 0.2 and 120 respectively. Therefore, the displayed -axis scale of the resulting Stokes and diagnostic null LSD profiles (Figures 12 and 13) corresponds to a line with a and .

The two sharp lines in the blue continuum of the Stokes LSD profiles are residuals from telluric lines adjacent to the HeI line. The emission bump present in the LSD profiles of the two first Narval observations are residuals from He emission, most likely of nebular origin because their centroid velocities and scaling match the nebular Balmer line emission.

We computed the traditional diagnostics – the detection probability of the individual profiles and the global longitudinal field – by integrating  km s around the line centre, for both the Stokes and profiles. The results are displayed in Table 8. Definite detections (column 8) are achieved for the three ESPaDOnS observations. Marginal detection is achieved only for one Narval observation. No signal is detected in the null profiles (column 11).

The global longitudinal fields for the Stokes LSD profiles (column 6) show a positive trend, hinting that the positive magnetic pole is located somewhere on the visible hemisphere. Although the longitudinal field approaches zero, the magnetic signal is still detectable given a sufficient s/n (for example 2007-03-06). This hints that we are looking nearly at the magnetic equator at some phases.

Assuming a reasonable range of possible rotation periods, defined by the and breakup velocity, the longitudinal field curve and the LSD profiles can be phased with various periods. Figure 14 shows the longitudinal field measurements (black dots), along with sinusoidal curves for three possible periods. As no other variability has been observed by other means (spectral or photometric), this is a good example of a case where the rotational phases of the observations are not known.

Given the range of global longitudinal field measurements, we extended the grid up to  kG. We set the parameter of the Jeffreys prior to 25 G, corresponding to two times the parameter grid step. We performed the analysis on both the Stokes and null profiles. The fit of the Stokes profiles is shown in Figure 12. Vertical dotted lines represent the velocity range of the fits.

The odds ratios for each individual observation are given in Table 9. The odds ratios for the null profiles (column 3) all favour the absence of a magnetic field ( model). The odds ratios for the Stokes observations all favour the oblique dipole model (), but vary from more than 10 orders of magnitude for the high s/n observations to less than one order of magnitude for the low s/n observations. As mentioned for realistic numerical tests (Section 4.3), it is possible to obtain an odds ratio in favour of the dipole model when the noise happens to have a shape similar to a magnetic signal. When combining all the observations together, we get an odds ratio strongly in favour of the magnetic model, by 73 orders of magnitude, as expected given the definite signal detections in the ESPaDOnS observations. However, if only the low s/n observations were available (2007-11-08, -09 and -10), the situation would be more ambiguous. None of these observations lead to a traditional definite detection when considered on their own. We therefore combined and analysed these three observations in pairs, and then all together. The three possible pairs of observations led to odds ratios in favour of by more than one order of magnitude, and the combination of the three observations lead to an odds ratio , i.e. more than 3 orders of magnitudes. Therefore, the Bayesian algorithm is able to recover the magnetic signal that would have been buried in the noise and undetected by traditional diagnostics.

In Figures 12 and 13, the grey dashed lines represent the model (). The best fits achievable by a dipole model for each observation taken individually (maximum of the individual likelihoods) are shown in red. The best fit produced by a single oblique dipole (maximum of the joint likelihood) is shown in green. The reduced s are indicated with corresponding colours. For the Stokes observations where odds ratios strongly favour the model, the reduced is much improved with the addition of the more complex model compared to the of the model. Not only are the data reproduced by the dipolar profiles, they can all be simultaneously reproduced by a single configuration, as shown by the similar for the individual fit (red) and the maximum of the joint likelihood (green).

However, the reduced remains high in one case (2010-02-23; =1.95), meaning that there is some extra variance in the observation that neither models are able to reproduce (a 3 deviation would correspond to a reduced of 1.65). The reduced s are low for all the null profile observations (Figure 13). This points toward a systematic deviation or a model-based effect for the Stokes observation of Feb. 2010 rather than an extra instrumental scatter or underestimated error bars.

When performing parameter estimation, the extra scatter is addressed by the noise scaling term (Eq. 17). It is therefore possible to extract the probability density function marginalised for the noise scaling parameter, and determine if the model is able to reproduce the observations down to the noise level. Figure 15 shows that the Bayesian estimate of the variance is larger than the assumed variance (i.e. error bars) for both Stokes and the null profiles, but by less than a factor of two.

The posterior probability density for the Stokes observations is shown in Figure 16. Dipole configurations with large inclination or obliquity are less favoured, as shown by the probability densities marginalised for and , because the observations mainly show either the positive pole or the magnetic equator on the visible hemisphere. This is more likely to occur if the and angles are small and the negative pole spends little or no time on the visible hemisphere. The angle values are interrelated as shown by the - 2D plane. If the inclination is small, the obliquity is more likely to be large, in order to display an equator-like magnetic signature. The probability density marginalised for the field strength shows a sharp lower limit. The high- tail of the distribution is attributable to the high-inclination (low obliquity) configurations, as shown in the 2D planes for - and -.

Figure 17 shows the probability densities for the null profile observations, which are, as expected, similar to the realistic simulation of pure noise (as shown in Figure 9). We also show in Figure 18 the probability densities we obtained when considering only the three low s/n observations. The constraints on the parameters, especially the angles, are worse than when we consider the full data set. Therefore, the constraints on the angles are mainly defined by the high s/n observations.

The credible regions of the probability density marginalised for are given in Table 10. We also present the maximum of the joint posterior probability density (MAP), the mode of the probability density marginalised for , as well as the median. For the Stokes observations, the MAP, mode and median values are all similar, therefore the fits to the data obtained for these values are nearly undistinguishable. The fit obtained with the MAP values is shown in Figure 12 (blue curve). Using the 68.3% credible region and the median, the dipole field strength of LP Ori is estimated to be  G. As was seen for the probability density in Figure 16, the lower limit on the field strength is sharp. The 68.3% credible region lower limit is 667 G, and 557 G for the 99.7% credible region. The distribution has a high- tail, and the 99.7% credible region extends up to 2.6 kG.

The second row of Table 10 illustrates the effect of the noise scaling parameter on the inferred field values. When considering the variance as a model parameter, the MAP, mode and median are shifted to slightly lower values. The credible regions are somewhat larger, and also slightly shifted to lower values.

The parameter estimation for the low s/n observations is less robust than that of the full dataset. The estimated field value from the median and the 68.3% credible region is  G. However, we do not have a good constraint on the lower limit, as the credible regions quickly go to zero. The 99.7% credible region extends to 2.7 kG as well.

For the null profile observations, the MAPs and modes are all 0 G. Given the shape of the probability distribution, the median is non-null. However, given the probability density in favour of the model, it makes more sense to express the field strength estimation in terms of the upper limit of the credible regions. The null profiles are a good representation of what our data would look like in the absence of a stellar magnetic field. For the whole dataset, the upper limit of the 95.4% credible region is 144 G, well below the inferred dipole strength ( kG) from the Stokes observations. For the 99.0% credible region, we can say that the probability that an undetected field would have a dipole strength of more than 300 G is only 1%. The noise scaling does not significantly change the credible regions. Interestingly, the dipole strength inferred from the low s/n Stokes observations is around 600 G, and the odds ratio favours the magnetic model. From the low s/n null profile observation, we can see that a field with a dipole strength of the order of the 95.4% credible region upper limit (502 G) can indeed be detected.

## 6 Conclusion

In this paper, we have described a method based on Bayesian statistics to infer the magnetic properties of stars observed spectropolarimetrically in the context of large surveys like the Magnetism in Massive Stars project. This approach is well-suited for stars for which the stellar rotation period, and therefore the rotational phases of the small number of observations, are not known.

The model used to predict the expected Stokes profiles is that of an oblique dipolar magnetic field, parametrised by the field strength at the pole, the inclination of the rotational axis, the obliquity of the magnetic axis with respect to the rotational axis and the rotational phase (which is allowed to take any value for each individual observation). In the present case, the calculations are performed under the weak-field approximation, although any polarised spectral synthesis code can in principle be used with the Bayesian algorithm.

The result of the analysis is a multidimensional posterior probability density that describes the relative likelihood of models spanning the parameter space of the dipolar field model. We have used synthetic observations to explore the behaviour of the Bayesian algorithm under ideal and realistic conditions. In the case of an ideal non-detection, the posterior probability density for the field strength has the form of a decreasing exponential, the extended tail of the field strength distribution being due to a specific family of dipole orientations. This tail becomes less extended with the addition of multiple observations, as some of these specific dipole orientations only present low-amplitude or null Stokes profiles over a restricted range of phases. However, the possibility that the dipole is aligned with the rotational axis and seen equator-on ( and ) always remains as this configuration never produces any circular polarisation at any phase. When a detectable field is present, the probability distribution shows a sharp lower limit on the dipolar field strength, and a slow decrease towards higher field strengths, producing an asymmetrical distribution. With only one observation, not much can be inferred about the dipole orientation. The longitudinal field indicates which pole is located on the visible hemisphere and the shape of the Stokes V profile indicates when some obliquity is necessary because the pole is located at a non-null rotational radial velocity.

A particularly useful quantity that can be computed from the posterior probabilities is the so-called “odds ratio”, which compares the relative compatibility between the observations and the magnetic dipole hypothesis versus the non-magnetic hypothesis. By adding a magnetic signal corresponding to the upper limit of the credible regions found for the ideal non-detection, we have explored the detection capability of the odds ratios. We find that fields corresponding to the 68.3% and 95.4% region are below detection, whereas those corresponding to the 99.0% and 99.7% regions are well detected with the odds ratios. In contrast, traditional diagnostics (detection probability and global longitudinal field significance) only detect fields corresponding to the 99.7% region.

By using a set of five realistic simulated observations, we have also shown that in the case of noise emulating the shape of a weak magnetic signal, it is possible to use the odds ratios to distinguish between noise and real signal by obtaining a small number of additional observations. Combining all the observations together, it is therefore possible to detect a weak magnetic field (with a strength corresponding to roughly the upper limit of the 95.4% credible region of the noise-only case) under the detection capabilities of the traditional diagnostics. We have therefore shown that the odds ratio is an powerful quantitative indicator of which undetected stars in a survey should be re-observed in priority.

In most applications where a field is indeed detected, the resulting probability density will generally be marginalised for the dipole strength, as only limited information can be obtained for the rotational inclination and the magnetic obliquity from a few observations (the magnetic geometry is recovered by the most probable inclination and obliquity, but the credible regions are quite extended). We have shown that the dipole strength probability distribution provides a reasonable estimate of the field strength.

We have applied our method to real spectropolarimetric observations of the magnetic B-type star LP Ori. The dataset consists of 3 ESPaDOnS and 4 Narval observations, of various signal-to-noise ratios. A magnetic signal is indeed detected by the odds ratios, even when only considering the low s/n observations where the traditional diagnostics do not detect the magnetic field. Using all the available spectra, we used the median of the marginalised posterior probability density, as well as the 68.3% credible region, to infer a dipolar field strength of  G. Although the probability density for the obliquity and inclination do not provide any tight constraints, geometries for which the negative magnetic pole spends little or no time on the visible hemisphere are preferred, since the dataset consists of positive or nearly null longitudinal field measurements. We also performed our analysis on the diagnostic null spectra, which resulted in a non-detection. The null profiles provide an useful verification of spurious signals and also provide an estimate of the field strengths that would be detectable with data of such quality. For example, when considering only the low s/n observations, the Stokes profile analysis yielded the detection of the  G dipolar field and the diagnostic null profile analysis yielded an upper limit of  G for the 95.4% credible region.

## Acknowledgments

VP acknowledges support from Fonds québécois de la recherche sur la nature et les technologies. GAW acknowledges support from the Discovery Grants programme of the Natural Science and Engineering Research Council of Canada. We thank M. Gagné and R. P. Breton for useful discussions, and the anonymous referee whose helpful comments led to the improvement of this paper.

## References

• Abt et al. (1991) Abt H. A., Wang R., Cardona O., 1991, ApJ, 367, 155
• Arregui & Asensio Ramos (2011) Arregui I., Asensio Ramos A., 2011, ApJ, 740, 44
• Asensio Ramos (2011) Asensio Ramos A., 2011, ApJ, 731, 27
• Aurière et al. (2007) Aurière M., Wade G. A., Silvester J., Lignières F., Bagnulo S., Bale K., Dintrans B., Donati J. F., et al., 2007, A&A, 475, 1053
• Babel & Montmerle (1997a) Babel J., Montmerle T., 1997a, ApJ, 485, L29
• Babel & Montmerle (1997b) Babel J., Montmerle T., 1997b, A&A, 323, 121
• Bagnulo et al. (2006) Bagnulo S., Landstreet J. D., Mason E., Andretta V., Silaj J., Wade G. A., 2006, A&A, 450, 777
• Bagnulo et al. (2002) Bagnulo S., Szeifert T., Wade G. A., Landstreet J. D., Mathys G., 2002, A&A, 389, 191
• Borra et al. (1981) Borra E. F., Fletcher J. M., Poeckert R., 1981, ApJ, 247, 569
• Borra & Landstreet (1980) Borra E. F., Landstreet J. D., 1980, ApJS, 42, 421
• Donati et al. (1997) Donati J.-F., Semel M., Carter B. D., Rees D. E., Collier Cameron A., 1997, MNRAS, 291, 658
• Gray (1992) Gray D. F., 1992, The Observation and Analysis of Stellar Photospheres.. Camb. Astrophys. Ser., Vol. 20
• Gregory (2005a) Gregory P. C., 2005a, ApJ, 631, 1198
• Gregory (2005b) Gregory P. C., 2005b, Bayesian Logical Data Analysis for the Physical Sciences: A Comparative Approach with ‘Mathematica’ Support. Cambridge University Press, Cambridge, UK.
• Hunter et al. (2008) Hunter I., Lennon D. J., Dufton P. L., Trundle C., Simón-Díaz S., Smartt S. J., Ryans R. S. I., Evans C. J., 2008, A&A, 479, 541
• Jaynes (2003) Jaynes E. T., 2003, Probability Theory: The Logic of Science. Cambridge University Press, Cambridge, UK
• Jeffreys (1998) Jeffreys H., 1998, Theory of Probability; 3 edition. Oxford University Press, Oxford, UK
• Kochukhov et al. (2010) Kochukhov O., Makaganiuk V., Piskunov N., 2010, A&A, 524, A5
• Kochukhov & Wade (2010) Kochukhov O., Wade G. A., 2010, A&A, 513, A13
• Kupka et al. (2000) Kupka F. G., Ryabchikova T. A., Piskunov N. E., Stempels H. C., Weiss W. W., 2000, Baltic Astronomy, 9, 590
• Landi degl’Innocenti & Landolfi (2004) Landi degl’Innocenti E., Landolfi M., 2004, Polarization in Spectral Lines. Vol. 307 of Astrophysics and Space Science Library, Kluwer Academic Publishers, Dordrecht, The Netherlands
• Landstreet & Borra (1978) Landstreet J. D., Borra E. F., 1978, ApJ, 224, L5
• Landstreet & Mathys (2000) Landstreet J. D., Mathys G., 2000, A&A, 359, 213
• Lanz & Hubeny (2007) Lanz T., Hubeny I., 2007, ApJS, 169, 83
• Maeder & Meynet (2005) Maeder A., Meynet G., 2005, A&A, 440, 1041
• Manoj et al. (2002) Manoj P., Maheswar G., Bhatt H. C., 2002, MNRAS, 334, 419
• Petit et al. (2008) Petit V., Wade G. A., Drissen L., Montmerle T., Alecian E., 2008, MNRAS, 387, L23
• Preibisch et al. (1999) Preibisch T., Balega Y., Hofmann K.-H., Weigelt G., Zinnecker H., 1999, New Astronomy, 4, 531
• Semel & Li (1996) Semel M., Li J., 1996, Sol. Phys., 164, 417
• Sharpless (1952) Sharpless S., 1952, ApJ, 116, 251
• Shore & Brown (1990) Shore S. N., Brown D. N., 1990, ApJ, 365, 665
• Townsend et al. (2010) Townsend R. H. D., Oksala M. E., Cohen D. H., Owocki S. P., ud-Doula A., 2010, ApJ, 714, L318
• Tuomi et al. (2011) Tuomi M., Pinfield D., Jones H. R. A., 2011, A&A, 532, A116
• ud-Doula et al. (2009) ud-Doula A., Owocki S. P., Townsend R. H. D., 2009, MNRAS, 392, 1022
• Wade et al. (2011) Wade G. A., Alecian E., Bohlender D. A., Bouret J., Cohen D., Duez V., Gagné M., et al., 2011, in Active OB stars Vol. 272 of IAU Symposium. p. 118
• Wade et al. (2001) Wade G. A., Bagnulo S., Kochukhov O., Landstreet J. D., Piskunov N., Stift M. J., 2001, A&A, 374, 265
• Wade et al. (2000a) Wade G. A., Donati J.-F., Landstreet J. D., Shorlin S. L. S., 2000a, MNRAS, 313, 851
• Wade et al. (2000b) Wade G. A., Donati J.-F., Landstreet J. D., Shorlin S. L. S., 2000b, MNRAS, 313, 823
• Wade et al. (2011) Wade G. A., Howarth I. D., Townsend R. H. D., Grunhut J. H., Shultz M., Bouret J.-C., Fullerton A., Marcolino W., Martins F., Nazé Y., ud-Doula A., Walborn N. R., Donati J.-F., the MiMeS Collaboration 2011, MNRAS, 416, 3160
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters