may dip
Jenny O'Hara: Research log for Summer 2011
Solar Physics REU at Montana State University, Bozeman, MT


Project Summary.

I will focus on understanding energy release in solar flares magnetic reconnection, which refers to magnetic field lines changing connectivities and thus bringing the magnetic field configuration to a lower energy state. The major approach will be analyzing satellite observations (such as SDO observations) of a solar flare. From the analysis, we will derive how much magnetic flux is reconnected, a key parameter in reconnection physics, and we will also investigate the temporal and spatial distribution of reconnection events during a flare, and use that information to infer how much energy is released during reconnection to heat flare loops, which then radiates prominently at EUV and X-ray temperatures.
Dates:01-03 June
For my introductory idl project I was given a UV image of a solar flare and a magnetogram of the flare region together with their coordinate information, which includes the physical size of a pixel, the field of view, and the centre coordinates of the images.
I then tried to figure out how to co-align these two images(doing math by hand). In IDL I tvscl the UV image and then contour the co-aligned magnetogram on top of the UV image.

The two images to begin with were:
idl pic

Using the number of arcsecs per pixal I scaled the magnetogram to have the same number or arcsecs per pixal as the UV image.
idl pic

I then used the relative distances of the two images from the sun's centre, in arcsecs, to coallign them. However, this did not quite work correctly as I used negatives in my code to translate them as they did not completely overlap which produced a slight error.
idl pic

To fix this I cut of the section of the contour magnetogram not overlapping the UV image and produced this better co-allignment:
idl pic

With the images coaligned I then went on to find pixels that are located in positive and negative magnetic field,respectively. Then Sum up all the brightness (counts) in positive region and those in negative region.

I summed the counts for more than 10%, or 30%, or 50% of the maximum counts in the image.

The code I used for this was:

new=congrid(, 101*1.98564/0.5 , 101*1.98654/0.5)
ndata=data(0: -129.1648+400, 19.76: 19.76+400)
newnew=new(129 :400, 0 :400)
ndata=data(0: -129.1648+400, 19.76: 19.76+400)
tvscl, alog(ndata>1.)
contour,newnew, pos=[0, 19.76, 271.83552, 420.76], /noerase, xst=1, yst=1, /dev, lev=[-100,100]
for k=0,2 do begin $
s2=where(newnew lt 0 and data gt ((((2*k)+1)/10.)*max(data))) & $
s1=where(newnew gt 0 and data gt ((((2*k)+1)/10.)*max(data))) & $
print, total(data(s2)), total(data(s1)) & $

The results I obtained were:
6.40346e+06    6.00581e+06
1.12959e+06    772647.
286980.        116040.

And in each case you can see that the total counts in the two kernels are not equal.
I have been told I will be using a more advanced program that follows a similar idea to this in my actual research project so it was good to get a feel for how it worked.

Date: 06 June
Today I began work on my actual project. To begin, I am looking at a flare from March 7 2011 which was an M1.7 classification . I used data of images already co-aligned with each other and which are also already 'congrided' to a smaller size.
I will be using UV observations to identify pixels that are brightned during the flare, and use this information to derive reconnection flux and reconnection rate.
I summed data to derive the uv light curve (scaled by factor of 1e13 in graph). The UV lightcurve gives the number read from the SDO (Solar Dynamics Observatory) satellite CCDs (which is related to photons of a certain waveband received by the satellite detector). This is plotted in blue.
I then compared this with GOES soft X-rays plots which shows the radiation flux at two wavebands at 1AU from the sun. (In W/m^2) I have plotted both of these against time in hours (UT).

UV lightcurve and goes

From this we can see the correlation, which supports the principle that magnetic flux is conserved along the reconnecting field line, so that energy released from reconnection is transported down the fieldlines to the lower atmosphere and gives off radiation signatures.

The 'GOES'data also allows me to plot the tempr (in MK) for the time period of the flare.


Date: 07 June

I also made a movie of 325 images to show the flare.

This movie is also on this link.

Date: 08 June

I have now used the program which produces the total flux reconnection and reconnection rates for the positive magnetic field and then for the negative. It uses a system of cuts in which it takes values greater than some cut*background.This is to identify which pixals are flaring pixals. I then plotted the total reconnection flux (in Maxwells) against time (in minutes) shown in the first graph for both positive (in white) and negative (in colour) magnetic field.
Below I have also plotted the reconnection rates (positive and negative) in Maxwell per second

flux graphs
Neither rates of flux or the totals of flux are equal for positive and negative but appear to follow same trends.

I then tried to scale all the graphs onto one:

bright blue - UV light curve - scaled by 1e10
white and red - GOES radiation - scaled by 1e22 and 1e23 respectively
turquoise -reconnection rate for neg B
light green - reconnection rate for pos B

all on one scaled

Here you can see correlation, though the GOES curves which measure radiation have maximum slightly later. My mentor Jiong was initially concerned with the dip in the rate of reconnection, that it may have been artificial, but now seeing the comparison with the GOES curves and UV lightcurve which both also show a dip at the same time it is more likely that it is a real feature and we are seeing two stages of reconnection/energy release in the flare.

Date: 13-15 June

I have spent the last few days going over a paper, called 'Highly Efficient modeling of dynamic coronal loops' by Klimchuk and Patsourakos (2008), that will be relevent to my project. The paper descibes a 0D model called EBTEL, (Enthalpy based thermal evolution of loops). This model describes the evolution of the average temperature, pressure and density along a strand in a coronal loop. Where each loop is not treated as one object but as hundred of smaller strands to be modelled individually. As this is a 0D model and uses averages (Has 1 value of temp, pressure and density at any given time in simulation) it is able to calculate the strands evolution quickly while still showing good agreement with more complicated 1D models such as Cargill, but takes 4 orders of magnitude less time.

I have also been using the output from the previous programs I was using to look at pixels that are recognized as flaring pixels in the UV images. We assume that these pixels reflect where the flare lower atmosphere is heated, and therefore are where coronal loops are rooted.

To begin with I have been looking at the pixels classified as flaring pixals for cut 6 (the lowest cut) so where they are more than 6*background. I have found 1611 brightened pixals.

I can also show a movie of the pixels that are brightened with time for this cut:

This movie is also on this link (hopefully without being cut off).

I now want to make a UV light curve of each of these brightened pixels. Assuming that each pixal can represent the base of a single strand where a coronal loop is made up of many strands being formed and heated individually and at different times. I will then use these individual lightcurves to construct the heating functions of each single 'strand' of a loop which is one of the inputs of the model in the paper. This will then hopefully allow me to compute coronal plasma evolution.

I first needed to find out how many pixels are brightened and then there coordinates so I can plot each pixals uv readings over time.

I have written a code that tells me the number of brightened pixels and then uses this to create an array of their coordinates, and then similarly to the total UV light curve I made last week I can then plot UV against time for each individual pixel.

An example of some of the light curves for 6 pixels is:

UV lightcurve and goes

This are the uv light curves for pixals:

(295.000, 238.000)
(296.000, 236.000)
(296.000, 237.000)
(297.000, 235.000)
(297.000, 236.000)

Date: 20-21 June

I have been looking at a program Jiong gave me to alter for the flare I'm looking at, in which I construct the heating function, which is an important input into the EBTEL model, from the brightened pixals. The heating function also uses different methods to model the rise time including gaussaian, linear (triangular) and eslope which we can then compare to see which models the data best.(still using data from cut6.)

This is a plot of the total of the heating function using Gaussian method:

heating functions

This is a plot of the total of the heating function using Linear method:

heating functions

This is a plot of the total of the heating function using eslope method:

heating functions

To compare all three heating functions can be seen together:

heating functions

I have tried to plot the heating function against the GOES plots but my timing scale seems a bit off?

heating functions

Date: 24 June

Jiong and I noticed some discrepencies in the heating function program and tried to alter the constraint of the earliest time that it will recognise the earlier peak and alter the times so it can better interpolate and I then plotted the new result. But still the heating function is not able to recognise the first peak in the UV counts, this maybe because it was from the same pixal as the other peak and the program only takes the maximum?

heating functions

The gaussian seems very similar and there is only a slight change on the downward slope. The looking at the gaussian fit of individual pixals there are still a similar number of discrepencies where the gaussian fit is to the right of the uv.

heating functions

This is one example of the gauss fit for one pixal that is off to the right still after modifications.

heating functions

The heating function is still not picking up first peak?
Date: 27 June

Jiong noticed my problem. I had not alterred the heating function in the rise phase but only in the decay phase. The heating function should now be able to model the ealier maximum we see as well as the later one (still using data from cut6).

These are the new individual plots of the total of the heating function using Gaussian, linear, and eslope method from ammended heating function (units is counts per second per pixel)

heating functions

heating functions

heating functions

I have now plotted two of the most useful heating function forms (gaussian and linear) with UV total counts light curve and Goes lightcurves.

heating functions

Below I have plotted two histograms of the rise time of the 'useful' pixals for gaussian and linear outputs of heating function. Where the 'useful' pixals are pixals where the rise time is greater than 0 and the peak is relatively close to what we expect.

heating functions

heating functions

Date: 28 June

I have now tried to realign the histograms and have plotted them together for better comparison.

The total number of 'useful' pixals for linear heating function is 1439
And the total number of 'useful' pixals for Gaussian heating function is 982

heating functions

I will now explain a bit more about the EBTEL model. The idea of which is to equate an enthalpy flux of evaporating or condensing plasma with any excesses or deficit in heat flux relative to transition region radiation loss rate. Where loops are modelled as individual flux tubes or 'strands' where the plasma is assumed to be approximately uniform in cross section. Then after somes derivation and integration over the length of the strands, the model becomes:
heating functions

heating functions

heating functions


n: electron number density
T: Temperature
P: Total pressure
v: Bulk velocity
F: Heat flux
Q: volumetric heating rate

From these see energy enters system only through coronal heating and leaves system only through radiation. (Thermal conduction and flows transport energy between corona and transition region but do not add or remove energy from system)

The real radiative loss rate from the corona would be the integral of the density squared and the radiative loss function but as we are using the averages this can be approximated by: heating functions

Another key assumption in the model is that the radiative loss rates of the transition region and the corona maintain a fixed ratio at all times, C1. The paper uses 4.0 and finds this works well for their examples but for our flare we have been using 2.0, our flare is slightly less powerful.

heating functions heating functions heating functions

Another form of the model also goes on to take into account the energy going into non thermal electron beam but we have disregarded this for now.

Using this model in a program which outputs the radiative coronal loss for each pixal at each time using interpolation.

To begin I have been using a constant strand length for each strand at a pixal which I calculated by taking the average distance between each two footprints of opposite polarity and then modelling the loops as semi circles. I am using the average strand length to be 5.8 10^9 cms.

Using the initial output of the program with this model I have plotted the calculated goes flux with the observed goes data.

In this program I have had to set initial temperature, density which set the default values to be 1mK (10^6 K) and 0.1 (10^9/cm^3) respectively unless otherwise specified.

heating functions

Date: 29 June

This looks like a very good initial output from the comparison with goes but could be improved even more by adjusting the scaling function of the heating function scl0 or tq from 1.5 our initial guess to 2.0 in the program. Result seen below

heating functions

For the 982 useful pixals specified previously I also wrote a short code to find the maximum Temperature and Density for each strand and have plotted these in 2 histograms:

heating functions

heating functions

Similarly, now we have some initial results to go on we can continue to refine the program including adjusting the loop length. So far we have used a crude approximation using constant length. However, we know the length should increase so for the time being we can assume the growth is linear and estimate the min and max loops lengths so we can then interpolate the loop length during the flare.

I have tried to write a program to do this and when I again compare with the goes, I obtain the following plot:

heating functions

Date: 05 July

I am now also going to try to compare the calculated output of the ebtel model with the aia (Atmospheric Imaging Assembly on SDO) data. Here is a plot of the recoreded AIA data at the 171 angstrom wavelength (that focuses on the quiet corona and the upper transition region) compared with what ebtel model computed it to be:

171 aia comp

This does not appear to show very good comparison but 171 picks up reasonably cool temperatures 1 mk and so the early readings may be from previously brightened footpoints from heated lower atmosphere.

I also made a movie of the loops in 171 wavelength to see if it can explain this discrepency.

This movie is also on this link.

Here's a better movie movie

Date: 06 July

I am also going to look at the comparison with the aia data at the 94 wavelength:

heating functions

The shape and timing of this looks good and I just need to explore the magnitude.

Can also see movie of 94 wavelength on this link.

To see more than one wavelength with goes comparisons for diferrent scaling factors in model:

scl0 = total heating scaler(tq) = 2.0
scl1 = non thermal scaler(nth) = 0.0
scl2 = transition region scaler (tr) = 2.0
t0 =1 mk (10^6K)
d0 =0.1 (10^9/cm^3)
heating functions

Date: 07-08 July

I then edited the program to enable me to run for loops to allow me to compute the ebtel outputs and plot the results as shown above for different scaler values to get the best fit with observed data (but keeping scl1=0, non thermal as zero) .
recall defenitions in model: heating functions heating functions

I will post a few of the results here.

for scl0=1.5, scl2=1.5
scl0=1.5 scl2=1.5 scl0=1.5 scl2=1.5
Changes to the value of scl2 (ratio of transition region radiation to coronal radiation) has a prominent influence on the computed aia flux at 171 wavelength. See this UV flux is very sensitive to this intial scaler. The graphs also show large radiation at the start that is not seen in the computed values. From looking at the video of the time we can infer that this is picking up the heated footprints that were heated before the heat from reconnection but the later stage should fit better.

for scl0=1.5, scl2=3.0
scl0=1.5 scl2=3.0 scl0=1.5 scl2=1.5

This is probably one of the best fits to the observed data overall.
for scl0=1.5, scl2=4.5
scl0=1.5 scl2=4.5 scl0=1.5 scl2=4.5

Similarly, taking a larger value of the scl2 shows how sensitive the 171 light curve is.

I have also plotted the histograms of the computed peak temperature and density for the flaring pixals for these different values of the scaler scl2.
From the equations of the ebtel model, smaller scl2 value also means smaller overall radiation and is imput into heating function so smaller overall energy. Therefore, I would expect a smaller value of scl2 (first graph) to produce a smaller peak density and smaller peak temperature. And this is what we observe from the histograms.
Date: 11 July

I have also plotted the different aia 171 computed results for a random individual pixal (number 300) to show how the sensitivity of the scaling factor between the transition and coronal radiation can affect the computational results.

ind aia

I will also show the affect of the scaling factor for the other wavelengths 94 and 335 as investigated above.

ind aia

ind aia

The two rises shown in aia 335 are very interesting an we are unsure whether this is because of a strange combination of the temperature and density or because the response function is sensitive to more than one temperature (has two peaks),

Response function at this wavelength looks like:
aia resp

I will also show the affect of the scaling factor on the two computed 'goes flux' for an individual pixal, so what we imagine as an individual loop.

ind aia

We can also look at this loop's individual density and temperature plots over time for these different scaler values to see how the evolution of the plasma in the individual loop changes for the different ratio of coronal to transition radiation

ind aia ind aia

Date: 13 July

Looking at above 3 examples and doing a number of repeated trials for different scaler values, I believe the best match parameters to be scl0=1.8 and scl2=2.5. These values produce:

new goes

Now that I have experimented with some of the scaling factors my next objective would be to go back and redo the observed goes data. This is because Courtney, who was looking at the same flare as me, discovered that very unusually there were two flares happening on the sun at one time. Although the aia observed data can give readings from only a specific area goes takes its data from much larger area and so she is now working on seperating the data from the two flares. I will then hopefully get a goes light curve just from the flare I am analysing which will show a better fit.

I have now used courtney's analysis of the two flares we know were happening at the same time, where she split the observed counts into ratios between the two flares, as shown in her graph:


I then used this infomation to try and adjust the observed goes data I have been trying to fit my model to so that it better represents the one flare we are focusing on.(still using scl0=1.8 ans scl2=2.5 in model)

new goes

Using this adjusted observed goes data in the comparison of total flux that I have been displaying gives:

new goes

Date: 18-19 July

I am now going to look at a different aspect of my project. Taking the original data from my flare I wrote a code to calculate and save the magnetic flux for each individual 'useful' flaring pixal which we ascertained in our heating function model. I then used the heating function for each of these pixals and integrated to find the area beneath then which would correspond to the total energy for that pixal.

I am continuing to think of each pixal as an individual strand which is one of hundreds within a loop.

We can then calculate the theoretical current in each of these strands as energy/mag flux, where the current is therefore a measure of the free energy in the strand before reconnection.




See similarity in distribution of magnitude of the current for the positive and negative flux tubes, with fairly similar peaks, as we would expect.

Date: 20 July

I also wanted to look at the distribution of the calculated current for strands with different times when they reach their maximum heating.
I have split the total observed pixals into 20 min sections.

5 times of current

This appears to show the flux tubes are most enegetic (more enegy per unit flux released) around the peak (the 2nd or 3rd time interval) as opposed to the rise phase.
In the very late phase we also see very low current which would support the general idea that these later, higher loops are not as twisted as others.

In simple terms we may think of the current as a measure of how twisted the field lines were, but it is in fact more specifically how much energy is released per unit flux as we don't know for sure they relax all the way back to potential (i.e becoming completely untwisted).

Date: 21 July

I am now looking back at what I have done so far and am trying to adjust certain parts to make things more accurate.

I have been trying to obtain new data for aia 171 curve with which I compare my calculated values. As previously there was a high reading before the flare from preexisting aterial or previously heated footprints. I have tried to take this out of the data by only taking data from a smaller area around pur specific flare to enable us to try and produce a better comparison with our calclated values.

For our original scaler values scl0=1.8 and scl2=2.5:

The observed 171 data shown in the top right below now has removed some of the earlier peak which we felt was likely to be showing earlier heated loops and this supports this. We believe that the large peak in the middle that the model does not really show is from the bright footpoint radiation, so would not be in our model of the coronal radiation.

I will now try to pin down the best parameters again to match the decay phase of this curve.

scl0=1.9, scl2=2.0

scl0=2.0, scl2=1.8

However, the actual 171 and computed aia data do show correlation in there rise times and I will go on to try and see visually whether we can see the change at these times.

The bright footpoints can be seen at around 14.30

Date: 26 July

To be accurate I am going to look at the results without background heating. We are particularly interested in the effect of the scl2 on the magnitude and even more the time of the last peak of aia171 in our model. One thought is that the scl2 factor may depend on the temperature, so i have also plotted histograms of the peak temperatures of the strands with there peak temperatures corresponding to the three time intervals of the 3 obvious peaks.

scl0=1.9, scl2=1.8 with no background heating
t t t

scl0=1.9, scl2=2.0 with no background heating
t t t

scl0=1.9, scl2=2.2 with no background heating
t t t

scl0=1.9, scl2=2.5 with no background heating

t t

recall defenitions in model: heating functions heating functions

We can also look at what exactly is going on when we see these different peaks in our observations and model for 171 by looking at the images in 171.

at time 14.16: we see mainly footpoint radiation.
at time 14.30 we see combination of loops and footpoints
at time 15.10 is mainly loops left
at time 15.35 larger post flare loops
t t t t