NoisyQUInclination

  |   Source

There is a recent debate about which is the distribution of magnetic field inclinations in the quiet Sun. Some people say that the distribution is close to isotropic and others sustain that the amount of what they call "horizontal" fields is around 5 times larger than the "vertical" fields. I show in this calculation which is the effect of the presence of noise on the inferred value of the inclination if the maximum likelihood solution is used. We recently demonstrated that the maximum likelihood value of the inclination in the weak-field regime is biased towards high values of the inclination and here I show some simple calculations.

Let's first start with a function to return the Stokes parameters in the weak-field regime. In this regime, we find

$$ V(\lambda) = -C\lambda_0^2 g_\mathrm{eff} B_\parallel \frac{\partial I(\lambda)}{\partial \lambda} \\ Q(\lambda) = -\frac{C^2}{4} \lambda_0^4 G_\mathrm{eff} B_\perp^2 \cos 2\phi \frac{\partial^2 I(\lambda)}{\partial \lambda^2} \\ U(\lambda) = -\frac{C^2}{4} \lambda_0^4 G_\mathrm{eff} B_\perp^2 \sin 2\phi \frac{\partial^2 I(\lambda)}{\partial \lambda^2} \\ $$

where $C=4.67 \times 10^{-13}$ G$^{-1}$ $\unicode{x212B}^{-1}$

If we define the following merit function $$ \chi^2 = \sum_i \frac{V(\lambda_i)-V^\mathrm{obs}(\lambda_i)}{\sigma_i^2} + \sum_i \frac{Q(\lambda_i)-Q^\mathrm{obs}(\lambda_i)}{\sigma_i^2} + \sum_i \frac{U(\lambda_i)-U^\mathrm{obs}(\lambda_i)}{\sigma_i^2} $$

Following Martínez González et al. (2012), one can find the maximum likelihood solution for a set of observed Stokes parameters by using:

$$ B_\parallel = -\frac{1}{C} \frac{\sum_i V(\lambda_i) I'(\lambda_i)}{\sum_i I'(\lambda_i)^2} \\ B_\perp^2 = \frac{1}{C^2} \frac{\sqrt{(\sum_i Q(\lambda_i) I''(\lambda_i))^2+(\sum_i U(\lambda_i) I''(\lambda_i))^2}}{\sum_i I''(\lambda_i)^2} \\ \tan \theta = \frac{B_\perp}{B_\parallel} $$
In [3]:
import numpy as np
import matplotlib.pyplot as pl
import pdb
%matplotlib inline

def stokesSynth(B, inclination, azimuth, stokesIp, stokesIpp):
	cte = 4.67e-13	
	Q = - cte**2 * (B * np.sin(inclination))**2 * np.cos(2*azimuth) * stokesIpp
	U = - cte**2 * (B * np.sin(inclination))**2 * np.sin(2*azimuth) * stokesIpp
	V = - cte * B * np.cos(inclination) * stokesIp
	
	return Q, U, V

Let us assume that we have a field of 200 G that is inclined 20$^\circ$ with respect to the line-of-sight. We assume a simple Gaussian spectral line with 40% absorption for the Fe I line at 630.2 nm.

In [24]:
np.random.seed(1234)

nPix = 300
B = 200.0
inclination = 20.0 * np.pi/180.0
azimuth = 10.0 * np.pi/180.0
sigma = 0.2
f = 0.4
cte = 4.67e-13
lambda0 = 6302.5

Lambda = 6302.5**2
gbar = 2.5
nLambda = 100

wavelength = np.linspace(-1.5,1.5,nLambda)
stokesI = 1.0 - f*np.exp(-wavelength**2 / sigma**2)
stokesIp = (2*f*wavelength / sigma**2 * np.exp(-wavelength**2 / sigma**2)) * Lambda * gbar
stokesIpp = (2*f*(sigma**2-2.0*wavelength**2) / sigma**4 * np.exp(-wavelength**2 / sigma**2)) * 0.25 * Lambda**2 * gbar**2

Q, U, V = stokesSynth(B, inclination, azimuth, stokesIp, stokesIpp)

fig, ax = pl.subplots(nrows=2, ncols=2, figsize=(8,8))
ax = ax.flatten()
ax[0].plot(wavelength, stokesI)
ax[1].plot(wavelength, Q)
ax[2].plot(wavelength, U)
ax[3].plot(wavelength, V)
pl.tight_layout()

Now we do the experiment of using the maximum-likelihood value of $B_\parallel$ and $B_\perp$ from Martínez González et al. (2012) when we add more and more profiles.

In [25]:
pl.close('all')
fig, ax = pl.subplots(nrows=3, ncols=4, figsize=(12,6))
	
noise = [1e-2,1e-3,1e-4,1e-5]
for loop in range(4):
    BparML = np.zeros(nPix)
    BperpML = np.zeros(nPix)
    thetaML = np.zeros(nPix)
    QObs = np.zeros(nPix)
    UObs = np.zeros(nPix)
    VObs = np.zeros(nPix)
    for i in range(nPix):
        QNoise = Q[:,np.newaxis] + noise[loop]*np.random.randn(nLambda,i+1)
        UNoise = U[:,np.newaxis] + noise[loop]*np.random.randn(nLambda,i+1)
        VNoise = V[:,np.newaxis] + noise[loop]*np.random.randn(nLambda,i+1)
        QObs = np.mean(QNoise,axis=1)
        UObs = np.mean(UNoise,axis=1)
        VObs = np.mean(VNoise,axis=1)
        
        BparML[i] = -1.0 / cte * np.sum(VObs*stokesIp) / np.sum(stokesIp**2)
        BperpML[i] = np.sqrt(1.0 / cte**2 * np.sqrt(np.sum(QObs*stokesIpp)**2 + np.sum(UObs*stokesIpp)**2) / np.sum(stokesIpp**2))
        thetaML[i] = np.arctan2(BperpML[i], BparML[i])
	
    ax[0,loop].plot(BparML, 'b-')    
    ax[0,loop].set_title(r'$\sigma=${0}'.format(str(noise[loop])))
    ax[0,loop].set_ylabel(r'B$_\perp$')    
    ax[0,loop].axhline(B*np.cos(inclination), color='r')
    ax[1,loop].plot(BperpML, 'r-')
    ax[1,loop].axhline(B*np.sin(inclination), color='b')
    ax[1,loop].set_ylabel(r'B$_\parallel$')
    ax[2,loop].plot(thetaML * 180.0 / np.pi)
    ax[2,loop].axhline(inclination * 180.0 / np.pi)
    ax[2,loop].set_ylabel('Inclination [deg]')
	
pl.tight_layout()

The results show that the line-of-sight component of the magnetic field is quite easy to get right and no bias is seeing, independently of the noise level. The only effect of the noise is in the uncertainty in the estimation of the field. However, this is not the case for the perpendicular component, which converges slowly to the correct value. Only when the noise is small, we find this component quite reliably. If the noise is large ($>10^{-3}$), we find that even adding 300 profiles leaves us far from the correct value of the perpendicular component.

Now, let us assume how an isotropic distribution of fields is inferred with this technique. We generate a set of Stokes profiles from an isotropic distribution and perturb it with noise. Applying the maximum-likelihood estimation of the inclination, we then plot the histograms of the original and inferred inclinations.

In [30]:
nPix = 1000
mu = 2.0*np.random.rand(nPix)-1.0
inclination = np.arccos(mu)[:,np.newaxis]
Q, U, V = stokesSynth(B, inclination, azimuth, stokesIp, stokesIpp)
In [39]:
noise = [1e-2,1e-3,1e-4,1e-5]
pl.close('all')
fig, ax = pl.subplots(nrows=3, ncols=4, figsize=(15,8))
for loop in range(4):
    VObs = V + noise[loop] * np.random.randn(nPix,nLambda)
    QObs = Q + noise[loop] * np.random.randn(nPix,nLambda)
    UObs = U + noise[loop] * np.random.randn(nPix,nLambda)
    
    BparML = -1.0 / cte * np.sum(VObs*stokesIp[np.newaxis,:],axis=1) / np.sum(stokesIp[np.newaxis,:]**2,axis=1)
    BperpML = np.sqrt(1.0 / cte**2 * np.sqrt(np.sum(QObs*stokesIpp[np.newaxis,:],axis=1)**2 + 
                                             np.sum(UObs*stokesIpp[np.newaxis,:],axis=1)**2) / np.sum(stokesIpp[np.newaxis,:]**2,axis=1))
    thetaML = np.arctan2(BperpML, BparML)
        
    ax[0,loop].plot(thetaML * 180.0 / np.pi, 'b.')
    ax[0,loop].plot(inclination * 180.0 / np.pi, 'g.')
    ax[0,loop].set_ylabel(r'$\theta$ [deg]')
    ax[0,loop].set_title(r'$\sigma=${0}'.format(str(noise[loop])))
    ax[1,loop].hist(inclination * 180.0 / np.pi, alpha=0.5, bins=30)
    ax[1,loop].set_xlabel(r'$\theta$ [deg]')
    ax[1,loop].set_ylabel('Original')
    ax[1,loop].set_xlim((0,180))
    ax[2,loop].hist(thetaML * 180.0 / np.pi, alpha=0.5, bins=30)    
    ax[2,loop].set_xlabel(r'$\theta$ [deg]')
    ax[2,loop].set_ylabel('Inferred')
    ax[2,loop].set_xlim((0,180))

    pl.tight_layout()

So, when the profiles are noisy, there is a tendency to give horizontal profiles. It is important to note that the results strongly depend on the specific width of the spectral line.

Another interesting experiment to carry out is the following. We generate many Stokes profiles from an isotropic distribution of fields and then start averaging the profiles in chunks of 4, 8, 10, 20, 50 and 100. If we infer the inclination from them and plot the histogram, we see that we rapidly tend to overestimate the inclined fields, even if the noise is very small.

In [61]:
def rebin(a, shape):
    sh = shape[0],a.shape[0]//shape[0],shape[1],a.shape[1]//shape[1]
    return a.reshape(sh).mean(-1).mean(1)

def histogramNormalized(x, nBins):
    hist, bins = np.histogram(x, bins = nBins)
    maxVal = max(hist)
    hist = [ float(n)/maxVal for n in hist]
    center = (bins[:-1]+bins[1:])/2
    width = 0.7*(bins[1]-bins[0])
    return center, hist, width    

nPix = 10000
mu = 2.0*np.random.rand(nPix)-1.0
inclination = np.arccos(mu)[:,np.newaxis]
Q, U, V = stokesSynth(B, inclination, azimuth, stokesIp, stokesIpp)
noise = [1e-2,1e-3,1e-4,1e-5]
pl.close('all')
fig, ax = pl.subplots(nrows=7, ncols=4, figsize=(18,13))
factors = [1,4,8,10,20,50,100]
for loop in range(4):
    VObsOrig = V + noise[loop] * np.random.randn(nPix,nLambda)
    QObsOrig = Q + noise[loop] * np.random.randn(nPix,nLambda)
    UObsOrig = U + noise[loop] * np.random.randn(nPix,nLambda)
    
    for i in range(7):
        VObs = rebin(VObsOrig, (nPix / factors[i], nLambda))
        QObs = rebin(QObsOrig, (nPix / factors[i], nLambda))
        UObs = rebin(UObsOrig, (nPix / factors[i], nLambda))
    
        BparML = -1.0 / cte * np.sum(VObs*stokesIp[np.newaxis,:],axis=1) / np.sum(stokesIp[np.newaxis,:]**2,axis=1)
        BperpML = np.sqrt(1.0 / cte**2 * np.sqrt(np.sum(QObs*stokesIpp[np.newaxis,:],axis=1)**2 + 
                                             np.sum(UObs*stokesIpp[np.newaxis,:],axis=1)**2) / np.sum(stokesIpp[np.newaxis,:]**2,axis=1))
        thetaML = np.arctan2(BperpML, BparML)
        
        center, hist, width = histogramNormalized(thetaML * 180.0 / np.pi, 30)        
        ax[i,loop].bar(center, hist, align = 'center', width = width)        
        ax[i,loop].set_xlabel(r'$\theta$ [deg]')
        ax[i,loop].set_ylabel('Inferred')
        ax[i,loop].set_xlim((0,180))
        center, hist, width = histogramNormalized(inclination * 180.0 / np.pi, 30) 
        ax[i,loop].bar(center, hist, align = 'center', width=width, alpha=0.3, color='g')
        if (loop == 0):
            ax[i,0].annotate('{0}'.format(factors[i]),(5,0.8))
    ax[0,loop].set_title(r'$\sigma=${0}'.format(str(noise[loop])))
    

pl.tight_layout()
In [ ]: