Emissivity vs RA
This material started off as a copy of the README in
~/Orion2004/Data/, which contains the programs for this comparison.
It describes my efforts to extract some graphs from the Orion SII and
H-alpha data, so as to compare them with the results of my numerical
simulations of flows from a plane density interface.
Contents:
Observational comparisons for the plane i-front hydro paper
We find an empirical fit to the [SII] densities:
n = n0 * [(R-R1)/(R2-R)]**A
where
R1 = 0.697116291016009 R2 = 2.33801218015265 n0 = 2489.37802355335 A = 0.99765316171648
A = 1 seems to work quite well too.
So, we now have nearly everything that we need to do the comparisons.
Things still to do....
I had the bright idea of trying to use the [S II] densities for the velocity range that coincides with the peak of the H II emission. Some aspects of this are documented in OrionSiiTere. However, this hasn't worked as well as I had hoped - what has happened is that the effective thickness comes out independent of position. As a result, I will probably go back to using the density from the peak of the [S II] lines instead.
Anyway, the velocity range I am using is +10 -> +20 km/s. This is not quite right since the peak Halpha velocity is around 20 km/s near the Trapezium and about 15 km/s out towards the bar (this from Henney & O'Dell 1999).
I am fitting core+power-laws as a function of radius to the emission measure and electron density data in gnuplot:
f(x) = a/(1.0+(x/b)**c)
For the H alpha surface brightness (which is a proxy for the EM, although I haven't worked out the constant of proportionality yet), I get the following parameters:
a = 119541 +/- 355.8 (0.2976%) b = 119.523 +/- 0.5136 (0.4297%) c = 1.66981 +/- 0.007015 (0.4201%)so we have EM ~ r^{-1.66} asymptotically.
For N_e from the +010+020 data, I get these parameters:
a = 3825.79 +/- 15.16 (0.3961%) b = 189.006 +/- 0.9592 (0.5075%) c = 1.84754 +/- 0.0159 (0.8605%)This has a similar sort of asymptotic dependence (-1.84754) but has a wider core. If I set the core radius the same as in the EM case (119 pixels), then I get
a = 5055.47 +/- 6.224 (0.1231%) c = 1.16651 +/- 0.004762 (0.4082%)which shows a shallower dependence on r.
Combining the two suggests that h ~ EM/N_e^2 ~ r^{-1.66981+2.33302=0.66321}, which would be fine if true. Actually, this seems quite a reasonable way to present things. Maybe it would be best not to restrict the lateral scale height to be the same in the two cases though. Also a better functional form might be
ne(x) = a/(1.0+(x/b)**2 ) **(c/2.0) em(x) = a_em/(1.0+(x/b_em)**2 ) **(c_em/2.0)which gives the following fit for the density, N_e:
a = 3829.5 +/- 12.37 (0.323%) b = 149.815 +/- 2.723 (1.818%) c = 1.43189 +/- 0.02789 (1.948%)and the following for the emission measure, EM:
a_em = 116398 +/- 266.1 (0.2286%) b_em = 93.4084 +/- 0.7598 (0.8134%) c_em = 1.38427 +/- 0.009497 (0.6861%)
This looks real good since if I divide em(x) by ne(x)**2, then I get a curve that looks just like the h-curves from the numerical models.
I have also plotted the density from the +010+020 data against the density determined from the -100+100 data and find that the best straight line fit through the origin has a slope of 1.21456. This is consistent with Baldwin's statement about the correction that needs to be made to the [S II] densities to get them in line with the [O II] densities. This is presumably because much of the [S II] emission comes from the partilly ionized region where the electron densities are lower.
So, the only outstanding issues are
In Teresa's data I am using, we have a pixel size of 0.31 arcsec but we have resampled to 0.535 arcsec to match up with Bob and Takao's Hα data. The position of theta1C is x0m = 188.0, y0m = 315.0.
We are using the Hα map but this requires calibration against either the HST image or the 20cm data. Then, we need to convert to emission measure.
OK, so I have converted the 20cm data from surface brightness to radio τ by using some routines I had lying about from the 3d reconstruction project. Then I have remapped this data onto the same grid as is used in the vcubes. This is done in radiotau.f90.
By plotting the corrected halpha data against the radio tau and fitting a straight line, we get a conversion factor of 6.89255e-06.
Hence conversion factor between Halpha line map and emission measure should be
6.89255e-06*3.0865e18/(8.24e-2*8900.0**(-1.35) * (20.0/30.0)**2.1) = 1.29833614515e20Multiplying this by the first coefficient of the fit to the line intensity versus radius gives
1.29833614515e20*116398 = 1.51123730623e25 cm^{-5}
Hence, we get the height on the axis to be
h = EM0/n0^2 = 1.51123730623e25/3829.5**2 = 1.03050123898e18 cm = 0.334 parsec
If we deal with the angle-averaged data as a function of radius from th1C then we get a very large value for h. This is not very realistic since the peak in the EM and density are offset from the star to the West. Furthermore, the point-by-point values for h are much smaller on the West side than on the East side.
In order to show this, I intend to calculate EM and density profiles
integrated over a wide slit, oriented EW. I have already done this by
sticking regions on the fits images with ds9. These are in
{em,ne}cut-wide-EW.dat. Now I will redo it with a fortran program so
I have more control over things and can add in the dispersion of
values for each position along the slit.
The idea is that these will be compared with the profiles that I had calculated for my models viewed at an angle.
Now that I have redone Teresa's SII reductions (see OrionSiiTere), I need to get back to working out the quantities I will need for comparing with the numerical simulations. First, I will try the one-d functional fits again.
Currently in eden-tere.f90 (which no longer uses Tere's
calibration...) I have the sum of 5 vel slices from +004 to +024. The
idea is that these are reasonably representative of the bulk of the
Halpha emission, and so their density should not need any
correction. Of course, the image of the SII emission in these ranges
does not really look like the image of the Halpha emission so we are
still not measuring exactly the same gas.
Rather than use these thin slices, I will make a single slice that does the whole 20km/s range, then I can more easily play around with different smoothing schemes. I'm also using only 2x undersampling instead of 4x.
printf '%s\n' KPNOsiil jw4-tere.csv 4 24 0.5 | src/wisomom_raw
printf '%s\n' KPNOsiis jw4-tere.csv 4 24 0.5 | src/wisomom_raw
ID=+004+024.wisomom-sum
# raw slits
printf '%s\n' {siil,siis,ratio}_${ID}-raw | src/isoratio
printf '%s\n' ratio_${ID}-raw 0.697 2.338 | src/isopgm
# default interpolation
printf '%s\n' siil_${ID}-raw 1 01 | src/wisomom_smooth
printf '%s\n' siis_${ID}-raw 1 01 | ./src/wisomom_smooth
printf '%s\n' {siil,siis,ratio}_${ID}-01 | src/isoratio
printf '%s\n' ratio_${ID}-01 0.697 2.338 | src/isopgm
# 2 arcsec Gaussian - slit weighting by flux
printf '%s\n' siil_${ID}-raw 3 1 2.0 8.0 2 02GF | ./src/wisomom_smooth
printf '%s\n' siis_${ID}-raw 3 1 2.0 8.0 2 02GF | ./src/wisomom_smooth
printf '%s\n' {siil,siis,ratio}_${ID}-02GF | src/isoratio
printf '%s\n' ratio_${ID}-02GF 0.697 2.338 | src/isopgm
# 2 arcsec Gaussian - slit weighting uniform
printf '%s\n' siil_${ID}-raw 3 1 2.0 8.0 1 02GU | ./src/wisomom_smooth
printf '%s\n' siis_${ID}-raw 3 1 2.0 8.0 1 02GU | ./src/wisomom_smooth
printf '%s\n' {siil,siis,ratio}_${ID}-02GU | src/isoratio
printf '%s\n' ratio_${ID}-02GU 0.697 2.338 | src/isopgm
Note that I have made some improvements to wisomom_smooth, which are
detailed in OrionSiiTere
The uniform-weighted maps look a lot better, since they avoid the discontinuous stripes at the KPNO/SPM border that are seen in the natural-weighted version. This problem with the natural weighting needs looking into in more detail sometime - just not now. The 2 arcsec smoothing is actually the σ of the gaussian, so the FWHM is 2*1.66 = 3.32 arcsec.
These really work very well but I am going to have to wait until tomorrow to describe them properly.
I have corrected a lot of problems with the SII calibrations. Some of
them seem to have been due to the fits data files getting out sync
with the table in jw4-tere.gnumeric. This can happen if I forget to
remove the files *-smooth.fits or *-shear-smooth.fits after
changing some variable like fluxfiddle for a particular slit. What
really needs to be done is that the densities need to be checked using
the horizontal slits, but that is something for Tere. Another
interesting thing I found is that I have to use less than the full
extinction correction, which implies that a significant fraction of
the SII emission comes from in front of most of the absorbing
dust. Further discussion in OrionSiiExinct.
So, I now have decent maps at about 3.3 arcsec resolution. I have to redo all the one-d fits.
We now get the following fits (see OrionFlowComparison above) for the density
a = 3635.46 +/- 8.62 (0.2371%) b = 103.988 +/- 1.059 (1.018%) c = 1.17742 +/- 0.01035 (0.8792%)and for the emission measure
a_em = 1.5e+25 +/- 3.308e+22 (0.2206%) b_em = 96.6924 +/- 0.7677 (0.794%) c_em = 1.43271 +/- 0.009844 (0.6871%)
However, I don't think the radial fits are the way to go. Far better to use...
If we take an EW slice at the declination of θ1C, then we should be able to compare it with the tilted views of my models. We can also compare it with the PaperBaldwin91 data. This is where the fun starts.
set term post eps color solid pix = 0.5*1.07*430*1.5e13 # vcube pixel size in cm x0=188.0 # y0=315.0 # star position set key spacing 1.4 set xlabel "RA" set output "em-x.eps" set ylabel "EM" plot [:] [0:] \ 'all-map-tere.dat' u (($1-x0)*pix):(abs($2-y0)<10?$4:1/0), \ 'BaldwinEden.dat' u ($6*430/500):($8) w lp set ylabel "n_e" set output "ne-x.eps" plot [:] [0:] \ 'all-map-tere.dat' u (($1-x0)*pix):(abs($2-y0)<10?$5:1/0), \ 'BaldwinEden.dat' u ($6*430/500):($7) w lp set ylabel "h" set output "height-x.eps" plot [:] [0:8.e18] \ 'all-map-tere.dat' u (($1-x0)* pix ):(abs($2-y0)<10?$4/$5**2:1/0), \ 'BaldwinEden.dat' u ($6*430/500):($8/$7**2) w lpFollowed by
! epstoimg -o ne-x.png ne-x.eps ! epstoimg -o em-x.png em-x.eps ! epstoimg -o height-x.png height-x.epsThese give the following:
There are two problems with these:
Strangely enough, the heights do match up OK because the two factors
cancel. This seems to imply the Teresa never did the SII ratio
calibration correctly, since one of the sources she was supposed to
use was precisely the Baldwin densities. I am going to have to check
this by doing the full velocity range: -100+100.
I resurrected my eden-tere-allv.f90 program (it would really make
more sense if there was just one program but never mind). This writes
out the ratios to all-map-tere-allv.dat, which I am plotting here,
multiplied by 1.17 in red:
In green I show the same data but for the +004+024 range (also
multiplied by 1.17. The blue line shows the Baldwin data (note that
the Baldwin distances have been multiplied by 430/500 to account for
the different distances we are assuming to the nebula). After
scaling our ratios, they follow the Baldwin ones pretty closely.
The y limits of the graph are the theoretical minimum and maximum
ratios - thankfully we are still well within these even after
rescaling.
You can see that the +004+024 ratios are generally higher than the
-100+100 ones, indicating higher densities. We can fit a straight
line through the origin to one density against the other and we get a
proportionality factor of
A = 1.14239 +/- 0.0004797 (0.04199%)which means that the
+004+024 densities are
higher than the -100+100 densities by a factor of 1.14. This is
similar but slightly lower than I found earlier (see the
discussion there). It is also a much smaller correction than that
employed by PaperBaldwin91, who assume that the [SII] density is
only 0.63 of the fully-ionized density (whereas we get 1/1.14 =
0.88). They get this from a constant-pressure Cloudy calculation - it
would be interesting to check this for a more modern Cloudy run.
It may be that Baldwin is right and that our +004+024 densities are
still underrepresenting the true density - its hard to see how this
can be however. The only way to be sure would be to get some cospatial
[OII] data. I should try and dig out the Jones thesis to see if that
has anything usable in it.
Another worry is that the red leak line may be artificially depressing
the density derived from the -100+100 range.
This needs to be sorted out. Our emission measures come from calibrating the extinction-corrected Halpha map by means of the radio map. I need to go through this again to make sure that it is OK. Baldwin's emission measures come from the infrared H 11-3 line. But they also have Halpha, so we could check that against our own values. I could also try calculating EM directly from the radio map (although then we wouldn't be able to cut out the high-velocity emission). Another thing to look at would be the Wen & O'Dell paper and see what EM they have there.
PaperPogge92 measured the [SII] densities using Fabry-Perot imaging. I have taken their Fig.~5 (lower panel) and extracted the data by sticking dots on it with xfig and then using a gnumeric spreadsheet to convert them to a list of N_e versus declination. This is for an EW strip passing through θ1C. I can then plot all three density datasets on the same graph (now with my correction of 1.17 to get our densities to agree with Baldwin's).
set term post eps color solid cubepix = 0.5*1.07 # vcube pixel size in arcsec x0=188.0 # vcube star position y0=315.0 # set key spacing 1.4 set xlabel "RA" # conversion ratio->eden R1 = 0.697116291016009 R2 = 2.33801218015265 n0 = 2489.37802355335 nefunc(R) = n0 * ((R-R1)/(R2-R)) ratscale = 1.17 # correction to our ratios set ylabel "n_e" set output "ne-delta.eps" plot [:] [0:] \ 'all-map-tere-allv.dat' u (($1-x0)*cubepix):(abs($2-y0)<10?nefunc(ratscale*$4):1/0), \ 'BaldwinEden.dat' u ($9):($7) w lp, \ 'Pogge/Pogge-Ne.dat' w lp set output
This is not good at all! There is general agreement with the densities, but the density peak in our data is just not there in Pogge.
It looks like there are just 3 bad slits: ##55,57,58. The only thing
to do seems to be to remove them. I am now reverting to jw4.gnumeric
being the canonical configuration file. So, the slits I am deleting
are:
Note that I found a 4th bad slit as well. There were many slits close by there, so it won't be missed. The KPNO slits on the other hand are in a sparsely populated region so they will be missed, but it can't be helped.
Things now do look quite a bit better:
However, we still have higher densities than Pogge near θ1C. I think we can live with this for the time being though. Teresa can check it again later.
I have also taken a PA150 slice through my ratio data (using
ds9). This is in ratiocut-allv-PA150.dat. The idea is to compare
this with Pogge's data for the bar region.
Tryiing to get some work done on this at the weekend during Torsten's visit.
First, I will recalculate the emission measure from the radio map itself. This will be somewhat higher than the value I obtained from the Hα data because it includes the high-velocity emission too. This also comes out pretty low as compared to the Baldwin values. I have double-checked the formula for working out the emission measure and i can't find any flaw.
Next thing I have done is to calibrate the HST Hα mosaic by
using the image cal3.fits that Bob gave me many years ago for the
proplyds. This is a PC frame that covers the Trapezium stars, and
which has been calibrated in units of photons cm-2
s-1 sr-1 according to O'Dell & Doi 199?.
I did the calibration by taking a slice across from θ1 D to θ1 B for both the full Hα map and for the calibrated PC frame. Then I brought them to a common spatial scale and found the scale factor between them by trial and error:
Although I have calibrated the HST mosaic (f656.fits), it would probably be better to do it with the vcube data. Never mind, I'll carry on with the HST data for now. Whichever I use, it needs to be corrected for extinction. First of all, though I can check it against the Hα values that Baldwin have.
Baldwin list the Hβ surface brightness in units of erg cm-2 s-1 arcsec-2 (Table 1), where relevant conversion factors are
erg/photon @ Halpha : (1/4 - 1/9) * 13.6 * 1.602e-12 = 3.026e-12
arcsec^2/sr = 206265**2 = 4.2545250225e10
=> (phot/cm^2/s/sr) / (erg/cm^2/s/"^2) = 4.2545250225e10/3.026e-12
= 1.40598976289e22
This is resonable since Baldwin's Hβ surface brightnesses are of order 1.e-12, which gives 1.e10 when converted to the O'Dell units.
After scaling the HST and Baldwin Hα surface brightnesses, I find a discrepancy that can be removed by shifting all the the Baldwin positions by 7 arcsec to the East.
plot [:] [0:1.e11] 'f656-EW-slice.dat' u (($1-13.5)*0.1):($2*2.8e7), \
'BaldwinEden.dat' u ($9-7):($10*$12*1.e-14*1.40598976289e22):($13/4) w xe
This adjustment would also make the electron densities agree much better between Pogge and Baldwin, at least on the western side. There would still be some dicrepancies near the brightness peak, but our data would then fall between the Pogge and Baldwin values, which is reassuring.
Next think to do is to use the Baldwin extinctions to correct their
Hα and then convert to emission measure. I have the E(B-V)
column in BaldwinEden.dat. The A(lambda) is 4.71 at Hα. This
gives reasonable results for the emission measure - here I plot the EM
calculated from the extinction-corrected Halpha against the emission
measure in Baldwin's paper (calculated from the H11-3 line).
The conversion factor of 9.342e-15 is (alpha/4pi) for Halpha.
plot [:] [0:] 'BaldwinEden.dat' \ u ($9-7):($10*$12*1.e-14*1.40598976289e22*exp($11*4.71)/9.342e-15/$8):($13/4) w xe,\ 1
There is a slight discrepancy between the 2 since the Ha/Hb ratio is somewhat anomalous compared with the rest of the H lines.
Using the image emdirect.fits, which is calculated directly fro mthe
radio free-free brightness, I have taken an EW slit that starts at the
position of th1C (emdirect-EW-slice.dat) and is the same width as
the Baldwin slit:
plot 'BaldwinEden.dat' u ($9-7):($8*2./3.), \
'emdirect-EW-slice.dat' u ($1*0.5):($2)
I have multiplied the optically-derived emission measure by 2/3 to take account of the back-scattered component in the line surface brightness due to reflection from the molecular cloud. This seems a bit of a large correction, but it does bring the two estimates into agreement.
Now that I have got the emission measure sorted out I can go back to creating the figures I need for the paper.
I have modified the program eden-tere.f90 in order to include the
correction factor of 0.17 for our 6731/6716 ratios to bring our
densities into line with Baldwin's and Pogge's. We use the velocity
range +004+024 to select the density of the fully ionized gas as
well we can. This density is somewhat higher than the density from the
entire line. The image of the density map is written to
nemap-tere.fits
Model scaling behaviour is discussed in OrionModelScaling
To fit the models, I have done the following:
Star-ifront offset, z_h,0, in cm (from [cde]_timeseries.dat):
half-width of EM profile (from [cde]_latstats_P0201.dat) in terms of
z_h,0:
z_p,0
7.6054E+17 8.2626E+17 8.6107E+17
n_p
Half width of EM peak = 3.06e17 cm
Hence, we can caclulate scale factors for the models by dividing the model width (in cm) by the observed width:
Note these factors also need to be multiplied by cos(inc).
EM peak is around 1.3e25 cm^-5. Fitting the simulation EM peaks gives the following scale factors for 15 deg:
Note that these are factors that I divide the model EM by in order to compare with the observations.
Hence scale factors for Q_H are given by EM*length^2
It seems pretty weird that these vary so much and also that the factors are so large, which would give a very low Q_H of about 3.e48 for Model D. Granted that this is only the fraction that is not absorbed by dust but even so.....
From Felli et al 1993 we have an independent estimate of the total ionizing luminosity (from single-dish radio flux), which gives 9.e48 s^{-1}
By adding up the volume emission measure in the maps I get 5.5963E+48 for emdirect, which is based on the 20cm radio brightness, and 3.3765E+48 for emmap, which is based on the Halpha brightness and covers a smaller area.
Doing the same for the simulations (model_vem.f90) I get
where the second number is the fraction of the stellar ionizing photon luminosity.
The qualitative trend is in agreement with the Q_H scale factors above: the more divergent models are more density-bounded, so would require a larger Q_H to give the same emission measure.
On the other hand, the actual values don't seem to be quite consistent: for instance, with Model D, then 1.0769E+49/3.3765E+48 = 3.19 rather than 4.34 and for Model E, 6.0985E+48/3.3765E+48 = 1.81 rather than 1.71 - actually that is pretty close. Hmm, maybe there isn't anything wrong after all.
Felli et al. have the peak EM being 5e6 pc/cm^6 = 1.5e25 cm^{25} with
a 43 arcsec beam. This is consistent with the peak EM from emdirect,
with the emmap values being about 10% lower. The total VEM from
these maps, then implies that only 30-50% of the total recombinations
appear on our maps.
This would tend to favour Model E since that is the only one that has at least half the recombinations missing from the map. On the other hand, Model E doesn't do very well at fitting the electron densities.
Scale factors for densities are sqrt(EM/length)
For viewing inclination of 15 degrees
Rescaled model distances z_h,0
Rescaled model distances z_p,0
Simulation Q_H is 1.24e49, so deduced values are
Although the observational profiles are for a 60 arcsec wide slice, we are taking a single line for the model profiles. It would be better to take a similarly wide slice for them too.
To do this, we will have to do a rotate-and-splat program. This is necessary anyway in order to produce moment maps of the emission lines.
Design of program discussed in FortranWrotsplat
Rotate and splat program now works perfectly. It gives the surface brightness of Hα directly.
Need to write a program to extract a strip for comparison with the observations.
That's all done - problem is that the densities now come out way too low. This is probably because we are not taking into account the fact that S+ is a minor constituent in the H+ gas. In SiiVersusHzero I calculate a functional form for the S+ fraction as a function of the neutral hydrogen fraction.
The Sii fraction calculation has now been incorporated into wrotsplat
and the predictions for the density now come out much better. I also
now take the same velocity range for the model as for the
observations.
I now get the following for the scale factors (unlike in the previous stuff, these are now factors that multiply the length, density, etc.):
Model length density EM (L D^2) Q_H (L^3 D^2) ------------------------------------------------- D 0.5 1.39 0.96605 0.2415125 C 0.43 1.92 1.585152 0.2930946 E 0.54 2.17 2.542806 0.7414822
So, these are not too different from what we had before - in particular, only Model E comes up with a reasonable value of Q_H.
Using the values from OrionFlowComparison above (and Q_H is 1.24e49) I get
Model z_h,0/1.e17 z_p,0/1.e17 Q_H/1.e49 n_p/1.e3 --------------------------------------------------------- D 4.41575 3.8027 0.2994755 7.068011 C 3.875289 3.552918 0.363437304 11.077632 E 4.931766 4.649778 0.919437928 14.115633
I should compare the difference in velocities of [O I] and Halpha since much of the variation in [O I] velocity is due to motion in the molecular gas.
Question is do I use the peak or the mean velocity? The peak is more robust and won't be influenced by high-velocity flows or the back-scattered component. On the other hand, it would be more of a pain to deal wit in the models - I would have to calculate a full velocity cube. I could always use the mean velocity in a restricted range, which would be almost as good and much easier.
In the original Henney & Arthur 1998 paper we got 1-2e49 for the ionizing photon luminosity in the absence of dust. However, the densities in that paper were 50% too high, whih would give more than a factor of 2 in (n^2 r), which should be proportional to the incident flux. On te other hand, is the T is a bit higher than 10^4 in the cusps, then that would put up the densities again.
A better way may be to use one of the well-observed proplyds such as LV2. My 2002 paper has
n0 = 3.0 +- 0.5 e6 r0 = 7.9 +- 0.2 e14
so that EM = omega n0^2 r0 = 8.532e26 cm^-5, where I have assumed omega = 0.12. This gives an incident flux of (2.21832e14/T_4), taking into account the T-dependence of alpha_B.
Projected distance from thC is 7.8 arcsec, giving a physical distance of 5.031e16/sin(i) cm.
Hence we get Q_H = 4 pi D^2 F = 7.0557420228e48 / (sin^2(i) T_4)
The inclination is determined as 50+-10 deg, so sin(i) = 0.766, and T_4=1.2 according to the static Cloudy models.
So, we get a final result of
Q_H = 1.00196722167e49 s^-1
If we take the errors seriously we have 0.025 for r0, 0.17 for n0, and xxx 0.15 for sin(i), so the total error is or 0.6. This seems too large an error really: the observations actually give n0**2 r0 to quite a good precision. 20% is probably more like it.
Dust optical depth is sigma N_H where sigma ~ 5.e-22 and N_H is the column density through the flow. For LV2 we have N_H = 2.37e21, whereas for the nebula as a whole we have about 1 - 2.e21 depending on how I calculate it. So EUV tau will be about unity in both cases.
We can't really use Bob's KPNO data for comparing with the EW trend in the kinematic profiles. This is because he has assumed that the mean velocity for each slit is the same. However, for our OI, SII, and SIII data we have the horizontal slits, so we can check the EW variation in the velocities.
With the [O I] velocities, there is some EW variation. The mean velocity shows a small gradient from around 26 near to th1C to around 21 to the E or W. The peak velocity is around 30 just to the E of th1C, changing to around 25 on the E and W sides. In the dark bay on the far E, it becomes about 20 but the line is very faint.
The [S III] velocities are similar to [O I] on the E side in the Dark Bay but become incresingly blue-shifted towards the W.
This is plotted by the program vewplot.f90 for the same 60 arcsec
strip as we used before. For this we need to find the position of th1C
on the maps. It looks like it is at (x0,y0) = (178,373)
We also need to check the EW trend against the horizontal slits. I have a feeling that the W side of the OI slits might be a bit too blue.
The OI velocities may be partly caused by motions in the molecular cloud, which are described in OrionMolecVelocities. The Wilson et al results have the CO7-6 lines going:
This more or less follows the trend in the [OI] lines except that we don't see any return to the red in the West - this may be because we aren't going far enough over. Also, we have a shift to the blue beyond about 20 arcsec E of the Trap. This is not mentioned in the CO paper. Part of it seems to coincide with the red apron that is seen in the high ionization lines. Part corresponds to the dark bay, which is also bluer in SIII. Interestingly, all these regions are faint in OI so it may be some spurious effect. Looking at it more closely (using the velocity slice images) I see that the faint regions have a blue mean velocity because there is a diffuse blue-shifted component (around +14km/s) that tends to dominate in regions where the main red component is weak. Comparing the peak velocities might turn out better. Also might be better to weight the mean velocities by the intensity.{27} -> {26} -> {25} ... {28} Trap -> Ridge-> West ....
No, the weighting by flux made very little difference, but the peak velocities are definitely better. More discussion of the OI and SIII data reduction in OrionOiSiii
Update 28 Dec 2004 The above is all rubbish because I was using
the wrong config file and hence not subtracting the sky lines (see
OrionOiSiii). I have also now checked the EW trend of the OI line
using the horizontal slit spec341-horiz.fits in ds9 (making a
projection region and using the cursor key to zoom it down the line)
and this reproduces the general trend that the line goes bluer in the
W (then goes back to red again, but this is beyond the range of the
vertical slits).
These 3 plots were made with vewplot.gnuplot and are for the same 60
arcsecond wide slit that we have been always using.
Check against horizontal slits now described in OrionOiSiii
Note that we should also compare our results with previous studies of the Orion kinematics such as Pakonin 1979 and Wilson 1997. Wilson shows a map of the radio recomb line velocities, which bears some resemblance to our SIII map but there are differences too. He gets the same blue tongue coming in above the bar and also gets the blue velocities out to the W. On the other hand he doesn't get the red ridges at the bars so well, probably because of a lack of spatial resolution. He has a much larger field of view than we do and he confirms that things go blue again to the East in the dark bay, reaching a maximum blue just below the position of M43. Then as you go even further E the velocity goes red again, as the brightness really drops away.
They want to explain this by saying it is the stellar wind that is causing the positive velocities near the Trapezium (they work in V_LSR). This is not totally out of the question but our data indicates that the positive velocities are generally associated with bars. To be fair, the red bit to the W of the Trap doesn't seem to be associated with a bar (this is the red apron that can be seen in the velocity cubes), so this could be a wind effect.
In our discussion we need to talk about the multiple axes of symmetry and that the geometry is more complicated than in our models. For instance, our model has the flow axis pointing to the E, whereas on a large scale it is obvious that it turns round and points down to the SW. We also have the flow from the bar, whose axis points NW. Maybe we should make a 3d model using blender or something in order to show these.
Now we have established that the mean velocity works pretty well for OI, then we have an easier job with the synthetic profiles.
The next thing we have to do is calculate the critical density for the
OI and SIII lines. We did something similar for the SII lines before,
so we can try and find that... it was in ~Orion2004/Data/Cloudy/. I
have adapted it for the oi and siii lines and run it for densities
from 1.e1 to 1.e5 for constant-T (1.e4 K), one-zone models. I have
plotted the (emissivity/n_i n_e) against n_e and get a constant +/- a
few %. No sign of any collisional deexcitation. O^0 always has trace
abundances in these models, but that shouldn't matter.
I've looked at the way that the lines are calculated in Cloudy. The
SIII line is done in coolsulf.c line 207 and is done using the
3-level atom routine atom_pop3.c. The 6312 line is the 3-2
transition, with the 1-3 transition corresponding to
3722Å. Hence I should use the 3722 wavelength in my 2-level
calculation. This is an excitation T of 38300 K, so there will be a
very strong T-dependence. The A-values are
2-1: 0.08 s^-1 3-1: 0.81 s^-1 3-2: 2.22 s^-1
and the collision strengths are
C12 = cs14+cs24+cs34 = (1.34+4.03+6.708)* (1.e4) **(-0.06) = 6.95 C13 = cs15+cs25+cs35 = (0.109+0.327+0.545)* (1.e4) **(0.02) = 1.8 C23 = cs45 = 0.501*(1.e4) **(0.11) = 1.38
The statistical weights are 9,5,1
Hence the level 2 will be affected by collisional deexcitation first and this might even result in extra population of the level 3, which maybe why the normalised emissivity rises with density.
The OI 6300Å calculation in cooloxyg.c line 97 et seq is even
more complicated. It even includes optical depth in the lines and
does the escape probability.
What to do now?
The first came out a bit strange so I am leaving it for now (it seems to depend very much on the hardening of the radiation field - see discussion in SiiVersusHzero).
The second is now done. The velocities look pretty good for SIII
(around -10km/s) but are rather high for OI (around -6km/s) Next thing
to do is to do the EW slices of the velocity moments. Will do a new
program vewmodplot.f90.
This is now done. I am now on doing the graph that compares the models with the obserations. The velocities look sort of OK but the model widths are much too small, althoughI haven't yet subtracted the thermal and instrumental widths.
The instrumental width can be estimated from the geocoronal component to the [OI] line. In OrionOiFit we had a width of 8.0 for the 150 micron slit and 6.4 for the 70 micron slit.
a **2 + b **2 = 8.0 **2 (a/2) **2 + b **2 = 6.4 **2 => (3/4)*a **2 = 8.0 **2 - 6.4 **2 => a = sqrt( (4/3)*(8.0 **2 - 6.4 **2) ) = 5.54 b = 5.77
This is the σ in 
so we need
to multiply it by 2*sqrt(log(2)) 1.6651=, giving an instrumental
FWHM of 9.22 km/s. Alternatively, we can look at the comparison lamp
spectra. Doing this implies an instrumental width of about 4-5 pixels,
which is 8-10 km/s for the 150 micron slit.
For the thermal width we have FWHM = 21.4 * sqrt(T_4/A), so
W_oi = sqrt(9.22 **2 + 21.4 **2 / 16) = 10.66 W_siii = sqrt(9.22 **2 + 21.4 **2 / 32) = 9.97
After doing all this, the observed widths are still about 5-10 km/s larger than the model widths. The observed reduced [OI] width is pretty constant at around 14 km/s everywhere except at the bright bar, where it falls to about 9 km/s. [SIII] shows more variation in its linewidth, showing increases near the Trap and to the E and W edges. However, the increase towards the W edge is due to HH202.
I had made a mistake with my wrotsplat runs - I was still using the
restricted velocity range. Now I have corrected it to use the ful
lvelocity range and things work out much better with the mean
velocities. The widths are still far too low, though.
I am working out what it would take to make the Alfvén velocity high enough. It seems that a 5.e-4 G field is needed in the [O I] emission zone.
P_mag = B**2 / (8*pi) v_a **2 = B**2 / (4*pi*rho) => rho * v_a **2 = 2 * P_mag whereas rho * c_i **2 = P_gas => (P_mag/P_gas) = 0.5 * v_a **2 / c_i **2
In our ZH007 Cloudy model we had B = 1.e-4 at the illuminated face, which gave v_a = 2 km/s and P_mag/P_gas = 1%.
Assuming gamma_m = 4/3, then the 3-fold density increase from the ionized peak to the OI zone implies a 2-fold increase in B, so we would need about 2.5e-4 G in the fully ionized bit. This would have v_a = 5 km/s in the fully ionized bit.
Problem with all this is that v_a **2 ~ B **2/rho ~ rho **(gamma_m-1) so v_a ~ rho **0.5*(gamma_m-1) ~ rho 1/6. Hence v_a would increase as the density increased in the neutral/molecular part.
So, two possibilities:
References to history of blisters, bright rims, etc discussed in BlisterFlows
Remind him about project to measure the rotation measure in Orion.
Around line 153 of tfidle.c we have
/* term with hi added June 4, 93, to account for warm pdr */ phycon.edensqte = (float)((dense.eden + dense.xIonDense[ipHYDROGEN][0]/1e4)/phycon.sqrte); phycon.cdsqte = (float)(phycon.edensqte*COLL_CONST);
This is then used in the n-level atom calculations to do the collisional terms. Why has he added 1.e-4 times the neutral H density to the electron density here? Presumably this is to account for neutral H collisions, but why the factor of 10^-4? If it were just due to the difference in velocity of electrons and atoms it would be sqrt(m_e/m_H) = 0.023. But there must also be some difference in the collisional cross-sections between electron-ion and atom-ion....
Then we have another line
dense.EdenHCorr = dense.eden + dense.xIonDense[ipHYDROGEN][0]*1.7e-4;that is similar but not exactly the same. This does have a reference to Drawin (1969) Zs Phys 225, 483. The variable
dense.EdenHCorr is used in
the hydrogenic and helium-like calculations, whereas phycon.cdsqte
seems to be used everywhere else.
According to Steenbock, W.; Holweger, H. 1984A&A...130..319S the Drawin formula is for collisions between identical particles. They give a generalization for collisions with other atoms (but only neutrals I think). Kiselman 2000 (arXiv:astro-ph/0010300) says this is at best an order of magnitude estimate.