

In [1]:

    
import numpy as np
import matplotlib.pyplot as pl
import seaborn as sn
import pymilne
import pylte
import scipy.integrate
%matplotlib inline

Index

Mechanisms for generating asymmetries

Effect of noise in asymmetries

Mechanisms for generating asymmetries

In this section we discuss the minimal mechanisms to produce asymmetries in Stokes V. Two possible asymmetries can be defined for Stokes V. Area asymmetries $\delta A$ are defined as the total area of the Stokes V profile normalized to the total area of the absolute value. It is zero for an antisymmetric profile in which both the positive and negative lobes have the same area: $$ \delta A = s \frac{\int V(\lambda) d\lambda}{\int |V(\lambda)| d\lambda} $$

The amplitude asymmetry $\delta a$ is defined as the amplitude difference between the two lobes normalized to the sum of the amplitudes: $$ \delta a = \frac{a_b-a_r}{a_b+a_r} $$

In a Milne-Eddington atmosphere, the Stokes V profiles are antisymmetric, so that both asymmetries are zero.

Let us test this. We synthesize a spectral line in the Milne-Eddington approximation and compute the area and amplitude asymmetries.



In [30]:

    
lambda0 = 6301.5080
JUp = 2.0
JLow = 2.0
gUp = 1.5
gLow = 1.833
lambdaStart = 6300.8
lambdaStep = 0.015
nLambda = 100

lineInfo = np.asarray([lambda0, JUp, JLow, gUp, gLow, lambdaStart, lambdaStep])

s = pymilne.milne(nLambda, lineInfo)


stokes = np.zeros((4,nLambda))

BField = 100.0
BTheta = 20.0
BChi = 20.0
VMac = 2.0
damping = 0.0
B0 = 0.8
B1 = 0.2
mu = 1.0
VDop = 0.085
kl = 5.0
modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes = s.synth(modelSingle,mu)
wl = np.arange(nLambda)*lambdaStep + lambdaStart



In [31]:

    
deltaA = scipy.integrate.simps(stokes[3,:], x=wl) / scipy.integrate.simps(np.abs(stokes[3,:]), x=wl)
ab = np.max(stokes[3,:])
ar = np.abs(np.min(stokes[3,:]))
deltaa = (ab-ar)/(ab+ar)
print "Area asymmetry = {0}".format(deltaA)
print "Amplitude asymmetry = {0}".format(deltaa)









    



Area asymmetry = -1.58761444822e-08
Amplitude asymmetry = 9.98237377772e-05

Generating amplitude asymmetries

Amplitude asymmetries can be generated easily when several components lie in the resolution element. For instance, imagine we have two different components with exactly the same magnetic field but with a different velocity.



In [85]:

    
BField = 100.0
BTheta = 0.0
BChi = 20.0
VMac = 0.0
damping = 0.0
B0 = 0.8
B1 = 0.2
mu = 1.0
VDop = 0.1
kl = 5.0
modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes1 = s.synth(modelSingle,mu)

BField = 800.0
BTheta = 180.0
BChi = 20.0
VMac = 10.0
damping = 0.0
B0 = 0.8
B1 = 0.2
mu = 1.0
VDop = 0.1
kl = 5.0
modelSingle = np.asarray([BField, BTheta, BChi, VMac, damping, B0, B1, VDop, kl])
stokes2 = s.synth(modelSingle,mu)

stokes = 0.5*(stokes1+stokes2)



In [87]:

    
f, ax = pl.subplots(ncols=2, nrows=2, figsize=(12,10))
ax = ax.flatten()
labels = ['I/I$_c$','Q/I$_c$','U/I$_c$','V/I$_c$']
for i in range(4):
    ax[i].plot(wl, stokes[i,:])
    ax[i].set_xlabel('Wavelength [$\AA$]')
    ax[i].set_ylabel(labels[i])
    ax[i].ticklabel_format(useOffset=False)
pl.tight_layout()
deltaA = scipy.integrate.simps(stokes[3,:], x=wl) / scipy.integrate.simps(np.abs(stokes[3,:]), x=wl)
ab = np.max(stokes[3,:])
ar = np.abs(np.min(stokes[3,:]))
deltaa = (ab-ar)/(ab+ar)
print "Area asymmetry = {0}".format(deltaA)
print "Amplitude asymmetry = {0}%".format(deltaa*100)









    



Area asymmetry = -1.5855565196e-06
Amplitude asymmetry = -4.4614204403%

Generating area asymmetries

