Jerome Thiebaut

Simon Prunet

Christophe Pichon

Eric Thiébaut

J Thiébaut

Probing magnetic fields with multi-frequency polarized synchrotron emission

published or not. The documents may come

INTRODUCTION

The problem of studying the magnetic field structure of our Galaxy using measurements of the synchrotron emission of high energy electrons in the Galactic magnetic field is an old one (Ginzburg & Syrovatskii 1965;Ruzmaikin et al. 1988;Beck et al. 1996). The fact that the emitting medium is itself magnetized induces a differential Faraday rotation of the different emission planes transverse to the line of sight, resulting in a well known depolarization effect of the integrated emission that depends strongly on the electromagnetic frequency. This effect, described in the first place by Burn (1966) in the case of a constant magnetic field, has been further studied in semi-analytically for given functional forms of the magnetic field; it has also been studied from the statistical point of view in some asymptotic regimes (see e.g. Sokoloff et al. 1998). In the present work, we want to consider the more ambitious problem of using this depolarization effect, together with the solenoidal character of the magnetic field, to reconstruct the magnetic field structure from a set of polarized maps of the synchrotron emission of an ionized medium at different electromagnetic frequencies. With the upcoming prospect of detailed Multiband spectroscopy in the radio wavelengths (Röttgering 2003;Furlanetto & Briggs 2004), developed in the context of Galactic and high-redshift extragalactic HI lines, this type of investigation should become possible.

A statistical inference of the measurement of the Galactic magnetic field correlator as a function of scale from multifrequency polarization measurements has already been successfully achieved by [START_REF] Vogt | [END_REF] in the case of the Faraday rotation of the polarized light from background objects by the intra-cluster magnetized plasma. In this case, there is no depolarization effect due to differential Faraday rotation, and the relationship between the measured polarization at a given frequency and the polarization of light in the source plane is linear in the (longitudinal) magnetic field strength. The linearity of the problem makes the statistical analysis tractable in the former case. In the case that we investigate, the emitting and the rotating medium are the same, which results in depolarization effects of the emitted light. Moreover, the synchrotron emissivity itself depends non-linearly on the field strength transverse to the line of sight. The reconstruction of the magnetic field structure from the polarization data is in this case a non-linear inverse problem. Finally, we must note that to address the full problem of reconstruction of the magnetic field from the depolarized synchrotron emission we need in principle knowledge of both the thermal electron spatial distribution ne and the spatial distribution of cosmic ray electrons nr, when, in comparison, the inference of the magnetic energy spectrum from the rotation measures of background sources only requires knowledge of the thermal electron distribution.

In a first attempt at reconstructing the magnetic field, and for the sake of clarity, we make the assumption that the fluctuations of the thermal and cosmic ray electrons can be neglected compared to the fluctuations in the magnetic field itself. This assumption, if physically unrealistic, allows us to show the specific influence of the magnetic field statistical properties on the quality of the reconstruction. In the first sections, we thus consider the electronic distributions (both thermal and relativistic) as constant, and discuss the reconstruction of the magnetic field using only the leading coupling coefficient in the equation of radiative transfer. In the (thin medium, strong rotativity) limit that we assume for this work, this leading term is the usual Faraday term, responsible for the rotation of the plane of polarization. We will assume that the Faraday coefficient is dominated by the thermal electrons, which is a reasonable assumption in non-relativistic astrophysical plasmas. Finally, in section 4, we relax the unrealistic assumption of a constant thermal electrons density, and show that our method can still be used to reconstructed the magnetic field when the electronic density is spatially varying but known a priori, using simulated data sets from a magneto-hydrodynamical (MHD) simulation.

This paper is organized as follows: in section 2 we discuss the fonctional dependence of the polarization of the synchrotron emission and its variation with electro-magnetic frequency on the underlying magnetic field. We present a discretized version of this functional dependence that will be useful in the context of the reconstruction from discrete polarization data. In section 3 we investigate the reconstruction of the magnetic field from simulated multifrequency polarized data, when the functional dependence on the magnetic field has been linearized around a "mean" field. Taking advantage of the linear nature of this approximate problem, we give a strategy for choosing the best electromagnetic frequencies of observation, and investigate the statistical anisotropy of the magnetic field reconstruction errors. Finally, in section 4, we investigate the validity of the linearization procedure used in the precedent section, as a function of the quality of our prior knowledge of the magnetic field structure. We show how the approximate, linearized inverse problem investigated in this work could be used as a building block of the fully non-linear reconstruction problem. We emphasize that any gradient-based non-linear minimization algorithm can be decomposed into linear sub-problems, thus justifying the study of the linearized problem. In this context, we investigate how the conditioning of the linearized problem varies with the properties of the reference magnetic field around which the problem is being linearized. In particular it is illustrated on a realistic reference field from a MHD simulation. Finally, using the same MHD simulation data, we show that our method can deal with a non-constant electronic density, provided it is known a priori. In section 5, we summarize the main results of the paper, recalling the main simplifying assumptions used to derive them (notably the assumed-known electronic density hypothesis) and discuss how this assumption could be possibly alleviated by additional data (e.g., Hα, free-free) or by using second-order coupling terms involving the circular polarization in the case of relativistic sources (see C). We conclude on how the different results of the paper could be used to tackle the fully non-linear reconstruction of the magnetic field.

POLARIZED EMISSION

Our objective is to recover the magnetic field given observed polarization maps at different wavelengths. We tackle this ill-posed problem by means of an inverse problem approach [START_REF] Tarantola | Inverse Problem Theory[END_REF] which involves recovering the magnetic field that gives a polarization consistent with the observations while obeying some a priori properties. These priors are strict constraints, such as ∇ • B = 0, to insure that the sought field is physically meaningful and a regularization to lever the degeneracies of the inverse problem while avoiding artifacts due to noise amplification. We first derive the direct model of the polarization given the magnetic field and then introduce the inverse problem approach in a Bayesian framework.

Direct model

We only consider here the Faraday rotation in the transfer equation, and neglect all other coupling terms. In this case, the transfer equation of the Stokes parameters of linear polarization (Q, U ) can be integrated formally. We assume here that the density of electrons is constant, or that its fluctuations are only important on scales that are not considered here.

Consider a slab of ionized magnetized medium of width L which is emitting synchrotron radiation. The polarized emission, as a function of frequency, integrated over the line of sight then reads (Sokoloff et al. 1998):

P ≡ Q + i U = Z ǫ(r) e 2 i ψ(r) dz , (1)
with Q and U are the usual Stokes parameters, ǫ(r) the synchronton emissivity which obeys:

ǫ(r) = A nr(r) |B ⊥ (r)| γ+1 2 ν -γ-1 2 , (2)
and ψ(r) the sum of the Faraday rotation and the primordial orientation:

ψ(r) = π/2 + arctan (By/Bx) + K ν 2 Z 0 z neBz dz ′ , (3)
where r ≡ (x, y, z) = (x ⊥ , z) is the coordinate in the slab, ν is the frequency, and B = (Bx, By, Bz) = (B ⊥ , Bz) is the magnetic field. In equation ( 3), K reads:

K = q 3 e 8 π 2 m 2 e c ǫ0 , (4)
while, in equation ( 2), A is given by

