Analytical maximum likelihood estimation of stellar magnetic fields
M. J. Martínez González, R. Manso Sainz, A. Asensio Ramos, L. Belluzzi
aa r X i v : . [ a s t r o - ph . S R ] A ug Mon. Not. R. Astron. Soc. , 1–12 (2002) Printed 13 April 2018 (MN L A TEX style file v2.2)
Analytical maximum likelihood estimation of stellarmagnetic fields
M. J. Mart´ınez Gonz´alez , , R. Manso Sainz , , A. Asensio Ramos , , and L. Belluzzi , Instituto de Astrof´ısica de Canarias, V´ıa L´actea s/n, 38205, La Laguna, Tenerife, Spain Departamento de Astrof´ısica, Universidad de La Laguna, 38205, La Laguna, Tenerife
ABSTRACT
The polarised spectrum of stellar radiation encodes valuable information on the condi-tions of stellar atmospheres and the magnetic fields that permeate them. In this paper,we give explicit expressions to estimate the magnetic field vector and its associatederror from the observed Stokes parameters. We study the solar case where specificintensities are observed and then the stellar case, where we receive the polarised flux.In this second case, we concentrate on the explicit expression for the case of a slowrotator with a dipolar magnetic field geometry. Moreover, we also give explicit formu-lae to retrieve the magnetic field vector from the LSD profiles without assuming meanvalues for the LSD artificial spectral line. The formulae have been obtained assumingthat the spectral lines can be described in the weak field regime and using a maximumlikelihood approach. The errors are recovered by means of the hermitian matrix. Thebias of the estimators are analysed in depth.
Key words: techniques: polarimetric – stars: magnetic field
The most precise measurements of stellar magnetic fieldsare based on the observation and the interpretation ofpolarisation in spectral lines. Most of the magnetic fieldsdetected in stars at different phases of evolution havebeen found by means of the Zeeman effect, which gener-ates linear and circular polarisation in the presence of amagnetic field (e.g., Wade et al. 2000; Bagnulo et al. 2002;Jordan, Werner & O’Toole 2005; Aznar Cuadrado et al.2004; O’Toole et al. 2005; Silvester et al. 2009; Leone et al.2011). Here we present analytical expressions for infer-ring the magnetic field vector from the observed Stokesprofiles induced by the Zeeman effect in the weak fieldapproximation.The weak field approximation is broadly applied for theinference of solar and stellar magnetic fields from the obser-vation of the Stokes profiles. It is an analytical solution tothe radiative transfer equation whose main basic assump-tion is that the magnetic field is sufficiently weak throughoutthe whole region of the atmosphere where a spectral line isformed. Although simple, this approximation is very usefulsince non-magnetic mechanisms usually dominate the shapeof spectral lines, in both the solar and the stellar case.For example, in the quiet Sun, which occupies thevast majority of the solar surface (far from the sunspotsthat harbour very strong fields), the spectral lines are welldescribed by the weak field approximation. This has allowedthe success of many synoptic magnetographs like those of Big Bear (Spirock et al. 2001; Varsik 1995). It is evenused to produce modern vector magnetograms like thoseobtained with the IMaX instrument (Mart´ınez Pillet et al.2011) onboard the Sunrise balloon (Solanki et al. 2010).The weak field approximation has also been usedto diagnose the chromosphere (Merenda et al. 2006;Asensio Ramos, Trujillo Bueno & Landi Degl’Innocenti2008) due to the enhanced thermal width of the spectrallines.In night-time spectro-polarimetry, the weak field ap-proximation is at the base of the least-squares deconvolu-tion (LSD; Donati et al. 1997), the most successful tech-nique used to detect and measure magnetic fields in solar-type stars or other stars in which the polarisation signalper spectral line is well below the noise level. Other workshave used this approximation to diagnose magnetic fieldsin a large variety of stellar objects. A limited selection in-cludes some recent works on central stars of planetary neb-ulae (Jordan, Werner & O’Toole 2005; Leone et al. 2011),white dwarfs (Aznar Cuadrado et al. 2004), pulsating stars(Silvester et al. 2009), hot subdwarfs (O’Toole et al. 2005),Ap and Bp stars (Wade et al. 2000) and chemically peculiarstars (Bagnulo et al. 2002).
The weak field approximation is an analytical solution to theradiative transfer equation. The fundamental assumption is c (cid:13) M. J. Mart´ınez Gonz´alez that the magnetic field vector is constant and its intensityis sufficiently weak throughout the whole region of the at-mosphere where an spectral line is formed. Additionally, theline-of-sight velocity and any broadening mechanism haveto be constant with height in the line formation region. Asa consequence, the magnetic field can be considered as aperturbation to the zero-field case. Quantitatively, the ap-proximation holds whenever ∆ λ B ≪ ∆ λ D , where ∆ λ B isthe Zeeman splitting and ∆ λ D is the dominant broadeningmechanism (thermal, rotation, etc.). From its definition, itis clear that the weak field regime occurs at different fieldstrengths for different spectral lines (depending on the sen-sitivity to the magnetic field, the local temperature and theatomic mass) and for different stellar objects (dependingon any non-magnetic broadening mechanism for the spec-tral lines, the rotation being the most efficient mechanismin cool stars).To first order in ∆ λ B , the intensity profile of a spectralline formed in a weak magnetic field is insensitive to themagnetic field. In other words, it fulfils the transfer equationin the absence of a magnetic field. At this first order, thecircular polarisation profile, i. e., the Stokes V profile, for agiven spectral line has the following expression: V ( x ) = − C Λ gB || ∂I ( x ) ∂x . (1)The symbol B || stands for the longitudinal component ofthe magnetic field, i. e., B || = B cos θ B , where θ B is theinclination of the magnetic field with respect to the ob-server’s line of sight ( ~ Ω; see Fig. 1) and B is the mag-netic field intensity. The symbol g represents the effectiveLand´e factor of the line, which quantifies the magneticsensitivity of the line. The effective Land´e factor only de-pends on the quantum numbers of the transition (see, e.g.,Landi Degl’Innocenti & Landolfi 2004). Equation 1 is ex-pressed in terms of the generic wavelength variable x . If x represents the wavelength λ , then the parameter Λ equalsthe central wavelength of the line, λ . However, it is cus-tomary in stellar spectropolarimetry to consider x as thevelocity in Doppler units. In such a case, we have Λ = cλ ,with c the speed of light. The constant C = 4 . × − G − ˚A − .At first order in ∆ λ B , both Stokes Q and U are zero. Inorder to obtain an expression for the Stokes profiles charac-terizing linear polarisation we have to expand the radiativetransfer equation to second order in ∆ λ B and assume thatthe spectral line is not saturated. Under these assumptions,the general formulae for the Stokes Q and U are: Q ( x ) = − C GB ⊥ cos 2 φ B ∂ I ( x ) ∂x ,U ( x ) = − C GB ⊥ sin 2 φ B ∂ I ( x ) ∂x . (2)The symbol G plays the role of the effective Land´e factorfor linear polarisation and quantifies the sensitivity of linearpolarisation to the magnetic field. Again, it is only a func-tion of the quantum numbers of the transition (see, e.g.,Landi Degl’Innocenti & Landolfi 2004). The symbol B ⊥ = B sin θ B is the component of the magnetic field perpendicu-lar to the line of sight. The angle φ B represents the azimuthof the magnetic field vector with respect to an arbitrary ref-erence direction ( ~e a , ~e b in Fig. 1 show the coordinates chosenfor the reference of Q > R ∼ T ∼ −
20 km s − . However, when observingmetal lines at spectral resolutions as high as R ∼ I ( x, µ ) = I ( x ) (cid:0) − u − v + uµ + vµ (cid:1) = I ( x ) f ( µ ) , (3)where I ( x ) is the intensity profile at disc centre ( µ = 1).The parameters u and v have values between 0 and 1 andare supposed to be constant along the spectral line (althoughthey can vary from line to line). Note that the values of u and v have to fulfil the condition that I ( x ) >
0. The CLVis given in terms of µ = cos Θ, where Θ is the astrocentricangle between the normal to a point in the stellar surface andthe line of sight. Using this law for the CLV (and neglectingrotation) the flux for Stokes I is: F I ( x ) = I ( x ) Z d Σ ′ f ( µ ) , (4)where the integral is computed over the visible surface inthe plane of the sky Σ ′ . For an arbitrary function h ( ρ, α ) ofthe polar coordinates of the stellar surface ρ and α , we have Z Z d Σ ′ h ( ρ, α ) = Z ρdρ Z π dαh ( ρ, α ) . (5)Note that the variable ρ = sin Θ, which means that µ = p − ρ . Plugging the quadratic expression considered forthe CLV, the final closed expression for the integrated Stokes I is found to be: F I ( x ) = πI ( x ) (cid:18) − u − v (cid:19) . (6) c (cid:13) , 1–12 aximum likelihood estimation of stellar magnetic fields θ Β φ Β e a e b Ω P ’ P ’ Σ ’ P B Θ ρ BP Σ αρα Figure 1.
Geometry of the stellar model. The symbol ~ Ω represents the line of sight. The vector ~B displays the magnetic field vector inthe stellar surface. Its geometry is described in terms of the inclination with respect to the line of sight θ B and the azimuthal angle φ B (the axis ~e a being the 0 reference for the azimuth). The position of a point P in the stellar surface (Σ) is defined only by the astrocentricangle Θ. Its projection (P’) on the plane of the sky (Σ ′ ) is represented in cylindrical coordinates by the modulus ρ and the angle α ,which is referring to the ~e a and ~e b axis. Integrating Eqs. (1) and (2) over the visible surface, we endup with the following expressions for the polarised flux: F V ( x ) = − C Λ g h B k i ∂ F I ( x ) ∂x , F Q ( x ) = − C G h B ⊥ cos 2 φ B i ∂ F I ( x ) ∂x , F U ( x ) = − C G h B ⊥ sin 2 φ B i ∂ F I ( x ) ∂x , (7)with the definitions: h B k i = 6 π (6 − u − v ) Z d Σ ′ B k f ( µ ) h B ⊥ cos 2 φ B i = 6 π (6 − u − v ) Z d Σ ′ B ⊥ cos 2 φ B f ( µ ) h B ⊥ sin 2 φ B i = 6 π (6 − u − v ) Z d Σ ′ B ⊥ sin 2 φ B f ( µ ) . (8)Therefore, the weak-field approximation for the integratedpolarised flux remains formally the same but the compo-nents of the magnetic field B k and B ⊥ appear weightedby the CLV law. In general, the components of the mag-netic field may depend on the position on the stellar sur-face in a complicated manner, so that the previous inte-grals might not have closed expressions. In the simple casethat the magnetic field is constant across the stellar surface,we recover h B k i = B k , h B ⊥ cos 2 φ B i = B ⊥ cos 2 φ B , and h B ⊥ sin 2 φ B i = B ⊥ sin 2 φ B . One of the simplest non-trivialconfiguration we can consider explicitly is that of a dipolarfield. The magnetic field vector at each surface point is then: B ( r ) = − H d e − e · r ) r ] , (9)where the unitary vector e defines the orientation of thedipole and the unit vector r indicates positions on the stel-lar surface. The quantity H d represents the magnetic fieldstrength at the poles of the dipolar field. From Eq. 7 andEq. 9 and after some algebra, it is possible to write the flux of the Stokes parameters as: F V ( x ) = − C Λ g
15 + u − u − v H k ∂ F I ( x ) ∂x F Q ( x ) = − C G − u − v − u − v H ⊥ cos 2 φ d ∂ F I ( x ) ∂x F U ( x ) = − C G − u − v − u − v H ⊥ sin 2 φ d ∂ F I ( x ) ∂x , (10)where H k = H d cos θ d and H ⊥ = H d sin θ d , θ d being theinclination of the axis of the dipole while φ d is the azimuth ofthe dipole (for a similar derivation see Landolfi et al. 1993).The previous formalism has demonstrated that theweak-field approximation leads to formally the same ex-pressions either if we consider specific intensities (that isapplied whenever spatial resolution is available) or if weconsider fluxes (whenever the object of interest cannot beresolved). In other words, the observed circular (linear) po-larised spectrum is proportional to the first (second) deriva-tive of the observed intensity through the longitudinal (or-thogonal) component of the magnetic field. The only differ-ence resides on the exact definition of the components of thefield which, in the spatially unresolved case, are an averageover the stellar surface weighted by the CLV. Therefore, forthe sake of simplicity, it is possible to combine all expressionson the following general ones: V = − C B k I ′ Q = − C B ⊥ cos 2 φ I ′′ U = − C B ⊥ sin 2 φ I ′′ . (11)The meaning of the newly defined variables is summarizedin Tab. 1. Once the model is set, our aim is to infer the magneticfield vector parametrised in terms of B k , B ⊥ , and φ from c (cid:13) , 1–12 M. J. Mart´ınez Gonz´alez
Table 1.
Variables in Eq. 11.Variable Resolved Dipole I I F I V V F V Q Q F Q U U F U I ′ Λ g ∂I ( x ) ∂x
110 15+ u − u − v Λ g ∂ F I ( x ) ∂x I ′′ Λ G ∂ I ( x ) ∂x − u − v − u − v Λ G ∂ F I ( x ) ∂x B k B k H k B ⊥ B ⊥ H ⊥ φ φ B φ d θ θ B θ d B B H d observations of a set of spectral lines ( I i , Q i , U i , V i ) with i = 1 . . . n lines . Assuming that the weak field approximationcan be applied for the set of observed spectral lines andthat observations are corrupted with uncorrelated Gaussiannoise, we can use a least-squares estimator (maximum likeli-hood) to retrieve the magnetic field vector (see the Appendixfor more details). The χ merit function is defined as: χ = X ij ( V ij − V i ; modj ) ( σ i V j ) + X ij ( Q ij − Q i ; modj ) ( σ i Q j ) + X ij ( U ij − U i ; modj ) ( σ i U j ) , (12)where the subindex j indicates the position along the coordi-nate x . The label mod refer to the model for Stokes profiles,and the superindex i labels the spectral lines. The previ-ous merit functions consider the quite general case in whichthe standard deviation of the corrupting Gaussian noise isdifferent for Stokes Q , U , and V and for each wavelengthpoint x j . However, in practice, we have the simpler situa-tion in which σ i V j ≈ σ i Q j ≈ σ i U j = σ ij . Moreover, althoughthe number of photons arriving in the line cores are smallerthan in the far wings, we make the approximation that thenoise variance is wavelength-independent. Furthermore, wealso consider that it is independent of the considered line,so that σ ij = σ . These simplifications lead to less clutteredexpressions for the inferred parameters. However, we pointout that the general expressions that emerge from the opti-mization of Eq. (12) can be found in the Appendix.In order to infer a certain parameter, we have to findthe global minimum of the χ . This is trivially obtainedby solving the non-linear system of equations obtained bysetting the derivatives of the χ function with respect tothat parameter to zero. By so doing (see Appendix), we canobtain the expression of the magnetic field vector in termsof the observables: B || = − C P ij V ij I ′ ij P ij ( I ′ ij ) B ⊥ = 1 C q ( P ij Q ij I ′′ ij ) + ( P ij U ij I ′′ ij ) P ij ( I ′′ ij ) φ = 12 arctan P ij U ij I ′′ ij P ij Q ij I ′′ ij + φ , (13) and the derived quantities:tan θ = B ⊥ B || B = q B k + B ⊥ . (14)The phase shift φ is used to set the correct quadrant forthe azimuth and depends on the sign of the numerator U = P ij U ij I ′′ ij and the denominator Q = P ij Q ij I ′′ ij as fol-lows. If Q = 0 then: φ = U > Q > π if U < Q > π if Q < Q = 0 then: φ = ( π if U > π if U < Q and U are zero, the angle φ is obviously undefined.The estimated errors can be computed from the covari-ance matrix (e.g., Press et al. 1986). In the simple case weconsider of equal wavelength-independent standard devia-tion for all spectral lines, the covariance matrix is diagonal(no correlation between the parameters), and the errors ata confidence level of 68.3% (one sigma) are expressed as: δ B || = ± σC qP ij ( I ′ ij ) δ B ⊥ = ± σ C B ⊥ qP ij ( I ′′ ij ) δφ = ± σ C B ⊥ qP ij ( I ′′ ij ) δθ = ± q B ⊥ δ B k + B k δ B ⊥ B k + B ⊥ δ B = ± vuut B k δ B k + B ⊥ δ B ⊥ B k + B ⊥ , (17) It is widely known that all maximum likelihood estimatorsmay suffer from biases. The bias is the difference between thevalue of the estimator and the true value of the parameter.Each individual parameter has to be studied separately us-ing analytical/numerical simulations to understand to whichextent the estimation of B k , B ⊥ , and φ are subjected to bias.For simplicity the simulations we present in the following re-fer to the resolved case in which we use specific intensities.However, the behaviour of the estimator is the same also forthe stellar case.In order to carry out the simulations, we focus on the Fe i line with central wavelength λ = 5250 . D − D . The values ofthe effective Land´e factor for circular polarisation is g = 3,while it is G = 9 for linear polarisation. We consider a con-stant magnetic field strength of B = 500 G, which is suffi-ciently weak so that the weak-field approximation can still c (cid:13)000
110 15+ u − u − v Λ g ∂ F I ( x ) ∂x I ′′ Λ G ∂ I ( x ) ∂x − u − v − u − v Λ G ∂ F I ( x ) ∂x B k B k H k B ⊥ B ⊥ H ⊥ φ φ B φ d θ θ B θ d B B H d observations of a set of spectral lines ( I i , Q i , U i , V i ) with i = 1 . . . n lines . Assuming that the weak field approximationcan be applied for the set of observed spectral lines andthat observations are corrupted with uncorrelated Gaussiannoise, we can use a least-squares estimator (maximum likeli-hood) to retrieve the magnetic field vector (see the Appendixfor more details). The χ merit function is defined as: χ = X ij ( V ij − V i ; modj ) ( σ i V j ) + X ij ( Q ij − Q i ; modj ) ( σ i Q j ) + X ij ( U ij − U i ; modj ) ( σ i U j ) , (12)where the subindex j indicates the position along the coordi-nate x . The label mod refer to the model for Stokes profiles,and the superindex i labels the spectral lines. The previ-ous merit functions consider the quite general case in whichthe standard deviation of the corrupting Gaussian noise isdifferent for Stokes Q , U , and V and for each wavelengthpoint x j . However, in practice, we have the simpler situa-tion in which σ i V j ≈ σ i Q j ≈ σ i U j = σ ij . Moreover, althoughthe number of photons arriving in the line cores are smallerthan in the far wings, we make the approximation that thenoise variance is wavelength-independent. Furthermore, wealso consider that it is independent of the considered line,so that σ ij = σ . These simplifications lead to less clutteredexpressions for the inferred parameters. However, we pointout that the general expressions that emerge from the opti-mization of Eq. (12) can be found in the Appendix.In order to infer a certain parameter, we have to findthe global minimum of the χ . This is trivially obtainedby solving the non-linear system of equations obtained bysetting the derivatives of the χ function with respect tothat parameter to zero. By so doing (see Appendix), we canobtain the expression of the magnetic field vector in termsof the observables: B || = − C P ij V ij I ′ ij P ij ( I ′ ij ) B ⊥ = 1 C q ( P ij Q ij I ′′ ij ) + ( P ij U ij I ′′ ij ) P ij ( I ′′ ij ) φ = 12 arctan P ij U ij I ′′ ij P ij Q ij I ′′ ij + φ , (13) and the derived quantities:tan θ = B ⊥ B || B = q B k + B ⊥ . (14)The phase shift φ is used to set the correct quadrant forthe azimuth and depends on the sign of the numerator U = P ij U ij I ′′ ij and the denominator Q = P ij Q ij I ′′ ij as fol-lows. If Q = 0 then: φ = U > Q > π if U < Q > π if Q < Q = 0 then: φ = ( π if U > π if U < Q and U are zero, the angle φ is obviously undefined.The estimated errors can be computed from the covari-ance matrix (e.g., Press et al. 1986). In the simple case weconsider of equal wavelength-independent standard devia-tion for all spectral lines, the covariance matrix is diagonal(no correlation between the parameters), and the errors ata confidence level of 68.3% (one sigma) are expressed as: δ B || = ± σC qP ij ( I ′ ij ) δ B ⊥ = ± σ C B ⊥ qP ij ( I ′′ ij ) δφ = ± σ C B ⊥ qP ij ( I ′′ ij ) δθ = ± q B ⊥ δ B k + B k δ B ⊥ B k + B ⊥ δ B = ± vuut B k δ B k + B ⊥ δ B ⊥ B k + B ⊥ , (17) It is widely known that all maximum likelihood estimatorsmay suffer from biases. The bias is the difference between thevalue of the estimator and the true value of the parameter.Each individual parameter has to be studied separately us-ing analytical/numerical simulations to understand to whichextent the estimation of B k , B ⊥ , and φ are subjected to bias.For simplicity the simulations we present in the following re-fer to the resolved case in which we use specific intensities.However, the behaviour of the estimator is the same also forthe stellar case.In order to carry out the simulations, we focus on the Fe i line with central wavelength λ = 5250 . D − D . The values ofthe effective Land´e factor for circular polarisation is g = 3,while it is G = 9 for linear polarisation. We consider a con-stant magnetic field strength of B = 500 G, which is suffi-ciently weak so that the weak-field approximation can still c (cid:13)000 , 1–12 aximum likelihood estimation of stellar magnetic fields Figure 2.
Bias of the maximum likelihood estimator of the magnetic field in the weak field approximation. Different color lines representdifferent values of the noise level added to the synthetic profiles of the Fe i line at 5250.2 ˚A . The color code is the following: yellowrepresents the inversion of the profiles without noise. Orange displays the results for a noise level of 10 − I c , pink for 5 × − I c , red for10 − I c , brown for 5 × − I c , dark green for 10 − I c , light green for 1 . × − I c , light blue for 2 × − I c , dark blue for 5 × − I c , light violet for 10 − I c , and dark violet for 5 × − I c . The I c is the continuum intensity. be considered. The inclination and the azimuth of the fieldare set to vary uniformly between 0 and 180 ◦ . Letting theazimuth vary in this interval, we avoid the 180 ◦ ambiguityin the azimuth present in the radiative transfer equation.For simplification, we consider a Gaussian intensity profilein the form: I ( λ ) = 1 − d c exp (cid:20) − ( λ − λ ) w (cid:21) , (18)where d c = 1 / w = 0 .
05 ˚A. The parameters for theGaussian profile have been fixed to fit a solar observationobtained with the IMaX instrument. From the IMaX obser-vational capabilities, it can be verified that the 5250 . ∼ P , i. e., the value of the parameter that contains 50% of the area of the distribution). To quantify the dispersion produced bynoise we use the percentiles 16 and 84 (which encompass onestandard deviation around the estimated value).Figure 2 displays the inversion of the synthetic profileswith different values of the standard deviation of the noise,from 0 to 5 × − in units of the continuum intensity, I c .It is clear from the upper and lower left panels that boththe longitudinal component of the magnetic field and theazimuthal angle are unbiased quantities. This means thatthe estimated value statistically coincides with the originalone. As expected, the dispersion grows with the increasingnoise level.Contrary to B k and φ B (note that the notation forthe magnetic field vector is the one associated to resolvedsources), the transversal component of the magnetic fieldand hence the inclination angle presents a non-zero bias.Except for the case in which there is no added noise, smalltransversal components of the field (smaller or at the levelof the noise amplitude) are overestimated. The fundamentalreason is that the expression for B ⊥ in Eq. (13) is not ro-bust against noise. It can be verified that if Q ij and U ij are atthe noise level (can be described as Gaussian distributionswith zero mean and variance σ ), B ⊥ is a random variable c (cid:13) , 1–12 M. J. Mart´ınez Gonz´alez following the probability distribution: p ( B ⊥ ) = 2 C (cid:16)P ij ( I ′′ ij ) (cid:17) σ B ⊥ exp " − C P ij ( I ′′ ij ) B ⊥ σ . (19)The value of this distribution for a given percentile c fulfils: B c ⊥ = − − c ) P ij ( I ′′ ij ) ! / √ σC . (20)The percentile 50 ( c = 0 .
5) correctly captures the value ofthe bias at small values of B ⊥ . Once this value is computed,if the inferred value of B ⊥ is similar to this value, one shouldbe aware that the correct value of B ⊥ might be smaller. Ifthis is not taken into account, the estimated inclination ofthe field is larger than the real one and artificially horizontalfields might be inferred. Note also that in the stellar case,the inferred inclination of the dipole axis would be moreinclined than the real one. Most detections of faint signals in stellar atmospheres havebeen evidenced by adding many spectral lines. The polari-metric signal per spectral resolution element is known to bewell below the noise level, so that line addition is a must tofight against photon noise. The most widely used and suc-cessful technique that combines the information of manyspectral lines is the Least-Squares Deconvolution (LSD)technique (Donati et al. 1997). The equations presented inthis paper allow us to retrieve the magnetic field vector fromthe polarised stellar spectra taking into account many spec-tral lines. Thanks to it, we can rewrite the equations todirectly extract information from the LSD profiles.The LSD technique is fundamentally based on the appli-cation of the weak-field approximation and on the assump-tion of a constant CLV for all spectral lines (in our case, u and v are set constant). This means that all Stokes profilesfor each spectral line can be computed with a single spectralprofile whose proportionality constant changes from line toline (see Donati et al. 1997): I i = η i ¯ I (21) V i = η i Λ i g i ¯ V (22) Q i = η i (Λ i ) G i ¯ Q (23) U i = η i (Λ i ) G i ¯ U , (24)where η i is the line depth and ¯ I , ¯ Q , ¯ U , and ¯ V are the LSDprofiles that can be computed using a least-squares pro-cedure as explained in Donati et al. (1997). We refer toKochukhov, Makaganiuk & Piskunov (2010) for anin-depth analysis of the assumptions and potentialproblems of LSD.
From the previous equations, it is easyto show that it is possible to infer the magnetic field vector using: B k = − CK P j ¯ V j ∂ ¯ I j ∂x P j (cid:16) ∂ ¯ I j ∂x (cid:17) B ⊥ = 1 C K ′ r(cid:16)P j ¯ Q j ∂ ¯ I j ∂x (cid:17) + (cid:16)P j ¯ U j ∂ ¯ I j ∂x (cid:17) P (cid:16) ∂ ¯ I j ∂x (cid:17) φ = 12 arctan P j ¯ U j ∂ ¯ I j ∂x P j ¯ Q j ∂ ¯ I j ∂x + φ , (25)where K = 1 and K ′ = 1 / K = 110 15 + u − u − v , K ′ = 14480 420 − u − v − u − v (26)for the stellar dipole. Note that it is not necessary to assumeaveraged atomic parameters for the LSD profile, treating itas a mean spectral line. In our case, the estimated magneticfield vector only depends on the observables, the atomic pa-rameters of each observed spectral line (which is necessaryto compute the LSD profile), and the assumed CLV coeffi-cients. The inference power of the expressions developed inthe previous sections are illustrated with the aid oftwo different examples. The first one consists of asimulated stellar dipole using a Milne-Eddington (e.g.,Landi Degl’Innocenti & Landolfi 2004) atmosphere and thesecond one is a particularly interesting observational exam-ple in which we can illustrate the effect of the bias in thetransverse component of the magnetic field with spatial res-olution.Considering the stellar dipole, we simulate a Milne-Eddington atmosphere with a source function that variesas: S ( τ ) = S (1 + βτ ) . (27)For simplicity, we choose β = 1 and consider a static atmo-sphere at all points of the stellar surface. We assume a spec-tral line centred at λ = 5000 ˚A with a Doppler broadening of0.04 ˚A . Both the circular and linear effective Land´e factorsare equal to 1. The CLV variation of the Milne-Eddingtonatmosphere (that we have forced to be wavelength indepen-dent) gives u = 2 /
3. We integrate the Stokes signals on thevisible stellar surface and we add Gaussian noise to the finalintegrated flux with a standard deviation of 5 × − in unitsof the intensity flux at the continuum, F c .Figure 3 shows the flux of the simulated Stokes param-eters coming from a dipole field with the following parame-ters: H d = 1500 G, θ d = 80 ◦ , and φ d = 25 ◦ . The black linesrepresent the synthetic fluxes, and the rhombs the syntheticobservations with added noise. We invert the noisy Stokesfluxes using Eqs. 13, 14 and 17 and obtain H || = 265 . ± . H ⊥ = 1501 . ± . φ d = 27 . ◦ ± .
0, and the derivedquantities H d = 1524 . ± . θ d = 80 . ◦ ± .
2. Allquantities are nicely recovered. In fact, the bias estimations c (cid:13) , 1–12 aximum likelihood estimation of stellar magnetic fields Figure 3.
Inversion of a synthetic stellar dipole. The black lines display the synthetic flux of a dipole with a strength H d =1500 G, aninclination θ d = 80 ◦ with respect to the line of sight and an azimth of θ d = 25 ◦ . The rhombs are the synthetic polarised flux includingan additive Gaussian noise with a standard deviation of 5 × − F c . The blue lines are the fluxes inferred from the inversion using Eq.13. for the perpendicular component H ⊥ for a the percentiles16, 50, and 84 (i. e. containing one sigma probability) are302.0, 426.5, and 543.8 G, respectively. This means that thecomputed H ⊥ is reliable. Now, we assume a weak dipolewith H d = 100 G, θ d = 20 ◦ , and φ d = 25 ◦ , and the samenoise level. The inferred parameters are H || = 99 . ± . H ⊥ = 333 . ± . φ d = 39 . ◦ ± . H d = 348 . ± . θ d = 73 . ◦ ± .
1. In principle, both the longitudinaland the transverse components (as well as the inclination)should be well recovered but the azimuth remains undeter-mined. However, the bias of H ⊥ for the percentiles 16, 50,and 84 are 293.2, 414.1, and 528.0 G, respectively. There-fore, the inferred perpendicular component of the dipole isconsistent with a bias. This implies that the inferred value,although having a small (Gaussian) error, it has to be con-sidered an upper limit. In this case, we know that the per-pendicular component of the dipole is very small (34 G).Consequently, we overestimate both the perpendicular com-ponent and the intrinsic strength of the dipole. Additionally,the dipole appears to be much more inclined that in reality.The second experiment consists of filter-polarimetricdata observed with the IMaX instrument on the Sunrisemission. This polarimeter observes the 5250.2 Fe i line (forwhich we have carried out the bias experiment in Section4). The data set consists of the four Stokes parameters ob-served at the quietest areas of the solar disc centre at a spa-tial cutoff frequency of about 0.15-0.18 ′′ ( ∼
120 km in thesolar surface, the best one at the moment for instrumentswith polarimetric capabilities). The noise level in circular and linear polarization is 10 − in units of the continuumintensity, I c . The left panel of Fig. 4 shows the continuumintensity image, with brighter areas associated to granularregions, where the plasma ascends to the photosphere. Darkareas are the intergranular lanes, where the motions are pref-erentially downflowing. The right panel of this same figuredisplays the estimated bias for the transversal component ofthe magnetic field using Eq. 20 for the percentile c = 0 . c (cid:13) , 1–12 M. J. Mart´ınez Gonz´alez
Figure 4.
Left: continuum intensity at 227 m˚A from the core of the 5250.2 ˚A. Right: Median value of the bias for the transversecomponent of the magnetic field computed with Eq. 20 and c = 0 . Figure 5.
Left (centre): inferred value of the longitudinal (transverse) magnetic field for the 5250.2 ˚A line. Right: transverse fieldcorrected from the bias due to pure noise. Central and right panels share the same colorbar. e.g., Orozco Su´arez et al. 2007; Lites et al. 2008; Sheminova2009; Ishikawa & Tsuneta 2009, 2010, 2011).Since we can characterize the bias by its median value,it is possible to distinguish real signals from false signals pro-duced by the presence of noise. The way we proceed is asfollows. We compute a conservative upper limit for a trust-ful field as the bias for the percentile 84 ( c = 0 .
84 in Eq.20). This value changes from pixel to pixel. Then, we force B ⊥ = 0 in those places where the inferred B ⊥ is smaller thanthe bias. The result is represented in the right panel of Fig. 5.Now, the real signals (coming from pixels with linear polar-ization clearly above the noise level) are much more evidentand most of the background has disappeared. Note that wehave only removed those signals that were produced by thepresence of noise (the expected value of B ⊥ = 0 should bezero). However, if B ⊥ = 0, the bias can still be important ifthe noise level is high, as shown in Fig. 2. We have shown that the weak field approximation (theStokes parameters are proportional to the first and secondorder derivative of the intensity) also holds for observed stel-lar fluxes in the case of slow rotators. We have used a max-imum likelihood estimator to infer the magnetic field vectorfrom the observed Stokes profiles. The main result of this papers is that we give explicit formulae for the componentsof the magnetic field vector in terms of the observables. Theformulae are general and hold for specific intensities, for in-tegrated fluxes and are slightly modified for LSD profiles.In the particular case of a stellar dipole, the orientation ofthe dipole axis and its intensity can be recovered from oneobservation of the full Stokes vector.We have also studied the bias of this maximum like-lihood estimator. The longitudinal magnetic field and theazimuthal angle are unbiased quantities. However, thetransversal component of the magnetic field and hence theinclination of the field are overestimated in the presence ofnoise. We derive the estimated value of the perpendicularcomponent of the magnetic field in the case that there isno linear polarization signal above the noise (the bias for B ⊥ =0). We propose to evaluate this bias prior to the infer-ence of the perpendicular component of the field. One shouldbe very cautious when the inferred B ⊥ is of the order of thebias. ACKNOWLEDGMENTS
We are grateful to F. Leone for helpful comments. Thiswork has been funded by the Spanish Ministry of Sci-ence and Innovation under projects AYA2010-18029 (So- c (cid:13) , 1–12 aximum likelihood estimation of stellar magnetic fields lar Magnetism and Astrophysical Spectropolarimetry) andConsolider-Ingenio 2010 CSD2009-00038. REFERENCES
Asensio Ramos A., 2011, ApJ, 731, 27Asensio Ramos A., Manso Sainz R., 2011, ApJ, 731, 125Asensio Ramos A., Trujillo Bueno J., Landi Degl’InnocentiE., 2008, ApJ, 683, 542Aznar Cuadrado R., Jordan S., Napiwotzki R., SchmidH. M., Solanki S. K., Mathys G., 2004, A&A, 423, 1081Bagnulo S., Szeifert T., Wade G. A., Landstreet J. D.,Mathys G., 2002, A&A, 389, 191Claret A., 2000, A&A, 363, 1081Cox A. N., 2000, Allen’s Astrophysical Quantities, 4th ed.Springer Verlag and AIP Press, New YorkDonati J.-F., Semel M., Carter B. D., Rees D. E., CollierCameron A., 1997, MNRAS, 291, 658Ishikawa R., Tsuneta S., 2009, A&A, 495, 607—, 2010, ApJL, 718, L171—, 2011, ApJ, 735, 74Jordan S., Werner K., O’Toole S. J., 2005, A&A, 432, 273Kochukhov O., Makaganiuk V., Piskunov N., 2010, A&A,524, A5+Landi Degl’Innocenti E., Landolfi M., 2004, Polarization inSpectral Lines. Kluwer Academic PublishersLandolfi M., Landi degl’Innocenti E., Landi degl’InnocentiM., Leroy J. L., 1993, 272, 285Leone F., Mart´ınez Gonz´alez M. J., Corradi R. L. M., Priv-itera G., Manso Sainz R., 2011, ApJL, 731, L33+Lites B. W. et al., 2008, ApJ, 672, 1237Mart´ınez Pillet V. et al., 2011, Sol. Phys., 268, 57Merenda L., Trujillo Bueno J., Landi Degl’Innocenti E.,Collados M., 2006, ApJ, 642, 554Orozco Su´arez D. et al., 2007, ApJL, 670, L61O’Toole S. J., Jordan S., Friedrich S., Heber U., 2005, A&A,437, 227Press W. H., Teukolsky S. A., Vetterling W. T., FlanneryB. P., 1986, Numerical Recipes. Cambridge UniversityPress, CambridgeSheminova V. A., 2009, Astronomy Reports, 53, 477Silvester J. et al., 2009, MNRAS, 398, 1505Solanki S. K. et al., 2010, ApJL, 723, L127Spirock T. J. et al., 2001, AGU Spring Meeting Abstracts,51Varsik J. R., 1995, Sol. Phys., 161, 207Wade G. A., Donati J., Landstreet J. D., Shorlin S. L. S.,2000, MNRAS, 313, 851This paper has been typeset from a TEX/ L A TEX file preparedby the author. c (cid:13) , 1–12 M. J. Mart´ınez Gonz´alez
APPENDIX A: DERIVATION OF THE MAXIMUM LIKELIHOOD SOLUTION OF THE STELLARMAGNETIC FIELD
We start from the general weak-field equations that hold both for the resolved (solar) case and for the integrated stellar dipole(Eq. 11): V = − C B k I ′ Q = − C B ⊥ cos 2 φ I ′′ U = − C B ⊥ sin 2 φ I ′′ , (A1)where the explicit expressions for V , Q , U , I ′ , and I ′′ , for each case can be found in Table 1 in the main text. We denotethe Stokes vector by S = [ I , Q , U , V ]. Assuming that the difference between the data and the model follows a Gaussiandistribution, the likelihood (the probability distribution of the data given the parameters) can we written as: L = n l Y i =1 n λ Y j =1 4 Y k =2 h π ( σ ijk ) i − / exp " − ( S ijk − S i ;mod jk ) σ ijk ) , (A2)where the index i refers to the spectral line, j refers to the wavelength, and k to the element of the Stokes vector. Theabbreviation “mod” stands for the model. The symbol σ stands for the variance of S − S mod . Taking the logarithm, we buildthe log-likelihood, ln L : ln L = − n l n λ π ) − n l X i =1 n λ X j =1 4 X k =2 " ( S ijk − S i ;mod jk ) σ ijk ) + 12 ln( σ ijk ) . (A3)In order to estimate the parameters that fit the data given the proposed model we have to maximise the likelihood or,equivalently, minimise ln L . Let p = (cid:2) B k , B ⊥ , φ (cid:3) denote the set of parameters of our model. Following the standard approach,to estimate the vector of parameters p , we must solve the following set of equations: ∂ ln L ∂p l = 0 = n l X i =1 n λ X j =1 4 X k =2 ∂∂p l " ( S ijk − S i ;mod jk ) ( σ ijk ) + ln( σ ijk ) , ∀ l. (A4)There is an important point to clarify given that, in our case, the model is not fully analytical because it depends on theobserved intensity profile (Asensio Ramos & Manso Sainz 2011). Therefore, assuming no correlation between the model andthe observable (the only difference being produced by uncorrelated Gaussian noise), the variance for each Stokes parameteris given by: ( σ ) = σ ( S ) + σ ( S mod2 ) = σ ( S ) + C B ⊥ cos φ σ ( I ′′ ) (A5)( σ ) = σ ( S ) + σ ( S mod3 ) = σ ( S ) + C B ⊥ sin φ σ ( I ′′ ) (A6)( σ ) = σ ( S ) + σ ( S mod4 ) = σ ( S ) + C B k σ ( I ′ ) . (A7)Thanks to the previous equations, the variances that appear in Eq. (A4) depend on the actual parameters, which makes theminimisation much harder. Luckily, in the weak field regime and for the observational spectral resolution of interest nowadays,it is easy to verify that σ ( S mod i ) σ ( S i ). In any case, this condition should be checked before carrying out any inversionusing the formulae derived in this work. Assuming a first order approximation to the derivative I ′ , we find: σ ( I ′ ) = 2 K σ ( S )∆ x , (A8)where K = 1 for the resolved case and K =
110 15+ u − u − v for the dipole case. Assuming σ ( S mod i ) σ ( S i ) and that the noisein intensity is the same as in the circular polarisation, it holds that:∆ x > √ C B k K Λ g, (A9)which means that the spectral sampling has to be larger than or of the order of the Zeeman splitting. For instance, for awavelength of 5000 ˚A, a Land´e factor of 1.5, and B || = 500 G, the spectral sampling has to be larger or equal to 12 m˚Afor K = 1. Luckily, this is the case in most of the observational cases. For example, for the same wavelength, the expectedsampling for two spectrographs with a resolving power of R = 60000 , x > p GK ′ cos 2 φ C B ⊥ Λ∆ x > p GK ′ sin 2 φ C B ⊥ Λ , (A10)where K ′ = 1 / K ′ = − u − v − u − v for the dipole case. For B ⊥ = 500 G, φ = 0 ◦ and G = 2 . x = 6 m˚A. c (cid:13) , 1–12 aximum likelihood estimation of stellar magnetic fields Taking the previous considerations into account, Eq. (A4) can be simplified to: ∂ ln L ∂p l = 0 = n l X i =1 n λ X j =1 4 X k =2 ∂∂p l " ( S ijk − S i ;mod jk ) ( σ ijk ) ≡ ∂χ ∂p l , (A11)where the well known χ merit function is defined as: χ = X ij ( V ij − V i ;mod j ) ( σ i V j ) + X ij ( Q ij − Q i ;mod j ) ( σ i Q j ) + X ij ( U ij − U i ;mod j ) ( σ i U j ) . (A12)Explicitly, the derivatives of the χ with respect to the parameters we want to infer are: ∂χ ∂ B k = 2 C X ij V ij I ′ ij ( σ i V j ) + 2 C B k X ij ( I ′ ij ) ( σ i V j ) (A13) ∂χ ∂ B ⊥ = 4 C B ⊥ cos 2 φ X ij Q ij I ′′ ij ( σ i Q j ) + sin 2 φ X ij U ij I ′′ ij ( σ i U j ) ! + 4 C B ⊥ cos φ X ij ( I ′′ ij ) ( σ i Q j ) + sin φ X ij ( I ′′ ij ) ( σ i U j ) ! (A14) ∂χ ∂φ = 4 C B ⊥ cos 2 φ X ij U ij I ′′ ij ( σ i U j ) − sin 2 φ X ij Q ij I ′′ ij ( σ i Q j ) ! + 4 C B ⊥ sin 2 φ cos 2 φ X ij ( I ′′ ij ) ( σ i U j ) − X ij ( I ′′ ij ) ( σ i Q j ) ! (A15)By forcing these derivatives to zero, we obtain the maximum-likelihood estimate of each parameter. For the longitudinalmagnetic field, using Eq. A13, we obtain a unique solution: B k = − C P ij V ij I ′ ij ( σ i V j ) P ij ( I ′ ij ) ( σ i V j ) . (A16)Dividing Eq. A14 and A15, we obtain a solution for the azimuth:tan 2 φ = P ij ( I ′′ ij ) ( σ i Q j ) P ij U ij I ′′ ij ( σ i U j ) P ij ( I ′′ ij ) ( σ i U j ) P ij Q ij I ′′ ij ( σ i Q j ) . (A17)Dividing Eq. A15 by cos 2 φ and using Eq. A17, after some algebra we obtain two solutions. One is B ⊥ = 0, which is not validsince this solution maximises the χ . We have tested this computing the second derivative. The solution that minimises the χ is given by: B ⊥ = 1 C s(cid:18)P ij ( I ′′ ij ) ( σ i U j ) P ij Q ij I ′′ ij ( σ i Q j ) (cid:19) + (cid:18)P ij ( I ′′ ij ) ( σ i Q j ) P ij U ij I ′′ ij ( σ i U j ) (cid:19) P ij ( I ′′ ij ) ( σ i U j ) P ij ( I ′′ ij ) ( σ i Q j ) . (A18)The errors associated with each parameter are computed assuming that the surface around the maximum likelihood isapproximately a multidimensional Gaussian. This is equivalent to assuming a parabolic approximation to the χ close to theminimum. The curvature close to the minimum is given by the Hessian matrix ζ , whose elements are: ζ kl = 12 ∂ χ ∂p k ∂p l . (A19) c (cid:13) , 1–12 M. J. Mart´ınez Gonz´alez
The matrix can be calculated using the following second derivatives: ∂ χ ∂ B k = 2 C X ij ( I ′ ij ) ( σ i V j ) (A20) ∂ χ ∂ B ⊥ = 8 C B ⊥ cos φ X ij ( I ′′ ij ) ( σ i Q j ) + sin φ X ij ( I ′′ ij ) ( σ i U j ) ! (A21) ∂ χ ∂φ = 8 C B ⊥ sin φ X ij ( I ′′ ij ) ( σ i Q j ) + cos φ X ij ( I ′′ ij ) ( σ i U j ) ! (A22) ∂χ ∂ B k ∂ B ⊥ = ∂χ ∂ B ⊥ ∂ B k = 0 (A23) ∂χ ∂ B k ∂φ = ∂χ ∂φ∂ B k = 0 (A24) ∂χ ∂ B ⊥ ∂φ = ∂χ ∂φ∂ B ⊥ = 8 C B ⊥ sin 2 φ cos 2 φ X ij ( I ′′ ij ) ( σ i U j ) − X ij ( I ′′ ij ) ( σ i Q j ) ! . (A25)The square root of the diagonal of the covariance matrix gives the error estimates for each parameter. Such covariancematrix is just the inverse of the Hessian matrix, C = ζ − . In our case: C = C C C C C , (A26)with C = " C X ij ( I ′ ij ) ( σ i V j ) − (A27) C = 14 C B ⊥ sin φ P ij ( I ′′ ij ) ( σ i Q j ) + cos φ P ij ( I ′′ ij ) ( σ i U j ) cos φ sin φ (cid:18)P ij ( I ′′ ij ) ( σ i Q j ) − P ij ( I ′′ ij ) ( σ i U j ) (cid:19) + P ij ( I ′′ ij ) ( σ i Q j ) P ij ( I ′′ ij ) ( σ i U j ) (A28) C = C = 14 C B ⊥ sin 2 φ cos 2 φ (cid:18)P ij ( I ′′ ij ) ( σ i U j ) − P ij ( I ′′ ij ) ( σ i Q j ) (cid:19) cos φ sin φ (cid:18)P ij ( I ′′ ij ) ( σ i Q j ) − P ij ( I ′′ ij ) ( σ i U j ) (cid:19) + P ij ( I ′′ ij ) ( σ i Q j ) P ij ( I ′′ ij ) ( σ i U j ) (A29) C = 14 C B ⊥ cos φ P ij ( I ′′ ij ) ( σ i Q j ) + sin φ P ij ( I ′′ ij ) ( σ i U j ) cos φ sin φ (cid:18)P ij ( I ′′ ij ) ( σ i Q j ) − P ij ( I ′′ ij ) ( σ i U j ) (cid:19) + P ij ( I ′′ ij ) ( σ i Q j ) P ij ( I ′′ ij ) ( σ i U j ) (A30)The error bars are thus given by: δp k = ± ∆ √C kk , (A31)with ∆ = 1 , . , , . , , and 3 .
89 for a confidence level of 68.3%, 90%, 95.4%, 99%, 99.73%, and 99.99%, respectively. Thiscomputation assumes that, to estimate the error of a parameter, we fix the values of the rest to the ones that maximize thelikelihood and compute the confidence levels for the one-dimensional probability distribution of this parameter (see Press et al.1986, for more details). Note that this approximation does not take into account the degeneracies between parameters sinceit does not integrate the probability distribution of the rest of the parameters but fixes a certain value (see Asensio Ramos2011, for a robust Bayesian inversion). Note also that the covariance matrix is diagonal if the standard deviation of Q and U are the same ( σ i Q j = σ i U j ). This is generally the case in solar observations in which the measurement efficiencies for linearpolarization are similar. c (cid:13)000