In order to generate area asymmetries, one needs to have variationa along the line of sight of the physical properties. The reason is that adding zero-area profiles, one cannot generate non-zero-area profiles.



In [2]:

    
atmos = np.loadtxt('hsra.model', skiprows=2)
lines = np.loadtxt('lines.dat')
nDepth, _ = atmos.shape
atmos[:,5] = 0.0



In [3]:

    
f, ax = pl.subplots(nrows=2, ncols=2, figsize=(10,6))
ax[0,0].plot(atmos[:,0],atmos[:,1])
ax[0,1].plot(atmos[:,0],atmos[:,3])
ax[1,0].plot(atmos[:,0],atmos[:,4])
ax[1,1].plot(atmos[:,0],atmos[:,5])
ax[0,0].set_xlabel('Optical depth')
ax[0,0].set_ylabel('Temperature [K]')
ax[0,1].set_xlabel('Optical depth')
ax[0,1].set_ylabel('Velocity [km/s]')
ax[1,0].set_xlabel('Optical depth')
ax[1,0].set_ylabel('Magnetic field strength [G]')
ax[1,1].set_xlabel('Optical depth')
ax[1,1].set_ylabel('Magnetic field inclination [G]')
pl.tight_layout()



In [4]:

    
wl = np.linspace(6301.0,6302.0,100)
pylte.initAtmos(atmos)
pylte.initLines(lines, wl)



In [5]:

    
out, continuum = pylte.synthLines(atmos)



In [6]:

    
f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
labels = ['I/I$_c$','Q/I$_c$','U/I$_c$','V/I$_c$']
loop = 0
for i in [0,3]:
    ax[loop].plot(wl, out[i,:] / continuum)
    ax[loop].set_xlabel('Wavelength [$\AA$]')
    ax[loop].set_ylabel(labels[i])
    ax[loop].ticklabel_format(useOffset=False)
    loop += 1
pl.tight_layout()



In [7]:

    
transition = 50
atmosNew = np.copy(atmos)
atmosNew[0:transition,3] = -1.0
atmosNew[transition:,2] = 0.0
atmosNew[transition:,4] = 0.0
f, ax = pl.subplots(nrows=2, ncols=2, figsize=(10,6))
ax[0,0].plot(atmosNew[:,0],atmosNew[:,1])
ax[0,1].plot(atmosNew[:,0],atmosNew[:,3])
ax[1,0].plot(atmosNew[:,0],atmosNew[:,4])
ax[1,1].plot(atmosNew[:,0],atmosNew[:,5])
ax[0,0].set_xlabel('Optical depth')
ax[0,0].set_ylabel('Temperature [K]')
ax[0,1].set_xlabel('Optical depth')
ax[0,1].set_ylabel('Velocity [km/s]')
ax[1,0].set_xlabel('Optical depth')
ax[1,0].set_ylabel('Magnetic field strength [G]')
ax[1,1].set_xlabel('Optical depth')
ax[1,1].set_ylabel('Magnetic field inclination [G]')
pl.tight_layout()
outNew, continuumNew = pylte.synthLines(atmosNew)



In [8]:

    
f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
labels = ['I/I$_c$','Q/I$_c$','U/I$_c$','V/I$_c$']
loop = 0
for i in [0,3]:
    ax[loop].plot(wl, out[i,:] / continuum)
    ax[loop].plot(wl, outNew[i,:] / continuum)
    ax[loop].set_xlabel('Wavelength [$\AA$]')
    ax[loop].set_ylabel(labels[i])
    ax[loop].ticklabel_format(useOffset=False)
    loop += 1
pl.tight_layout()
deltaA = scipy.integrate.simps(outNew[3,:], x=wl) / scipy.integrate.simps(np.abs(outNew[3,:]), x=wl)
ab = np.max(outNew[3,:])
ar = np.abs(np.min(outNew[3,:]))
deltaa = (ab-ar)/(ab+ar)
print "Area asymmetry = {0}%".format(deltaA*100)
print "Amplitude asymmetry = {0}%".format(deltaa*100)









    



Area asymmetry = 50.7670554525%
Amplitude asymmetry = 58.5973703925%