A = √ 3 E γ 0 q 3 e 16 π ǫ0 me c " 3 qe 2 π m 3 e c 4 « γ-1 2 Γ " 3 γ -1 12 « Γ " 3 γ + 1 12 « ,
where E0 is the energy scale of the relativistic electron spectrum, me and qe stand for the mass and the charge of the electron, ne and nr are the thermal and relativistic electron densities supposed constant, while the exponent γ stands for the spectral index of the cosmic ray electrons, c is the speed of light, ǫ0 is the electric permittivity and Γ is the Euler gamma function. The lengths are in kilo-parsec (kpc) and so the density in kpc -3 , the magnetic fields in micro-Gauss (µG) and the frequencies in giga-Hertz (GHz). Reexpressing the intrinsic polarization phase in terms of powers of the magnetic field components, we get the following expression for the polarization:

P (x ⊥ , ν) = A ν -γ-1 2 Z 0 -∞ nr(x ⊥ , z) `B2 x + B 2 y ´γ-3 4 (x ⊥ , z) × `B2 x -B 2 y + 2 i Bx By ´(x ⊥ , z) × exp " 2 i K ν 2 Z 0 z (neBz)(x ⊥ , z ′′ ) dz ′′ « dz . (5)
As real data come in discrete form, let us discretize this expression by replacing all integrals with sums, assuming a regular discretization grid that will be defined more precisely below. Equation (5) then reads

P (x ⊥ , ν) = A h ν -γ-1 2 X z nr(x ⊥ , z) `B2 x + B 2 y ´γ-3 4 (x ⊥ , z) × `B2 x -B 2 y + 2 i Bx By ´(x ⊥ , z) × exp 2 i K h ν 2 X z ′ θH(z ′ -z) (ne Bz)(x ⊥ , z ′ ) ! . (6)
c 0000 RAS, MNRAS 000, 000-000

Probing magnetic fields with multi-frequency polarized synchrotron emission 3

Here θH is the Heaviside function (θH (x) = 1 for x ≥ 0 and 0 elsewhere), and h the discretization length along z. Equation ( 6) is formally a function of B ≡ {(Bx, By, Bz)} r where we use bold symbols to represent the discretized vector fields and r is a triple index spanning the magnetized volume on a regular cubic mesh with cell size h.

The solution to the inverse problem will be obtained by means of minimization of some merit function (as explained in what follows), we therefore need to compute the partial derivatives of the polarization with respect to the magnetic field. Let us first compute the derivatives with respect to the transverse components of the field:

∂P (x ⊥ , ν) ∂Bx(r ′ ) = δD(r -r ′ ) A nr(r ′ ) h ν -γ-1 2 (B 2 x + B 2 y ) γ-7 4 (r ′ ) × » 1 + γ 2 B 3 x + 7 -γ 2 B 2 y Bx + i (γ -1) B 2 x By + 2 i B 3 y - (r ′ ) × exp 2 i K h ν 2 X z ′′ θH(z ′′ -z ′ ) (ne Bz)(x ′ ⊥ , z ′′ ) ! , (7)
with r = (x ⊥ , z), r ′ = (x ′ ⊥ , z ′ ) and δD Dirac's delta function. The derivative with respect to By follows closely, with the square bracket term becoming:

» - 1 + γ 2 B 3 y - 7 -γ 2 B 2 x By + i (γ -1) Bx B 2 y + 2 i B 3 x - (8)
which corresponds to a π/2 rotation in the plane perpendicular to the LOS. We see that in both cases the phase term is unaffected since it is only a function of the longitudinal magnetic field component Bz. Finally let us compute the derivative with respect to Bz:

∂P (x ⊥ , ν) ∂Bz(r ′ ) = δD(x ⊥ -x ′ ⊥ )2iKAh 2 ν -γ+3 2 × X z nr(x ′ ⊥ , z)(B 2 x + B 2 y ) γ-3 4 (B 2 x -B 2 y + 2iBxBy)(x ′ ⊥ , z) ×θH(z ′ -z) exp 2iKh ν 2 X z ′′ θH(z ′′ -z)( neBz)(x ⊥ , z ′′ ) ! . (9)
We note that here the phase term, not the emissivity layer term, is involved. The case γ = 3 is detailed in Appendix A and leads to a simplification of the above equations.

Maximum A Posteriori formulation

From the direct model, we can express the observed data as:

dm = P `(x ⊥ , ν)m, B ´+ em , (10)
with m an index which spans the mixed frequency position-on-thesky cube, (x ⊥ , ν)m the corresponding coordinates, B the actual magnetic field and em an error term which accounts for noise and model approximations. Using vector notation, equation (10) simplifies to: d = P(B) + e with d = {dm} the vector collecting all the observations, P(B) = {P `(x ⊥ , ν)m, B ´} and e = {em}. Our inverse problem is to recover the magnetic field vector, B, given some noisy measurements of the polarization, d. Due to the unknown errors in equation ( 10) and to possible strict degeneracies of the direct model, there is not a unique magnetic field that yields a polarization consistent with the observations. We therefore need some means to select a unique solution and, hopefully, the best one given the data.

Probabilities provide a consistent framework to define such a solution; we thus define the sought magnetic field as being the most likely given the observations. It is the one which maximizes the posterior probability:

BMAP = arg max B P(B|d) , (11)
and which is termed as the maximum a posteriori (MAP) solution (see e.g. [START_REF] Pichon | [END_REF]. By Bayes' theorem, P(B|d) = P(d|B) P(B)/P(d), and since P(d) does not depend on the sought parameters B, this amounts to maximizing P(d|B) P(B).

The term P(d|B) is the likelihood of the data given the model, while the term P(B) accounts for any a priori knowledge about the magnetic field. We can anticipate two types of priors: (i) the strict constraint that, to be physically meaningful, the field should be solenoidal: ∇B = 0; (ii) some so-called regularization constraint to overcome the ill-conditioning of the inverse problem and to enforce the unicity of the solution. Without loss of generality, we state that the probabilities writes:

P(d|B) = κ1 exp `-1 2 L(B) ´, (12)
P(B) = ( κ2 exp `-1 2 R(B; µ) ´, if ∇B = 0, 0 otherwise. ( 13
)
where the factors κ1 and κ2 do not depend on B and µ accounts for parameters to tune the regularization. Finally, taking the logprobabilities and discarding constants, the maximum a posteriori magnetic field writes:

BMAP = arg min B,∇B=0 Q(B) , (14)
with:

Q(B) = L(B) + R(B; µ) , (15)
which is the objective function. Before going into the details of the expressions of L(B) and R(B; µ) we can already note that the solution BMAP will depend on the data d and on the regularization parameters µ. The value of µ can be chosen, e.g., to provide the best bias-variance compromise on the sought solution [START_REF] Wahba | Spline models for observational data[END_REF]Golub et al. 2000).

Likelihood

Assuming Gaussian statistics for the noise and model errors, the likelihood of the data is the so-called χ 2 and writes:

L(B) = `d -P(B) ´⊤ • C -1 n • `d -P(B)
´( 16)

with Cn the covariance matrix of the errors. There is a slight issue here because we are dealing with complex values. Since complex numbers are just pairs of reals, complex valued vectors such as d, P(B) and e can be flattened into ordinary real vectors (with doubled size) to use standard linear algebra notation. This is what is assumed in equation ( 16). Under these conventions, the covariance matrix of the errors writes Cn = e • e ⊤ with ⊤ to denote transposition.

Regularization

The regularization term R(B; µ) implements loose constraints to avoid over-fitting the data and enforce local unicity of the solution (see section 4.3). Requiring that the magnetic field be as smooth as possible (while being consistent with the data) matches these requirements and is supported by physics since the magnetic field should have no discontinuities. To simplify further computations, we choose the following particular expression of the regularization R to favor the smoothness of the field:

R = µs ∆ α/4 B 2 ∝ µs X k |k| α | B| 2 , ( 17
)
which scales as the integrated norm of the spatial Laplacian of the field to the power α/4. For a periodic field, this generic smoothing penalty is diagonal in Fourier space. In addition, if the model B is Gaussian and scale invariant, then α may be chosen to be the power law index of the power spectrum | B|2 k of the field. In this case, choosing the specific value of the hyperparameter, µs = 1/| B| 2 k=1 , the MAP solution correspond to the minimal variance Wiener filtered data.

Imposing ∇B = 0

For simplicity, we assume here that the magnetic field is multiperiodic, with period L in all three directions. We may then rewrite the magnetic field as:

B = F -1 • ( B⊥1 e ⊥1 + B⊥2 e ⊥2 ) ≡ Π • B , (18)
where F -1 = F † /N 3 r and F is the forward DFT operator, (e ≡ k/|k|, e ⊥1 , e ⊥2 ) form a spherical basis in Fourier space, while B⊥,i , i=1,2 are the projections over that basis of the Fourier component, B ≡ F • B of the field. Equation ( 18) defines the projector

Π = F -1 • (e ⊥ ⊗ e ⊥ ) • F. Such a field satisfies by construction k • B ≡ 0 , which implies ∇ • B ≡ 0 . (19)
In fact, there is a slight complication at the Nyquist frequencies where only one component of the field is free, see appendix B. Note that the divergence free condition could also be imposed by other means (see e.g. [START_REF] Nocedal | Numerical Optimization[END_REF]. For instance, by adding a quadratic penalty term like P r (∇B) 2 r to the total penalty Q(B). We however found that, in practice, the projector Π led to a better conditioned reconstruction problem.

Implementation

Given equations ( 16) and ( 17) the objective function writes:

Q = (P -d) † • C -1 n • (P -d) + µs ∆ α/4 B 2 . ( 20
)
To minimize Q(B), we used a variable metric limited memory optimization method with BFGS updates (Nocedal 1980) called VMLM and implemented in OptimPack 1 [START_REF] Thiébaut | Optimization issues in blind deconvolution algorithms[END_REF]. Finding the optimal solution, equation ( 14), involves computing the gradient of equation ( 20) with respect to B. Now differentiating equation ( 16) with respect to a magnetic field components we get

∂χ 2 ∂Bi = 2 Re » (P -d) † • C -1 n • ∂P ∂Bi - , (21)
where ∂P/∂Bi for i = x, y, z are given by equations ( 7) and ( 9). Similarly, differentiating equation ( 17) with respect to B yields

∂R ∂Bi = µs F -1 • Bi|k| α . ( 22
)
The VMLM algorithm is a quasi-Newton method which proceeds by solving successive linear problems. Let us therefore first consider in the next section a linearized version of our inverse problem,

1 OptimPack is freely available at http://www-obs.univ-lyon1.fr/labo/perso/eric.thiebaut/optimpack.html.

which may correspond to a physically motivated problem when a good first guess for the magnetic field is known. Note finally that equations ( 7) and ( 9) imply that ∂χ 2 /∂Bi = 0 at Bx = By = 0. Note also that if (Bx, By, Bz) is a solution to equation ( 5), so is (-Bx, -By, Bz). Consequently we expect that the χ 2 will be strongly multivalued as a function of B 2 . The smoothing penalty should in part prevent a pixel-by-pixel flip of the x and y component. It remains nonetheless to be shown that the zero divergence condition is sufficient to avoid flipping the field in regions bound by zeros of these two components, if such regions exist. Addressing these issues will be the topic of another paper.

LINEARIZATION

Let us first consider the situation when a fairly good guess for the overall magnetic field, B0, is known, on the basis, say of a first large scale investigation, or via some modelling of the field as a function of the underlying density (e.g. Cao et al. 2006;Kachelrieß et al. 2007). Let us then seek the departure from this guess. It is then legitimate to assume B = B0 + δB, with, possibly (if the prime guess is accurate enough) δB/|B0| ≪ 1, so that equation (5) becomes:

δP ≡ " ∂P ∂B « B 0 • δB, (23)
where the tensor ∂P/∂Bi is given by its components, equations ( 7), ( 8) and ( 9), while δP ≡ P -P(B0). Now equation ( 23) is likely to be a much better behaved equation as the linearity warrants convexity of the objective function, hence the formal unicity of the solution.

In this paper, we will address two linear problems in turn, one of academic interest, to understand the properties of the inverse problem at hand, while the second one should allow us to carry realistic reconstructions, in the regime when a fair reference field is known. Specifically, we will first assume that the (noise free) data is in the image of (∂P/∂B) B 0 :

d ≡ δPL = (∂P/∂B) B 0 • δB + e , linear problem (I),
while for the second problem (the so called Gauss-Newton approximation) d ≡ δPPL ≡ P -P(B0) + e , pseudo linear problem (II).

We investigate the linear problem in this section and the pseudolinear problem in section 4.

Linear reconstruction

Let us illustrate our method on a problem of realistic scales. This first simulation is carried on a N 3 r grid (Nr = 64) with Nν = 64 frequencies. The reference field B0 is chosen constant and set to 1 µG everywhere for each component, the power spectrum of the perturbation field δB has a power law index α = 2 and its RMS is 0.01. Data are simulated linearly (see section 3) and are noised with a SNR= 20. The y component and the x direction are not plotted since very close to the x component and the y direction. One can see that the z component is better reconstructed than the x or y components which is consistant with the variance measurements and the global conditioning of the problem (see section 3.3). The reconstruction is carried on a Nr = 64 grid with Nν = 64 frequency channel. The data are generated linearly (see section 3) with a SNR= 20. Bottom left: maps of |B| for a transverse section after smoothing of the fields. The green images represents the input field while the superposed white contours show the recovered one. Bottom right: power spectra of the input field (solid line) and the recovered one (crosses). As expected, the recovered power spectrum is damped at higher frequencies because of the regularization. To illustrate this we added the power spectrum of a reconstruction with SNR=200.

the results for the y component and the x direction are similar to the x component and the y direction, they are not plotted. Here, the solid lines stand for the input field and the dashed lines for the recovered one. It is clear that the two fields are very similar and that the z component is the best recovered (see section 3.3). The bottom left panel shows a map of |B| for a transverse section after smoothing. The smoothing is made by convolving the field with a four pixels full-width at half maximum (FWHM) gaussian. The green features represent the input field and the reconstructed one is shown in the superposed white contours. The bottom right panel shows the power spectra of the input field (solid line) and the re-covered one (crosses). Finally, figure 2 represents the field lines of the input field (top) and the recovered one (bottom). These figures show that, if the frequencies are correctly sampled (see section 3.2), the linear inverse problem (I) recovers qualitatively well the underlying field. The local and global properties of the field can be reconstructed provided that the linearization remains valid which will be investigated in section 4.

It is of interest to study the conditioning of the linear problem for two reasons (i) to understand the spatial spectral feature of the solution; in particular the biases of the eigenvectors of the linearized problem which induces anisotropy in the distribution of errors around the solution; (ii) to constrain the best sampling strategy in order to recover B. Eventually it will also have an impact on our ability to carry out the non linear reconstruction.

The requirements to set up a good conditioning of the global inverse problem can be formulated in steps. First a necessary condition is to make a proper choice of the (electromagnetic) frequency sampling, which can be achieved by looking at a smaller subproblem on a given LOS; however, this optimal sampling does not warrant a good global conditioning; we therefore investigate the quality of the global linear reconstruction by looking at different elements of the reconstruction covariance matrix in (spatial) frequency space. In particular, we will show that the quality of the reconstruction is anisotropic and depends on the components of the field, B, which is confirmed by looking at the eigenvectors of the covariance matrix for a low dimensional problem.

Conditioning of a line of sight and frequency sampling

One can see easily that in the relation between polarization and magnetic field (equation ( 5)), each line of sight is independent of the other. The link between them is provided by the solenoidal condition. In this subsection we will not consider this condition and the matrix (∂P/∂B)B 0 becomes block-diagonal. Moreover, the three components can be separated leading to three different matrices, (∂P/∂Bx), (∂P/∂By) and (∂P/∂Bz). The field B0 is taken constant and its modulus set at 1 µG. In this case, all blocks are the same and the study of the conditioning is reduced to the study of three Nν × Nr matrices with Nν the number of frequencies and Nr the number of pixels in the z direction.

Numerical investigations show that the conditioning of (∂P/∂Bx) depends mainly on the ratio K h ne Bz/ν 2 leading to the conclusion that the conditioning is dominated by the exponential term of equation ( 7). It follows that (∂P/∂By) has the same behavior as (∂P/∂Bx) since the exponential terms are the same in both equations ( 7) and ( 8), which is confirmed numerically.

Recall that since in this section the reference field is chosen constant, so is Bz; therefore the best sampling for the frequencies is to have ν -2 n -ν -2 n+1 constant, that is a constant step for the squared wavelength; hence: λ 2 n ≡ λ 2 0 + (n -1) ∆λ 2 with n = 1, . . . , Nν the index of the frequency/wavelength. So that the complex expo-nential becomes

e 2 i K ne Bz h m λ 2 n /c 2 = e 2 i K ne Bz h m λ 2 0 /c 2 e 2 i K ne Bz h m (n-1) ∆λ 2 /c 2 (24)
with m = 1, . . . , Nr the pixel index along the line of sight. The value of ∆λ 2 must be chosen in such a way that the frequency dependent complex exponentials are uniformly sampled on the complex circle. Hence K ne Bz h Nr (n -1) ∆λ 2 /c 2 must be a multiple of π for any n. With L = h Nr the maximum probed depth and taking the smallest multiple, this yields:

∆λ 2 = π c 2 K ne Bz L . ( 25
)
With this particular choice, the matrices (∂P/∂Bx) and (∂P/∂By) take the following form:

(∂P/∂B x/y )n,m = C x/y e Nriβ " λ 2 0 + (n -1)π K neBzL « γ-1 4 × " e -iβ e -2iπ(n-1) Nr « Nr-m , ( 26
)
where β = 2 K ne Bz h λ 2 0 and C x/y is a different constant in the x and y directions. If the factor `λ2 0 + (n -1)π/(K neBzL)

´γ-1 4 is set to 1, the matrix is a unitary Vandermond matrix and its conditioning is 1 (Cordova et al. 1990).

Accounting for this factor impairs the conditioning but it stays close to unity. The elements of the last matrix, (∂P/∂Bz) are just geometrical series of the elements of (∂P/∂Bx). Thus, they read:

(∂P/∂Bz) n,m = Cze Nriβ " λ 2 0 + (n -1)π K neBzL « γ+3 4 × 1 -exp(-iβ(Nr + 1 -m)) exp " -2iπ Nr (n -1)(Nr + 1 -m) « 1 -exp " -iβ exp(- 2iπ Nr (n -1)) « .
where Cz is yet another constant. At this stage, there is only one free parameter left, the first frequency λ0. The conditioning of `∂P/∂B x/y ´being always close to unity, the value of λ0 must be chosen in order to minimize the conditioning of (∂P/∂Bz).

Figure 3 (top panel) represents the conditioning of (∂P/∂Bz) as a function of λ0 for different grid sizes. The curves are very similar in shape and the best conditioning is represented by the red dots. In the bottom panel the wavelength providing the best conditioning for (∂P/∂Bz) is plotted as a function of the grid size. It appears that λ0 ∝ √ Nr and the precision on λ0 is not really important since the minimum of the curves are not really marked. These particular choices of λ0 give a conditioning of 1.29 for `∂P/∂B x/y

´, whatever grid size.

Conditioning of CMAP and a posteriori variances

Let us now investigate the a posteriori variances of different spatial frequencies of the reconstructed field. This covariance matrix can be written as

CMAP = (A T • C -1 n • A + C -1 B ) -1 , (27)
where A ≡ (∂P/∂B0) • Π with Π the projector that cancels the divergence (cf. equation ( 18)) and C -1 noise and the signal respectively 3 . Here we seek ĈMAP, the Fourier transform of CMAP as we want to understand the relative error in the amplitude of the spatial modes of B. Because of the potential high dimensionality of our problem, the covariance matrix, ĈMAP is not computed directly. We chose instead to compute the selected values by solving for B the following equation with a conjugate gradient method (CGM, Shewchuk 1994;[START_REF] Nocedal | Numerical Optimization[END_REF]:

n and C -1 B ≡ µs F -1 diag(|k| α )F
ŴMAP • B = Bref . (28)
Here, ŴMAP = Ĉ-1 MAP and the solution, B, found by the CGM is

B = ĈMAP • Bref . ( 29
)
The reference field, Bref , is equal to 1 or ±i for the chosen k frequency and its opposite -k in order to have a real field, and 0 elsewhere. The elements Bk and B-k of the solution are combinations of the covariance of k and -k and the variance of k. It allows us to determine the a posteriori variance of the chosen spatial frequency

3 Throughout this section (unless stated otherwise) we assume that α is given by minus the powerspectrum index of the sought magnetic field, and choose µs = 1/P(k = 1), which corresponds to the minimum variance solution.

k. To check this method, the same variances were also computed by the iterative VMLM method. One can check that:

( Bin -Bout) • ( Bin -Bout) † = ĈMAP, (30)
where † denotes conjugate transposition, Bin and Bout stand respectively for the input field and the reconstructed one in Fourier space. As expected, the higher the number of iterations, the closer the two estimates of the variance. Figure 4 represent the evolution of the a posteriori variance of different spatial frequencies k for the different components of the field in different directions (along a LOS or transverse to it) as a function of the SNR. The size of the box is Nr = 16 and the number of frequencies is Nν = 16. Figure 5 shows the evolution of the a posteriori variances of the same frequencies as of figure 4, but as a function of the spectral index, α, of the sought field (for a SNR= 20). As expected, the variance decreases as the index increases. In Figure 4 the SNR is defined as

SNR = RMS(data)/σn, (31)
with σ 2 n standing for the noise variance. The results for the By and Bz fields in the x direction are not plotted because there are exactly the same as those in the y direction. First note that the variances, σ 2 k for the Bz component of the field are much smaller in amplitude relative to the other components. For the Bx and By fields, at low SNR, the Wiener prior is important in the reconstruction, explaining the separation of the three curves corresponding to three different scales. In Fourier space, ĈB = µ -1 s diag(|k| -α ) with α the spectral index of the power spectrum of the input field. If the regularization dominates, ĈMAP ∼ |k| -α , which corresponds to the values on the figures when the SNR is low.

For the transverse frequencies (bottom panels), the behaviour of the variances is well understood. At low SNR, the Wiener prior dominates the reconstruction for the Bx and By components but not for the Bz one. Increasing the SNR implies increasing the relative weight of the data compared to the prior. So equation ( 27) becomes

CMAP ∼ (A T • C -1 n • A) -1 , when SNR → ∞ . (32)
If we assume a Gaussian white noise, Cn = σ 2 n I with I the identity matrix, equation (32) becomes

CMAP ∼ σ 2 n (A T • A) -1 , (33)
so CMAP ∝ σ 2 n or given equation ( 31), CMAP ∝ SNR -2 which is the slope of these curves. Finally, note that there is no symmetry breaking between the x and y directions and between the x and y components of the field or between the sine and cosine modes in CMAP. Now, consider the x and y components of the field along a LOS (top panels). At low SNR, the Wiener prior still dominate, providing the same value as in the transverse direction. Then, the variance decreases as SNR -2 but reaches a threshold and stagnate. It is clear on the figures that there is a symmetry breaking between the x and the y components of the field and a separation between the sine and cosine modes. At first it may be surprising that the variances reach a threshold since the frequencies have been chosen to provide the best possible conditioning for (∂P/∂B0) along a LOS (see section 3.2). In fact this is a consequence of the solenoidal condition. Recall that for the global inverse problem, the relevant linear model is A = (∂P/∂B0) • Π, where Π is the projector given by equation (18). This projector changes the matrix (∂P/∂B0) and adds off-diagonal terms to the block diagonal matrix considered in the previous subsection. In effect, the solenoidal condition degrades the global conditioning relative to the one LOS problem (but recall that without it we have an ill posed problem). In turn this changes the eigen structure of ĈMAP and therefore its projection in Fourier space.

Indeed, let us compute directly the whole matrix ĈMAP for a smaller, more tractable Nr = 8 constant reference magnetic field with Nν = 8 frequencies sampled following the procedure defined in section 3.2 4 . Figure 6 shows the global conditioning of the covariance matrix CMAP as a function of the SNR. One can see that the mixing of the LOS has a significant effect on conditioning, even 4 As expected the curves of the variance as a function of the SNR found previously are recovered exactly with this direct calculation. though the frequencies were chosen optimally. Figure 6 also shows that at realistic SNR, the global conditioning remains bounded and could be improved, e.g. for the purpose of numerical convergence, by artificially increasing the hyperparameter µs. Note finally that even though the global conditioning increases with the SNR, the variances all decrease, as expected.

Eigenspace analysis

In order to understand the plateau on figure 4, let us also explicitly diagonalize WMAP for the smaller above-described Nr = 8 problem with a SNR= 20. The corresponding spectrum is plotted on figure 7, bottom right panel. The global conditioning of WMAP is about 10 5 (consistently with what was shown on Figure 6 for CMAP), but note importantly that there is a cluster of eigenvalues followed by a gap. This gap is consistent with the plateau seen on figure 4. When increasing the SNR, one expects to filter out less and less eigen modes, and therefore to access more and more eigenvectors (corresponding to decreasing eigenvalues) in the reconstruction. However, when reaching the gap, although the SNR increases, no more eigenvalues are available for a while. The lower eigenvectors, encoding informations on higher frequencies, are not within reach, and the a posteriori variance of these frequencies stagnate, as seen in figure 4. If the SNR increases further, these eigenvalues (and therefore their associated eigenvectors) will be sampled, and we expect that the σ 2 k variances will decrease again 5 . The modulus of the first eigenvector (associated to the highest eigenvalue) is plotted on the top panels in the x-y (left) and x-z (right) planes. It is clear on these figures than the x and y directions are isotropic while the z one is anisotropic for this eigenvector. Moreover, the component of the power spectra in the bottom left panel show that the Bz component clearly differ from the other two components.

However, all of the main eigenvectors do not behave in the same way. Some of them clearly break the symmetry between the x and y directions or/and between the x and y components leading to the differences in the curves of figure 4. Finally note that the main eigenvectors are fairly high frequencies fields. So, the a posteriori variances be smaller for high frequencies than for low ones, which is reflected by the top panels of figure 4.

5 in other words, the plateau seen in the variance per mode in the top panels reflects the fact that those modes have non zero contributions from the low signal to noise eigen modes (i.e. eigen modes of Log(C G )

C -1/2 A • C B • C -1/2 A with low eigen values, where C -1 A ≡ A T • C -1 n • A).
Figure 6. Global conditioning of the (a posteriori) covariance matrix C MAP as a function of the SNR. The higher the signal to noise, to more difficult the inversion, but the smaller the covariance a posteriori. The 3D matrix, A = (∂P/∂B 0 )•Π appears to be more poorly conditioned than its 1D counterpart even though the sampling in electromagnetic frequency was the same as in section 3.2. It remains bounded and within reach of double precision calculation.

VALIDITY OF THE LINEAR APPROXIMATION

Linear and pseudo linear inversion

Let us first carry out a linear inversion of the same pertubative field δB, with RMS(δB) = 10 -3 µG, while considering both the linear (I) and the pseudo linear (II) data sets (see section 3). We work here on a Nr = 64 grid, with Nν = 64 frequencies, a constant reference field of module 1 µG and SNR=20. Recall that for the linear minimum variance solution, the hyperparameter µs = 1/P(k = 1) (see section 2.2.2), while for the the pseudo linear data set it may be tuned. Figure 8 top panel shows the input z component for the input field (solid line) along a given LOS and the output ones (dotted line for the linear data, δPL and dashed line for the pseudo linear, δPPL) while the bottom panel shows the different power spectra. As previously, the field recovered from linearized data sets fits quite well the input one. The recovered pseudo linear field, though . Top panels: maps of the modulus of the field corresponding to the first eigenvector of W MAP in the x-y (left) and x-z (right) plans for a 8 3 constant reference magnetic field with Nν = 8 frequencies sampled as explained in section 3.2 with a SNR= 20. The first eigenvector appears to be isotropic in x and y and anisotropic in the z direction. Bottom left: power spectra of the three components of this eigenvector. The anisotropy of the z component is clearly visible and in good agreement with the results found in section 3.1 (figure 1) and 3.3 (figure 4). Bottom right: spectrum of the eigenvalues of W MAP . somewhat different from the linear one, remains fairly close to the original field. The corresponding powers pectra are also shown on Figure 8 and confirm that the recovered field in setting (II) is quantitatively redder.

Second order residuals

Let us now study the second order residuals to quantify the domain of validity of the linearization. For this purpose, we subtract to the total polarization its zero and first order expansion to obtain (P-P0 -(∂P/∂B) B 0 ∝ δB 2 ) and we divide this quantity by the first order term (P -P0 ∝ δB). Figure 9 represents the average of this quantity as a function of RMS(δB). Here the perturbation consist of a single frequency and single component field. The solid lines represent the results obtained with a Bx component along the LOS at the lowest mode, while the dashed lines correspond to the lowest transverse mode of the Bz component. The dark curves rep-resent the real part, Q, of the polarization while light ones stand for the imaginary part U (see equation ( 1)). At very low RMS(δB), numerical noise dominate but decreases as the RMS increases. After reaching a minimum, note that the quantity plotted increase as RMS(δB) since ∝ δB 2 /δB and thus ∝ δB. As expected, the lower the RMS(δB), the better the linear approximation and the better the reconstruction. Note also the significant amplitude difference between the Bz and Bx components; we interpret this as a difference between the second derivatives of the field, which in turn, impairs the accuracy of the linearization for the z component. This should not be a limitation when carrying the non linear reconstruction using a method such as VMLM, as the amplitude of the subsequent changes in the magnetic field will be scaled by the inverse second derivatives. . Top: Bz along a LOS for the input field (solid line) and for the recovered fields with linear data, δP L (dashed line) and pseudo linear ones, δP PL (dotted line) (see section 3). Bottom: power spectra of these three fields. Note that the power spectrum of the reconstructed field from the pseudo linear data set is steeper.

Towards the non linear problem

Up to now, we have only considered the situation where B0 was assumed to be constant. What happens to the conditioning when we add spatial frequencies to B0 or/and over ne? It is easy to see that adding transverse frequencies to the x or the y component of B will not change the conditioning of a LOS. Indeed, according to equations ( 7) and ( 26), only the constants C x/y are modified and vary for each LOS, but remain constant along each of them, which has no effect on conditioning. On the contrary, if the modulation is along a LOS, C x/y is no longer constant, and varies for every pixel along a LOS. However, given that the conditioning is dominated by the exponential terms in the Vandermond approximation, it doesn't change dramatically. Hence the choice of λ0 and the sampling frequency remain the same but the conditioning increases slightly; it can reach 3 for `∂P/∂B x/y ´and 40 for (∂P/∂Bz). The situation is a priori more dramatic for the z component of the field or for the electronic density ne. Indeed, the addition of a transverse modulation has significant consequences, as the value of Bz (or/and ne) in equation ( 25) becomes different for each LOS. Therefore, the value of ∆λ 2 should in principle be different for each LOS to conserve the best conditioning. In practice it is simplest to take the average of Bz (or/and ne) as a guess. However the conditioning per LOS increases signicantly and the quality of the reconstruction should be affected.

However, it appears that the global conditioning of CMAP does not change dramatically compared to the constant reference field value, whatever the frequency and the amplitude of the added modulation. The solenoidal condition appears to be very effective. In fact, the repetition of the spectral analysis carried in section 3.4, shows that the main difference will be in the gap seen on figure 7. Adding modulation on a constant field induces earlier, deeper gaps. At fixed SNR, the number of useful eigenvalues for the reconstruction decreases with the modulation. The inversion can still be carried, but will be more biased by the lack of resolved eigenmodes. As a final illustration, figure 10 shows an implementation of the linear inversion on a more realistic reference field, B0 which is extracted from a magneto-hydrodynamical simulation (Kowal & Lazarian 2007), perturbed by a power-law fluctuation with a power spectrum of α = 2 and a relative amplitude of 10 -2 from a virtual data set of SNR=20. Note that for this more realistic illustration the electronic density ne is not constant but extracted from the same simulation. Both the shape of the correction and its power-spectrum are well recovered for this relative amplitude reflecting that although non constant model and electronic density impair the conditionning, reconstructions remain possible.

CONCLUSION AND PERSPECTIVES

We investigated the problem of reconstructing the threedimensional spatial structure of the magnetic field of a given simulated patch of our Galaxy, using multi-frequency polarized maps of . left top panel: map of a slice (of width 0.047kpc) of input reference magnetic field, B 0 ; right top panel: map of the same slice but for the known electronic density ne; left bottom panel: the input reference magnetic field, the input perturbation and the recovered one along a LOS. The perturbation field is a power-law fluctuation with a power spectrum of α = 2 and a relative amplitude of 10 -2 from a virtual data set of SNR=20; right bottom panel: input and recovered power spectra of the perturbation field.

the synchrotron emission at radio wavelengths. When starting from a fair approximation of the magnetic field, we were able to obtain a good estimate of the underlying field by using a linearized version of the inverse problem considered, up to a 64 3 grid size. The spectral analysis of the strictly linear problem (with a constant reference field, and the simulated data obtained through a linearized model) allowed us to specify the best sampling strategy in electromagnetic frequency, and predict a spatially anisotropic distribution of posterior errors. The best sampling strategy is in equal ∆λ 2 ; it follows from the shape of (∂P/∂B0) along one LOS, which can be approximately recast into a unitary Vandermond matrix when this particular sampling is used. The errors on the reconstructed Bx and By components of the field are shown to be larger than the error on the Bz component. This anisotropy can be traced back to the shape of the posterior covariance, and ultimately of the linearized model which is highly anisotropic, as only the z component of the field induces Faraday rotation.

We considered in turn three more realistic cases: (i) a pseudo linear model (linear reconstruction of non-linearly simulated data), (ii) a varying reference model B0, and (iii) a varying reference model B0 and a (known) varying electronic density ne. We found that for these reconstructions, the global conditioning of the minimum variance solution remained tractable. Finally, we investigated the case where the reference field is given by the outcome of a magneto-hydrodynamical simulation, and is perturbed by an additional fluctuating component of known power spectrum. We showed that even in this case the linear reconstruction quality is reasonable. This leads us to claim that a full non-linear reconstruction, based on a Gauss-Newton sequence of linear sub-problems of varying reference field, should be achievable.

Possible extensions of this work, beyond the scope of this paper, involve investigating systematically the degeneracies of the non-linear inversion. It would be worthwhile to construct specific estimators for the (possibly anisotropic) local power spectrum of the field (see e.g. Lazarian & Pogosyan 2006). Finally, from a modelling point of view, one of the main limitations of the present method is that we had to assume known thermal and relativistic electronic densities, in order to obtain a well posed inverse problem from synchrotron emission data alone. However, we could in principle relax this assumption by adding extra data constraining the electronic densities (e.g. Hα data, see Haffner et al. 2003) or emission measures of pulsars, and attempt a joint reconstruction of the magnetic field and the electronic densities. Any prior statistical information (e.g. extracted from MHD simulations) of possible correlation between B and ne could be used in this context. Another possibility would be to use the extra information given by the circular polarization of synchrotron emission (see Appendix C); this circular polarization, if negligible in the case of low energy sources (like our Galaxy), is measurable in the case of relativistic radio sources (see e.g. Jones & Odell 1977), and opens a way to constrain the electronic density together with the magnetic field structure of the source.

Next, consider the frequencies that have no conjugate, i.e. the frequency ki = 0 (constant) and ki = Ny (Nyquist frequency) where the index i represents the Cartesian coordinates. Let us define F1 as the set of these two particular values, i.e. F1 = [0, Ny], and F2 the set of all the other values, i.e. for a vector of dimension N ,

F2 = [-(N/2 -1), -(N/2 -2), ..., -1, 1, ...N/2 -2, N/2 -1].
When the three components of k belong to F1, the reality condition of the field is merely Im B = 0. After putting this imaginary part to 0, the field can be projected into the Cartesian basis. The difficulty arises when one or two components belong to F1. For example, consider the frequency k = (kx, ky, kz) with kx ∈ F1, ky and kz ∈ F2. In this case, condition (iii) become Bk = B-k⋆ where k = (kx, -ky, -kz) is the "opposite" of k. The problem is that in this case, e k⊥1/2 = -e -k⊥1/2 and the above discussed method can no longer apply. Fortunately, the combination of condition (ii) and Bk = B-k⋆ leads to the following set:

k • B ≡ 0, and Bkx = 0. (B3)
So, the trick is to put the faulty component to 0 and to generate the other two as previously but in 2D space. Now, if k2D = (ky, kz), we generate B = B⊥2D e ⊥2D , where e 2D ≡ k2D/|k2D| and e ⊥2D form a polar basis in Fourier space. As previously, the vectors e ⊥2D are chosen in such a way that e k 2D ⊥2D = -e -k 2D ⊥2D . In this 2D representation, conditions (ii) and (iii) lead to:

Bk 2D ⊥2D = -B⋆ -k 2D ⊥2D
, and Bk 2D 2D = 0. (B4)

Here we have only one degree of freedom left, thus, for these frequencies, we must generate one complex field B⊥2D with the desired power spectrum, and then apply equation (B4). When ky or kz belongs to F1, a similar procedure applies.

In the last case, two component belong to F1. For example, k = (kx, ky, kz) with kx ∈ F2, ky and kz ∈ F1. In this case, condition (iii) become Bk = B-k⋆ where k = (-kx, ky, kz) is the "opposite" of k. Again, e k⊥1/2 = -e -k⊥1/2 and the combination of condition (ii) and Bk = B-k⋆ leads to equations (B3). Consequently, the same procedure follows for these frequencies. After inverse Fourier transform, one can check that the field is real, solenoidal and with the right power spectrum up to the Nyquist frequency.

APPENDIX C: CIRCULAR POLARIZATION

Since the rotating term depends on the density field of thermal electrons ne in the medium, we cannot separate, with the Faraday rotation only, ne from Bz. One way to tackle this problem is to pick up the next coupling term of the Stokes parameters in the (optically thin medium, strong rotativity limit) assumption that describes our medium. This next term is a factor of conversion between linear and circular polarization, that can be considered together with the synchrotron emissivity of circular polarization (Jones & Odell 1977). Following the notations of Sazonov (1969), we write the transfer equation of the polarization tensor I αβ as follows:

dI αβ (z) dz = E αβ (z)-i(Tασ(z)δ βτ -δασT * βτ (z))Iστ (z) , (C1) with I αβ = " I + Q U + iV U -iV I -Q «
, and E αβ (z) is an emissivity term. In the assumption of a thin, strongly rotating medium, we can retain only the rotating terms (the Hermitian part) of T αβ . Defining

T = " h q + if q -if -h
« we can show that the transfer equation can be reexpressed in terms of the (Q, V, U ) "vector" as:

d dz 2 4 Q V U 3 5 = 2 4 EQ EV EU 3 5 - 2 4 h f q 3 5 × 2 4 Q V U 3 5 . (C2)
The fact that this differential equation involves multiplication by a non-Abelian group element -in SO(3) -prevents us from writing a formal solution to the equation in terms of exponentials. However, since we are in the end working on a discretized mesh, we can still write a formal solution to the discrete problem in terms of (finite) sums of (finite) rotations products as we will see below. One important point to notice, linked to the tensor nature of equation C1, is the transformation law of these "vectors" under rotation of the coordinate axes in the plane perpendicular to the line of sight. In this respect, the vector (h, q, f ) behaves the same way as the vector (Q, V, U ), i.e. the (Q, U )and (h, f ) subvectors are rotated by 2ψ when the coordinate axes are rotated by ψ. In the case of a homogeneous medium, this allows Sazonov (1969) and Jones & Odell (1977) to choose the coordinate axes used to measure Q and U so that the V Stokes parameter couples only to U (this is achieved when q is set to 0). In this reference frame, the projection of the (constant) magnetic field is aligned with the second coordinate axis.

In the case of a fluctuating magnetic field, such a scheme is not possible anymore, and we need to rotate the coupling coefficients (best expressed in the reference frame given by the local projection of the magnetic field) in a common, constant, reference frame. Thus, in an inhomogeneous medium, the equation C2 in the common reference frame takes the form: where (Q, U, V ) are measured in the common reference frame, and all other quantities are defined in the frame of the local magnetic field. In the applications we will consider in this paper, the rotation coefficients are dominated by the contribution of cold (thermal) electrons of the medium. In this context, (h, f ) take the following form (Sazonov 1969): h = -q 4 e ne B 2 ⊥ 4 π 2 m 3 e c 3 ν 3 , f = q 3 e ne B π m 2 e c 2 ν 2 .

It is interesting to note that both the frequency dependence, and the dependence on the magnetic field are different in the coupling terms. We note that C3 involves the multiplication of the Stokes "vector" by an element of a non-Abelian group (SO(3)), which precludes finding a formal solution to this differential equation. However, the linearity of the equation in the Stokes parameters, allows us to write a formal solution in the discretized case in terms of sums of products of rotations on the source terms. This equation is very similar to the rigid body type equations encountered in mechanics, with the (major) difference that it is linear. For simplicity, we will consider here a first-order discretization of the problem (i.e. we consider the different fields to be piecewise constant). The solution to the homogeneous Stokes transfer equation can be written as:

2 4 Q V U 3 5 (z) = nz -1 Y i=0 exp(-∆zMi) 2 4 Q V U 3 5 (z = 0) ,
where Mi is the skew-symmetric matrix corresponding to the vector (hi cos(2ψi), fi, -hi sin(2ψi)). This discretized solution ensures the exact conservation of the polarization degree in the abc 0000 RAS, MNRAS 000, 000-000

Probing magnetic fields with multi-frequency polarized synchrotron emission 15 sence of internal sources. It corresponds to the simplest possible case of integration of an equation on the SO(3) Lie Group (e.g. Celledoni & Owren 2001). By linearity, we can find the discrete solution to the transfer equation with sources (C3):

2 4 Q V U 3 5 (z) = nz -1 X i=0 nz -1 Y j=i exp(-∆zMj) 2 4
EQi cos(2ψi)

EV i -EQi sin(2ψi) 3 5 , with ψ = arctan(Bx/By)6 . This expression generalizes equation (6) and can be used to infer the Frechet derivatives of the polarization field with respect to the magnetic field, as was done in section 2 from the integral solution. The source terms in the frame attached to the local transverse magnetic field read (Jones & Odell (1977)):

EQ = C nr B γ+1 2 ⊥ ν -γ-1 2 , EV = -D nr B B γ 2 ⊥ ν -γ 2 (C4)
where nr is the distribution of high-energy electrons in the medium, and C and D are constants that depend on the energy distribution of relativistic electrons. Note that in equation (C4), nr is weighted differently in the expressions of EQ and EV , hence we can in principle disentangle B from nr. Here different assumptions can be made, namely assuming either that nr is related to the distribution of thermal electrons ne, or that it is constant, or that it is related to the magnetic field pressure locally (see Beck et al. (2003) for a discussion of the different assumptions). Another possible path is to add external constraints on either ne (coming for instance from Hα observations (see Haffner et al. 2003) or from dispersion measurements of pulsars), or eventually on the relativistic electron distribution nr with diffuse gamma-ray measurements.

Figure 1 .

 1 Figure1. Top: input (solid lines) and recovered (dashed lines) x and z components of the field along a LOS (left) and along the y transverse direction (right). The y component and the x direction are not plotted since very close to the x component and the y direction. One can see that the z component is better reconstructed than the x or y components which is consistant with the variance measurements and the global conditioning of the problem (see section 3.3). The reconstruction is carried on a Nr = 64 grid with Nν = 64 frequency channel. The data are generated linearly (see section 3) with a SNR= 20. Bottom left: maps of |B| for a transverse section after smoothing of the fields. The green images represents the input field while the superposed white contours show the recovered one. Bottom right: power spectra of the input field (solid line) and the recovered one (crosses). As expected, the recovered power spectrum is damped at higher frequencies because of the regularization. To illustrate this we added the power spectrum of a reconstruction with SNR=200.

Figure 2 .

 2 Figure 2. Field lines of the input (left) and the recovered (right) fields for a 64 3 reconstruction with Nν = 64 frequencies. The fields correspond to a reconstruction with a SNR of 200.

Figure 3 .

 3 Figure 3. Top: conditioning, Cz, of (∂P/∂Bz) as a function of λ 0 for different grid sizes. The red dots represent the best conditionings. Bottom: λ 0 giving the best conditioning as a function of the grid size, Nr. It appears that λ 0 ∝ √ Nr.

Figure 4 .

 4 Figure 4. A posteriori variance of different spatial frequencies k ∈ (1, 2, 3) for the different components of the field in different directions (along a LOS or transverse to it) as a function of the SNR. The size of the box is Nr = 16 and the number of frequency is Nν = 16. The top panels correspond to the variation of σ 2 k for three different values of kz while the bottom panels correspond to varying ky. The cosine mode (thick line) and sine mode (dashed line) are both shown. All variances decrease with increasing SNR as expected, although at different rate, see the main text. Note the different amplitude in σ 2 k for the bottom right panel which shows that the Bz component of the field is better recovered compared to the other components. This reflects the anisotropy of the model A which induces anisotropic reconstruction errors.

Figure 5 .

 5 Figure 5. Same as Figure 4 but as a function of the spectral index, α, of δB for a SNR= 20. As expected, the smoother the expected field, the larger α, the smaller the posterior variances.

  Figure7. Top panels: maps of the modulus of the field corresponding to the first eigenvector of W MAP in the x-y (left) and x-z (right) plans for a 8 3 constant reference magnetic field with Nν = 8 frequencies sampled as explained in section 3.2 with a SNR= 20. The first eigenvector appears to be isotropic in x and y and anisotropic in the z direction. Bottom left: power spectra of the three components of this eigenvector. The anisotropy of the z component is clearly visible and in good agreement with the results found in section 3.1 (figure1) and 3.3 (figure4). Bottom right: spectrum of the eigenvalues of W MAP .

  Figure8. Top: Bz along a LOS for the input field (solid line) and for the recovered fields with linear data, δP L (dashed line) and pseudo linear ones, δP PL (dotted line) (see section 3). Bottom: power spectra of these three fields. Note that the power spectrum of the reconstructed field from the pseudo linear data set is steeper.

Figure 9 .

 9 Figure 9. Average second order of the polarization divided by the first order as a function of RMS(δB). Here RMS(δB) is a single component and single mode field. Results are for a the lowest longitudinal mode for the x component (solid lines) and the lowest transverse mode for the z component (dashed lines). Dark curves represent the real part Q of the polarization while light ones are for the imaginary part U (see equation (1)).

  Figure10. left top panel: map of a slice (of width 0.047kpc) of input reference magnetic field, B 0 ; right top panel: map of the same slice but for the known electronic density ne; left bottom panel: the input reference magnetic field, the input perturbation and the recovered one along a LOS. The perturbation field is a power-law fluctuation with a power spectrum of α = 2 and a relative amplitude of 10 -2 from a virtual data set of SNR=20; right bottom panel: input and recovered power spectra of the perturbation field.

c 0000 RAS

c 0000 RAS, MNRAS 000, 000-000

For instance a magnetic loop close to the z axis (where Bx and By ∼ 0) and its mirror image by symmetry along the z axis have the same χ 2 and almost zero gradient. c 0000 RAS, MNRAS 000, 000-000

Beware that this angle corresponds to πψ in the notations of section 2 c 0000 RAS, MNRAS 000, 000-000

Acknowledgments

We thank Jean Heyvaerts, Martin Lemoine and Guy Pelletier for fruitful comments on the early stages of this work. Special thanks to Alex Lazarian for providing us with his interstellar magneto hydrodynamics simulations.

APPENDIX A: THE CASE γ = 3

For γ = 3, equation ( 5) takes a particularly simple expression

while equation ( 8) simplifies to:

Note that for this value of γ the two derivatives with respect to the transverse magnetic field are thus related:

APPENDIX B: SOLENOIDAL FIELDS WITH FIXED POWER SPECTRUM.

The generation of solenoidal (divergence free) fields with fixed power spectra up to the Nyquist frequency is a tricky problem. The field must obey the three following conditions:

(i) fixed power spectrum:

Given conditions (i) and (ii), the field is best generated in Fourier space. Since the field is multi periodic and we may write

where e ≡ k/|k|, e ⊥1 and e ⊥2 form a spherical basis in Fourier space, while B⊥,i , i=1,2 are the projection over that basis of the Fourier componant of the field. The vectors e ⊥1 and e ⊥2 are chosen in such a way that e k⊥1/2 = -e -k⊥1/2 . The spherical basis is direct for k and indirect for -k. In this representation, conditions (ii) and (iii) become, Bk⊥1/2 = -B⋆ -k⊥1/2 , and Bk = 0. (B2) So, the first step is to generate two complex fields B⊥1 and B⊥2 with the sought power spectrum and then apply equation (B2).