In [11]:

    
transition = 50
atmosNew = np.copy(atmos)
atmosNew[0:transition,3] = -3.0
atmosNew[transition:,2] = 0.0
f, ax = pl.subplots(nrows=2, ncols=2, figsize=(10,6))
ax[0,0].plot(atmosNew[:,0],atmosNew[:,1])
ax[0,1].plot(atmosNew[:,0],atmosNew[:,3])
ax[1,0].plot(atmosNew[:,0],atmosNew[:,4])
ax[1,1].plot(atmosNew[:,0],atmosNew[:,5])
ax[0,0].set_xlabel('Optical depth')
ax[0,0].set_ylabel('Temperature [K]')
ax[0,1].set_xlabel('Optical depth')
ax[0,1].set_ylabel('Velocity [km/s]')
ax[1,0].set_xlabel('Optical depth')
ax[1,0].set_ylabel('Magnetic field strength [G]')
ax[1,1].set_xlabel('Optical depth')
ax[1,1].set_ylabel('Magnetic field inclination [G]')
pl.tight_layout()
outNew, continuumNew = pylte.synthLines(atmosNew)



In [12]:

    
f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
labels = ['I/I$_c$','Q/I$_c$','U/I$_c$','V/I$_c$']
loop = 0
for i in [0,3]:
    ax[loop].plot(wl, out[i,:] / continuum)
    ax[loop].plot(wl, outNew[i,:] / continuum)
    ax[loop].set_xlabel('Wavelength [$\AA$]')
    ax[loop].set_ylabel(labels[i])
    ax[loop].ticklabel_format(useOffset=False)
    loop += 1
pl.tight_layout()
deltaA = scipy.integrate.simps(outNew[3,:], x=wl) / scipy.integrate.simps(np.abs(outNew[3,:]), x=wl)
ab = np.max(outNew[3,:])
ar = np.abs(np.min(outNew[3,:]))
deltaa = (ab-ar)/(ab+ar)
print "Area asymmetry = {0}%".format(deltaA*100)
print "Amplitude asymmetry = {0}%".format(deltaa*100)









    



Area asymmetry = 0.4539848856%
Amplitude asymmetry = 12.5993842107%

Effect of noise in asymmetries

Let us assume that we have the previous green profile and we compute the asymmetries using the same methods as before but adding noise.



In [178]:

    
nSigma = 20
sigmaNoise = np.linspace(0.0,0.01,nSigma)
deltaA = np.zeros((nSigma,300))
deltaa = np.zeros((nSigma,300))
for i in range(nSigma):
    for j in range(300):
        stokesV = outNew[3,:] / continuum + sigmaNoise[i] * np.random.randn(len(wl))
        deltaA[i,j] = scipy.integrate.simps(stokesV, x=wl) / scipy.integrate.simps(np.abs(stokesV), x=wl)
        ab = np.max(stokesV)
        ar = np.abs(np.min(stokesV))
        deltaa[i,j] = (ab-ar)/(ab+ar)



In [189]:

    
cmap = sn.color_palette()
f, ax = pl.subplots(ncols=2, nrows=1, figsize=(12,6))
ax[0].fill_between(sigmaNoise, np.mean(deltaA, axis=1)-np.std(deltaA, axis=1), np.mean(deltaA, axis=1)+np.std(deltaA, axis=1))
ax[0].plot(sigmaNoise, np.mean(deltaA, axis=1), color=cmap[1])
ax[0].set_xlabel('Noise standard deviation')
ax[0].set_ylabel('Area asymmetry')
ax[1].fill_between(sigmaNoise, np.mean(deltaa, axis=1)-np.std(deltaa, axis=1), np.mean(deltaa, axis=1)+np.std(deltaa, axis=1))
ax[1].plot(sigmaNoise, np.mean(deltaa, axis=1), color=cmap[1])
ax[1].set_xlabel('Noise standard deviation')
ax[1].set_ylabel('Amplitude asymmetry')
pl.tight_layout()
pl.savefig('asymmetries_noise.png')

The explanation for the area asymmetry comes from the bias introduced by the computation of the sum of the absolute value of positive and negative numbers.



In [199]:

    
f, ax = pl.subplots()
x = 0.0 + 0.2*np.random.randn(2000,50)
ax.hist(np.sum(x,axis=1), bins=50)
ax.hist(np.sum(np.abs(x),axis=1), bins=50)
pl.savefig('sumAbs.png')

The explanation for the amplitude asymmetry comes from the bias introduced by computing the peaks with a max/min operation.



In [205]:

    
f, ax = pl.subplots()
for i in [4,20,50]:
    x = 1.0 + 0.2*np.random.randn(2000,i)
    ax.hist(np.max(x,axis=1), bins=50)
pl.savefig('maxStat.png')



In [ ]: