- Derive velocity structure functions for H II regions of different sizes:
- Lagoon, Carina, Tarantula
- Compare
- It would be good to have a figure that
- Shows images of each H II region in our sample
- Shows the field used for calculating structure function
- Maybe even shows some velocity maps
- Orion
- Lagoon
- Carina
- 30 Dor
- Hub V
- Hub X
- NGC 595
- NGC 604
- I found that this is given by
\[
L = \frac{ 1} { m \, (ln 2)1/m} \, Γ\left(\frac{1}{m}\right) \, r_0
\]
m gamma(1/m) 1 / m (ln2)1/m L / r_0 1.44 / m Err 1/2 1.00 4.16 4.16 2.88 -0.31 2/3 0.89 2.60 2.31 2.16 -0.06 1 1 1.44 1.44 1.44 0.00 1.1 1.06 1.27 1.35 1.31 -0.03 1.2 1.13 1.13 1.28 1.20 -0.06 1.3 1.20 1.02 1.22 1.11 -0.09 5/3 1.49 0.75 1.12 0.86 -0.23 2 1.77 0.60 1.06 0.72 -0.32 - So L/r_0 is of order 2, but varies with m (see above table), falling from 2.3 (m = 0.667) to 1.3 (m = 1.2).
- We could combine [O III] and Ha to get better s/n and fewer gaps
| -1 | 7.9 | 5.9 | -.1 | -3.1 | 2.5 | ||||||||||||||||||
| -.6 | 2.4 | 5.3 | 4.7 | -3.1 | 2.6 | ||||||||||||||||||
| -6.3 | -4 | -2.1 | -.2 | -1.6 | -1.3 | 1.9 | 2.6 | 3.9 | 1.4 | ||||||||||||||
| 3.4 | -2.1 | -5.8 | -3 | -1.2 | 3.3 | -1.7 | -.4 | 1.6 | 2.3 | 4.6 | |||||||||||||
| -.7 | 1.0 | -1.1 | 0 | -1.4 | .1 | -.4 | 1.9 | .8 | 1.2 | 2.4 | 6.9 | ||||||||||||
| -2.8 | -1 | -2.1 | -.3 | -.6 | .2 | 1.8 | -1.9 | -0.5 | 1.0 | .4 | -6.1 | 6.4 | 5.9 | ||||||||||
| -1 | -.9 | -.7 | .3 | -.9 | -.3 | -2.7 | -2.5 | -2.5 | -2.2 | -3.2 | -3.8 | -.6 | 2.8 | 7.1 | 7.3 | ||||||||
| -.3 | 0.0 | -0.8 | 0.9 | .8 | -1.3 | -3.3 | -2.9 | -3.6 | -2.7 | -5.4 | -6 | -5.1 | -3.4 | -5.1 | |||||||||
| -1.5 | -3.9 | 1.5 | -.3 | -2.8 | -1.4 | -1.2 | -2.8 | -1.6 | -3.8 | -0.4 | -2.9 | -4.6 | -4.1 | -.3 | -2.1 | -1.1 | |||||||
| -1.1 | -1.1 | -1.9 | .6 | -.3 | -.2 | -2.9 | -.3 | -4.9 | -2.3 | -.6 | -1.9 | -3.8 | -1.8 | -2.6 | -5.1 | -4.2 | -5.3 | -2.8 | -8.1 | .9 | |||
| -5.6 | -.7 | -4 | -1.4 | -1.5 | 1.4 | .8 | -3.8 | -.5 | -4.7 | -1.6 | -4.8 | -2 | -1.8 | -2.1 | -3.5 | -3 | .6 | 2.9 | |||||
| 1.1 | .2 | .1 | -.5 | -1.2 | -4.3 | -2.2 | -1.7 | -5.5 | -3.1 | -4.8 | -1.2 | -5.3 | -2 | -.8 | -.7 | -3.2 | -2.9 | -.6 | |||||
| .2 | -.9 | -.3 | -2.1 | -4.5 | -4.1 | -4.7 | -3 | -.6 | -2.6 | -2.9 | -2.7 | -2.1 | -1.9 | -1.1 | -3.5 | -5 | -.1 | .2 | 2.9 | -1.1 | |||
| -1.3 | -3 | -1.3 | -3.7 | -5.1 | -2.4 | -2.2 | -4.5 | -3.6 | -3.6 | -1.1 | -.1 | -.8 | -2.4 | -3.5 | -2.7 | .8 | 5.9 | .9 | |||||
| .3 | -.9 | 1.5 | -1.1 | -4.1 | -5.7 | -3.6 | -3.4 | -2.8 | -4.9 | -2.3 | -1.6 | -2.4 | -2.4 | -.8 | -2.3 | -1.9 | 4.8 | -6.1 | -3.1 | ||||
| 2 | .6 | -3.9 | -2.4 | 2.9 | -3.2 | -3.2 | -1.6 | -2.9 | -2.5 | -1.1 | -3.8 | -1.8 | -2.4 | -.5 | -2.1 | -.1 | |||||||
| -3.3 | -5.1 | 2.5 | -2.2 | 2.3 | -2.2 | -3.6 | -4.1 | 2.4 | -2.9 | -1 | -1.4 | 0 | -.9 | -2.6 | 9.9 | -.1 | -2.4 | ||||||
| .2 | .1 | .9 | -2.1 | -3.1 | .2 | -4.1 | -2.1 | -.8 | -.3 | -.9 | 1.7 | 1.2 | |||||||||||
| 2.9 | -1.1 | -2.4 | .4 | .6 | -2.6 | .9 | -.1 | 8.9 | |||||||||||||||
| -.1 | -1.1 | -3.1 | -5 | -5.5 | -2.5 | -.2 | 2.1 | ||||||||||||||||
| 2.9 | -2.1 | -.1 | .9 | 2.9 | |||||||||||||||||||
| -5.1 | -4.1 | -3.1 | -.1 |
import numpy as np
from pathlib import Path
from astropy.io import fits
from astropy.wcs import WCS
ny, nx = len(TAB), len(TAB[0])
data = np.empty((ny, nx))
for j, row in enumerate(TAB):
for i, cell in enumerate(row):
if cell:
data[ny - j - 1, i] = float(cell)
else:
data[ny - j - 1, i] = np.nan
w = WCS(naxis=2)
w.wcs.cdelt = -1/60, 1/60
w.wcs.crval = 83.8185491, -5.3898078
w.wcs.crpix = 11, 12
w.wcs.ctype = "RA---TAN", "DEC--TAN"
outfile = Path("data/Orion") / f"orion-haenel-{ID}.fits"
fits.PrimaryHDU(data=data, header=w.to_header()).writeto(outfile, overwrite=True)- https://sites.google.com/view/phangs/home/data
- Closest galaxies are at 5 Mpc, which is 5x further away than NGC 595+604 in M33, and 10x further away than Hubble V and X
- So we can only look at very large regions
- They have gas and star kinematics, so we could remove the rotation curve by taking the difference
- Also, there are ALMA observations of CO that we could use
- Headlight cloud in NGC 628 has L(Ha) = 6.3e39 erg/s
- Herrera:2020i
- E(Ha) = h c / 6563 ang = 3.02673387238e-12 erg
- => Q(H) = 6.3e39 / 3e-12 = 2e51 phot/s, so similar to 30 Dor
- Also has R = 140 pc and σ = 50 km/s, supposedly
- D = 9.6 Mpc, so 1 arcsec = 46 pc => R = 3 arcsec
- That is going to make it very hard to measure the structure function since seeing = 0.7 arcsec, so only 16 independent pixels
- Once we have found the correlation length, L_0, then we can look at the large-scale (> L_0) POS variation of the small-scale (< L_0) fluctuations.
- In particular variation in
- Magnitude of fluctuations: σ(< L_0)
- Slope of structure function
- And even in L_0 itself
- Although for this, we need to be working on scales biggest than the largest local L_0
- This should be easy to do with the pandas table method of doing the structure function, since we already have all the information - we just need to partition it up
- This is an old idea that I had - explained in ~/Dropbox/Org/garrelt-simulations.org
- Might be fun to apply to structure functions. Could also be used on brightness fluctuations
- I wonder if there is prior art on something similar - I imagine that the Planck people have looked at anisotropy of density structures, as compared with magnetic field
- Yes, I have found Planck-Collaboration:2016a
- Planck intermediate results. XXXII. The relative orientation between the magnetic field and structures traced by interstellar dust
- https://ui.adsabs.harvard.edu/abs/2016A%26A…586A.135P/abstract
- Section 4 explains what they did - some sort of Hessian method
- I still think there might be room for the Stokes parameters approach
- They also mention “anisotropic wavelet techniques”, which might be interesting
- Yes, I have found Planck-Collaboration:2016a
- We could call it polfluctuarization or something!
- [2021-01-26 Tue] This is the first thing I am going to do in order to reactivate the project
- The plan is to use it for the Javier paper, which will compare structure functions and other statistics for a range of H II regions with different sizes and luminosities
- [X] Fix my multi-threaded strucfunc implementation
- This is in
- Now it works again, and faster than ever!
- Multi-threaded version on my laptop is now twice as fast as it ever was on iMac at work
- Time to calculate structure function for 400x400 array is 27.2 s
import sys
from pathlib import Path
import json
import numpy as np
from astropy.io import fits
sys.path.append("../muse-strucfunc")
import strucfunc
try:
LINEID = sys.argv[1]
except:
LINEID = "ha"
fitsfilename = {
"ha": "GAUS_Ha6562.8_060_Will.fits",
"nii": "GAUS_NII6583.45_060_Will.fits",
}
datadir = Path("data/Tarantula/MUSE_R136toWill")
hdulist = fits.open(datadir / fitsfilename[LINEID])
n = None
sb = hdulist[1].data[:n, :n].astype(np.float64)
vv = hdulist[2].data[:n, :n].astype(np.float64)
ss = hdulist[3].data[:n, :n].astype(np.float64)
# Replace spurious values in the arrays
m = ~np.isfinite(sb*vv*ss) | (sb < 0.0)
if LINEID == "nii":
# Remove bad patch from the [N II] map
m = m | (sb > 6e4)
sb[m] = 0.0
vv[m] = np.nanmean(vv)
sb /= sb.max()
rslt = strucfunc.strucfunc_numba_parallel(vv, wmap=sb, dlogr=0.05)
good = (~m) & (sb > 0.001)
rslt["Unweighted mean velocity"] = np.mean(vv[good])
rslt["Unweighted sigma^2"] = np.var(vv[good])
v0w = rslt["Weighted mean velocity"] = np.average(vv, weights=sb)
rslt["Weighted sigma^2"] = np.average((vv - v0w)**2, weights=sb)
class MyEncoder(json.JSONEncoder):
def default(self, obj):
if isinstance(obj, np.integer):
return int(obj)
elif isinstance(obj, np.floating):
return float(obj)
elif isinstance(obj, np.ndarray):
return obj.tolist()
else:
return super(MyEncoder, self).default(obj)
jsonfilename = f"tarantula-strucfunc-{LINEID}.json"
with open(jsonfilename, "w") as f:
json.dump(rslt, fp=f, indent=3, cls=MyEncoder)
print(jsonfilename, end="")time python tarantula-strucfunc.pyResults are written to:
import json
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
def bfunc(r, r0, sig2, m):
"Theoretical structure function"
C = 1.0 / (1.0 + (r/r0)**m)
return 2.0*sig2*(1 - C)
data = json.load(open("tarantula-strucfunc-ha.json"))
sns.set_color_codes()
fig, ax = plt.subplots(figsize=(5, 5))
figfile = "tarantula-strucfunc-plot-ha.pdf"
pixscale = 0.2 # arcsec
pixscale *= 0.242 # parsec
r = pixscale * 10**np.array(data["log10 r"])
B = np.array(data["Unweighted B(r)"])
sig2 = data["Unweighted sigma^2"]
B_w = np.array(data["Weighted B(r)"])
sig2_w = data["Weighted sigma^2"]
# Plot fit to unweighted strucfunc
rgrid = pixscale * np.logspace(0.0, 2.7)
r0 = np.interp(sig2, B, r)
m = 1.22
flabel = rf"$m = {m:.2f}$, $r_0 = {r0:.1f}$ pc, $\sigma^2 = {sig2:.0f}$ (km/s)$^2$"
ax.fill_between(
rgrid,
bfunc(rgrid, r0, sig2, m - 0.1),
bfunc(rgrid, r0, sig2, m + 0.1),
color="k", alpha=0.1,
)
ax.plot(rgrid, bfunc(rgrid, r0, sig2, m), color="k", label=flabel)
# Plot points from unweighted strucfunc
ax.plot(r, B, 'o', label="Unweighted")
# Plot fit to weighted strucfunc
r0_w = np.interp(sig2_w, B_w, r)
m_w = 1.30
flabel_w = rf"$m = {m_w:.2f}$, $r_0 = {r0_w:.1f}$ pc, $\sigma^2 = {sig2_w:.0f}$ (km/s)$^2$"
ax.fill_between(
rgrid,
bfunc(rgrid, r0_w, sig2_w, m_w - 0.1),
bfunc(rgrid, r0_w, sig2_w, m_w + 0.1),
color="k", alpha=0.05,
)
ax.plot(rgrid, bfunc(rgrid, r0_w, sig2_w, m_w), lw=0.5, color="k", alpha=0.5, label=flabel_w)
# Plot points from weighted strucfunc
ax.plot(r, B_w, 'o', ms=3, alpha=0.5, label="Flux-weighted")
melnick_r = np.array([2.5, 7.5, 12.5, 17.5, 22.5, 27.5])
melnick_B = np.array([2.0, 2.1, 2.2, 2.2, 2.25, 2.25]) * 18.2**2
ax.plot(melnick_r, melnick_B, 's', label="Melnick+ (2020)", color="y", zorder=-10)
ax.axhline(sig2, color="k", ls="--")
ax.axhline(sig2_w, color="r", ls=":")
ax.legend(title=r"30 Doradus H$\alpha$")
ax.set(
xscale = "log",
yscale = "log",
ylim = [0.5, 1500],
xlabel = "Separation, pc",
ylabel = r"$B(r)$, (km/s)$^2$",
)
fig.tight_layout()
sns.despine()
fig.savefig(figfile)
fig.savefig(figfile.replace(".pdf", ".jpg"))Same but for the [N II] line
import json
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
def bfunc(r, r0, sig2, m):
"Theoretical structure function"
C = 1.0 / (1.0 + (r/r0)**m)
return 2.0*sig2*(1 - C)
data = json.load(open("tarantula-strucfunc-nii.json"))
sns.set_color_codes()
fig, ax = plt.subplots(figsize=(5, 5))
figfile = "tarantula-strucfunc-plot-nii.pdf"
pixscale = 0.2 # arcsec
pixscale *= 0.242 # parsec
r = pixscale * 10**np.array(data["log10 r"])
B = np.array(data["Unweighted B(r)"])
sig2 = data["Unweighted sigma^2"]
B_w = np.array(data["Weighted B(r)"])
sig2_w = data["Weighted sigma^2"]
rgrid = pixscale * np.logspace(0.0, 2.7)
# Plot fit to unweighted strucfunc
r0 = np.interp(sig2, B, r)
m = 0.95
flabel = rf"$m = {m:.2f}$, $r_0 = {r0:.1f}$ pc, $\sigma^2 = {sig2:.0f}$ (km/s)$^2$"
ax.fill_between(
rgrid,
bfunc(rgrid, r0, sig2, m - 0.1),
bfunc(rgrid, r0, sig2, m + 0.1),
color="k", alpha=0.1,
)
ax.plot(rgrid, bfunc(rgrid, r0, sig2, m), color="k", label=flabel)
# Plot points from unweighted strucfunc
ax.plot(r, B, 'o', label="Unweighted")
# Plot fit to weighted strucfunc
r0_w = np.interp(sig2_w, B_w, r)
m_w = 1.05
flabel_w = rf"$m = {m_w:.2f}$, $r_0 = {r0_w:.1f}$ pc, $\sigma^2 = {sig2_w:.0f}$ (km/s)$^2$"
ax.fill_between(
rgrid,
bfunc(rgrid, r0_w, sig2_w, m_w - 0.1),
bfunc(rgrid, r0_w, sig2_w, m_w + 0.1),
color="k", alpha=0.05,
)
ax.plot(rgrid, bfunc(rgrid, r0_w, sig2_w, m_w), lw=0.5, color="k", alpha=0.5, label=flabel_w)
# Plot points from weighted strucfunc
ax.plot(r, B_w, 'o', ms=3, alpha=0.5, label="Flux-weighted")
melnick_r = np.array([2.5, 7.5, 12.5, 17.5, 22.5, 27.5])
melnick_B = np.array([2.0, 2.1, 2.2, 2.2, 2.25, 2.25]) * 18.2**2
ax.plot(melnick_r, melnick_B, 's', label="Melnick+ (2020)", color="y", zorder=-10)
ax.axhline(sig2, color="k", ls="--")
ax.axhline(sig2_w, color="r", ls=":")
ax.legend(title="30 Doradus [N II]")
ax.set(
xscale = "log",
yscale = "log",
ylim = [0.5, 1500],
xlabel = "Separation, pc",
ylabel = r"$B(r)$, (km/s)$^2$",
)
fig.tight_layout()
sns.despine()
fig.savefig(figfile)
fig.savefig(figfile.replace(".pdf", ".jpg"))Same but in parsec.
Distance = 50 kpc => 1 arcsec = 0.242 pc
import json
import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns
data = json.load(open("tarantula-strucfunc.json"))
fig, ax = plt.subplots(figsize=(5, 5))
figfile = "tarantula-strucfunc-plot-parsec.jpg"
pixscale = 0.2 # arcsec
pixscale *= 0.242 # parsec
r = pixscale * 10**np.array(data["log10 r"])
B = np.array(data["Unweighted B(r)"])
sig2 = data["Unweighted sigma^2"]
B_w = np.array(data["Weighted B(r)"])
sig2_w = data["Weighted sigma^2"]
ax.plot(r, B, 'o')
ax.plot(r, B_w, 'o')
ax.axhline(sig2, color="k", ls="--")
ax.axhline(sig2_w, color="r", ls=":")
ax.set(
xscale = "log",
yscale = "log",
xlabel = "Separation, parsec",
ylabel = r"$B(r)$, (km/s)$^2$",
)
fig.tight_layout()
sns.despine()
fig.savefig(figfile)
fig.savefig(figfile.replace(".jpg", ".pdf"))What we got from test with pure python version with n = 20
{'log10 r': array([0. , 0.15, 0.3 , 0.45, 0.6 , 0.75, 0.9 , 1.05, 1.2 ]), 'Sum dv^2': array([ 286.176572 , 1816.08667744, 4456.64384807, 10247.69199721,
15614.41156876, 17353.22440817, 39015.36434385, 40118.31285852,
2608.48845663]), 'Sum weights': array([ 27.48855598, 50.72646818, 100.8550699 , 152.64575999,
181.63544682, 149.66641101, 179.43422474, 99.73101336,
13.14206985]), 'Sum w * dv^2': array([ 188.07159422, 1119.13494059, 2684.97040556, 6123.48035994,
8857.18554824, 9056.58693937, 18557.14879415, 17924.77263585,
1209.60779833]), 'N pairs': array([ 52, 98, 191, 283, 345, 301, 382, 223, 28]), 'Unweighted B(r)': array([ 5.50339562, 18.53149671, 23.33321386, 36.21092579,
45.25916397, 57.65190833, 102.13446163, 179.90274824,
93.16030202]), 'Weighted B(r)': array([ 6.84181426, 22.06214982, 26.62206678, 40.11562693,
48.76352993, 60.51182011, 103.42034147, 179.73117922,
92.04088944])}
- So far I have had medium-size data sets, which have allowed the use of inefficient algorithms using pandas
- Sizes of datsets
- Damiani:2016a Carina
- 866 spatial points
- Damiani:2017b Lagoon
- 1177 spatial points
- Castro:2018a Tarantula (30 Dor)
- 649 x 649 pixels => 421201 points
- This is 400 x what I was doing before
- Estimate of memory requirement to store all the pairs (assume 8 bytes per number)
N points pairs = N^2 / 2 Memory (GB) 1000 500000 0.0037 421000 88620500000 660.2742 - So that is not feasible - need to use the more efiicient algorithm.
- Based on my tests so far, I estimate that the numba parallel algorithm should be able to do the Tarantula structure function in 27.2 (649/400)**4 = 188.5 s, or 3 minutes
- So that is fine
- Damiani:2016a Carina
import sys
from pathlib import Path
from matplotlib import pyplot as plt
import seaborn as sns
import json
import numpy as np
from astropy.io import fits
sys.path.append("../muse-strucfunc")
import strucfunc
try:
LINEID = sys.argv[1]
except:
LINEID = "ha"
try:
METHOD = sys.argv[2]
except:
METHOD = "standard"
USE_COLDEN = "colden" in METHOD
USE_DEPTH = "depth" in METHOD
fitsfilename = {
"ha": "GAUS_Ha6562.8_060_Will.fits",
"nii": "GAUS_NII6583.45_060_Will.fits",
}
wav0 = {"ha": 6562.8, "nii": 6583.45}
atm_wt = {"ha": 1.0, "nii": 14.0}
fs_var = {"ha": 10.233, "nii": 0.0}
# Assume 1e4 K for thermal broadening
thermal_var = 82.5 / atm_wt[LINEID]
datadir = Path("data/Tarantula/MUSE_R136toWill")
hdulist = fits.open(datadir / fitsfilename[LINEID])
n = None
sb = hdulist[1].data[:n, :n].astype(np.float64)
vv = hdulist[2].data[:n, :n].astype(np.float64)
ss = hdulist[3].data[:n, :n].astype(np.float64)
# optionally use column density, instead of surface brightness
if USE_COLDEN:
dd = fits.open(datadir / "Density.fits")["DATA"].data[:n, :n].astype(np.float64)
sb /= dd
if USE_DEPTH:
dd = fits.open(datadir / "Density.fits")["DATA"].data[:n, :n].astype(np.float64)
sb /= dd**2
# Convert sigma to km/s
ss *= 3e5 / wav0[LINEID]
# Subtract instrumental width and thermal width
ss = np.sqrt(ss**2 - 48.0**2 - fs_var[LINEID] - thermal_var)
# Replace spurious values in the arrays
m = ~np.isfinite(sb*vv*ss) | (sb < 0.0)
if LINEID == "nii":
# Remove bad patch from the [N II] map
m = m | (sb > 6e4)
m = ~m # invert mask
# additional mask for bright pixels
# BRIGHT_THRESHOLD = 0.1*np.max(sb[m])
BRIGHT_THRESHOLD = np.median(sb[m])
mb = sb > BRIGHT_THRESHOLD
# Brightness-weighted average sigma
AV_SIG = np.average(ss[m], weights=sb[m])
NBIN = 100
BMAX = np.max(1.2*sb[m])
BMIN = BMAX / 1000.0
if USE_COLDEN:
BMAX = 5*BRIGHT_THRESHOLD
BMIN = BMAX / 100.0
if USE_DEPTH:
BMAX = 10*BRIGHT_THRESHOLD
BMIN = BMAX / 500.0
SMIN, SMAX = 0.0, 90.0
VMIN, VMAX = 220.0, 330.0
GAMMA = 1.5
vlabel = "Centroid velocity, km/s"
slabel = "RMS line width, km/s"
blabel = "log10(Surface brightness)"
if USE_COLDEN:
blabel = "log10(Column density)"
if USE_DEPTH:
blabel = "log10(LOS depth)"
fig, axes = plt.subplots(2, 2)
linestyle = dict(lw=0.7, ls="--", color="r", alpha=0.5)
# I - sigma
xmin, xmax = np.log10(BMIN), np.log10(BMAX)
ymin, ymax = SMIN, SMAX
H, xedges, yedges = np.histogram2d(
np.log10(sb[m]), ss[m],
bins=[NBIN, NBIN],
range=[[xmin, xmax], [ymin, ymax]],
)
axes[0, 0].imshow(
(H.T)**(1.0/GAMMA),
extent=[xmin, xmax, ymin, ymax],
interpolation='none', aspect='auto',
origin='lower', cmap=plt.cm.gray_r,
)
# Show brightness thereshold
axes[0, 0].axvline(np.log10(BRIGHT_THRESHOLD), **linestyle)
# Show average sigma
axes[0, 0].axhline(AV_SIG, **linestyle)
axes[0, 0].set(
xlabel=blabel,
ylabel=slabel,
xlim=[xmin, xmax],
ylim=[ymin, ymax],
)
# I - V
xmin, xmax = np.log10(BMIN), np.log10(BMAX)
ymin, ymax = VMIN, VMAX
H, xedges, yedges = np.histogram2d(
np.log10(sb[m]), vv[m],
bins=[NBIN, NBIN],
range=[[xmin, xmax], [ymin, ymax]],
)
axes[1, 0].imshow(
(H.T)**(1.0/GAMMA),
extent=[xmin, xmax, ymin, ymax],
interpolation='none', aspect='auto',
origin='lower', cmap=plt.cm.gray_r,
)
# Show brightness thereshold
axes[1, 0].axvline(np.log10(BRIGHT_THRESHOLD), **linestyle)
axes[1, 0].set(
xlabel=blabel,
ylabel=vlabel,
xlim=[xmin, xmax],
ylim=[ymin, ymax],
)
# V - sigma
xmin, xmax = VMIN, VMAX
ymin, ymax = SMIN, SMAX
H, xedges, yedges = np.histogram2d(
vv[m & (~mb)], ss[m & (~mb)],
bins=[NBIN, NBIN],
range=[[xmin, xmax], [ymin, ymax]],
)
# Show average sigma
axes[0, 1].axhline(AV_SIG, **linestyle)
axes[0, 1].imshow(
(H.T)**(1.0/GAMMA),
extent=[xmin, xmax, ymin, ymax],
interpolation='none', aspect='auto',
origin='lower', cmap=plt.cm.gray_r,
)
axes[0, 1].set(
xlabel=vlabel,
ylabel=slabel,
xlim=[xmin, xmax],
ylim=[ymin, ymax],
)
# V - sigma but bright only
xmin, xmax = VMIN, VMAX
ymin, ymax = SMIN, SMAX
H, xedges, yedges = np.histogram2d(
vv[m & mb], ss[m & mb],
bins=[NBIN, NBIN],
range=[[xmin, xmax], [ymin, ymax]],
)
axes[1, 1].imshow(
(H.T)**(1.0/GAMMA),
extent=[xmin, xmax, ymin, ymax],
interpolation='none', aspect='auto',
origin='lower', cmap=plt.cm.gray_r,
)
# Show average sigma
axes[1, 1].axhline(AV_SIG, **linestyle)
axes[1, 1].set(
xlabel=vlabel,
ylabel=slabel,
xlim=[xmin, xmax],
ylim=[ymin, ymax],
)
fig.tight_layout()
plotfile = f"tarantula-I-sigma-hist-{LINEID}.png"
if USE_COLDEN:
plotfile = plotfile.replace(".", "-colden.")
if USE_DEPTH:
plotfile = plotfile.replace(".", "-depth.")
fig.savefig(plotfile, dpi=200)
print(plotfile, end="")
python tarantula-I-sigma-hist.py hapython tarantula-I-sigma-hist.py niipython tarantula-I-sigma-hist.py ha coldenpython tarantula-I-sigma-hist.py ha depth- General definition:
\[
\mathrm{Re} = \frac{u L}{ν}
\]
- u is velocity
- L is size
- ν is kinematic viscosity in cm^2/s
- ν = λ_mfp u_therm
- viscosity is dominated by electrons with u_therm ≈ 550 T_41/2 km/s
- u_therm ≈ (m_p/m_e)1/2 c_s
- from kappa paper: λ_mfp ≈ 1.3e12 T_4^2 n_e^-1 cm
- So Re = (u / u_therm) (L / λ_mfp)
- Relation to Knudsen number and Mach number:
- Kn = λ_mfp / L
- Ma = u / c_s
- => Re ≈ (m_e/m_p)1/2 Ma / Kn ≈ 0.023 Ma / Kn
- Relation to ionization parameter
- U = Q / 4 π R_s^2 n c
- Q is ionizing luminosity; R_s is Strömgren radius
- From the kappa paper, I have:
- Kn = 3.76e-12 U^-1
- Also U = 0.0006 (Q_49 n)1/3
- Hence
- Re ≈ 6.1e9 Ma U
- Re ≈ 3.7e6 Ma (Q_49 n)1/3
- Typically, Q_49 n ≈ 1e4 => U ≈ 0.013
- Re ≈ 8e7 Ma
- U = Q / 4 π R_s^2 n c
- Taylor scale Reynolds number
- This is Re* in Elsinga:2020a
- They have an equation (2.4) that relates Re* to the global Re
- we use it with their default parameters
- D = 0.5
- α = 0.010
- b = 0.67
- To get
- Re* = 9.65 Re1/3 [0.67 + 0.33 0.01 Re]1/6
- For Re >> 1000, this becomes 3.7 Re1/2
- So Re = 1e8 -> Re* = 3.7e4
- If we substitute this again (full equation this time), we get
- Re** = 717 for sub-layers
- Re*** = 104 for sub-sub layers
- So fact that sub-sub layers have Re*** < 150 means that they are not themselves turbulent.
- we use it with their default parameters
- “Significant shear layers” from Elsinga:2020a
- They have thickness pf 4 λ_T
- And velocity difference of order U, which is the large-scale velocity dispersion
- So this velocity difference is much larger than the RMS b = |v-v’| at scale λ_T
- because the volume filling factor of the shear layers is small, of order λ_T/L
- This is repeated for the sub-layers and sub-sub layers
- they all have velocity differences of order U
- but smaller and smaller volume fractions
- Wikipedia equations for Taylor microscale
- This gives λ_T / L = sqrt(10 / Re)
- For Re = 1e8, this gives λ_T = 3e-4 L
- This is consistent with the result Re* = 3.7e4 Re, assuming same U at all scales
- Also, η / L = Re-3/4 = 1e-6 L
- But this is the Kolmogorov length at the mean dissipation rate
- The Elsinga:2020a theory is that the shear layers have enhanced dissipation and smaller and smaller Kolmogorov scales with each nesting
- So if the driving scale is 1 pc, then we have a Taylor microscale of 0.3 mpc, or 0.15 arcsec at Orion
- And Orion has a smaller driving scale, so Taylor microscale would be tough to observe
- Hierarchy
- L -> L* -> L** -> L***
- L = 1 pc -> L* = 0.3 mpc -> L** = 0.005 mpc (1 AU) -> L*** = 0.0006 mpc (0.12 AU)
- η = 0.001 mpc (0.2 AU) -> η* = 0.0001 mpc (0.02 AU) -> η** = 0.00003 mpc (0.006 AU) -> stop there!
- This gives λ_T / L = sqrt(10 / Re)
- Conclusion
- Typical Reynolds number should be of order 1e8 on scale of Strömgren radius
- If injection scale is smaller than that, then will be reduced accordingly
- From the analysis of Elsinga:2020z this means that the flow should support multiple nested shear layers, up to sub-sublayers
- The Reynolds number at the Taylor scale should be smaller, but only by a factor of l_T/L
- This is because the velocity difference is the same at all nested layers
- Typical Reynolds number should be of order 1e8 on scale of Strömgren radius
- https://en.wikipedia.org/wiki/Integral_length_scale
- This has a definition in terms of the autocorrelation C(l):
- L = ∫_0^∞ C(l) dl
- We had been using C(l) = [1 + (l/l_0)^m]^-1
- This has the problem that L does not converge for m <= 1
- If we fudge this by using an outer scale of A l_0, then we get
- L/l_0 = A 2F1(1, 1/m, (1 + m)/m, -A^m)
- With A = 10, we get L/l_0 = [3.0, 2.3, 1.6] for m = [0.7, 1.0, 1.8]
- With A = 100 we get L/l_0 = [8.5, 5, 2] for m = [0.7, 1.0, 1.8]
- What if we used C(x) = [1 + (1/a) x^m]^-a with x ≡ l/l_0 ?
- Small x expansion is still
- C(x) = 1 - (a/a) x^m = 1 - x^m
- Take case a = 2
- Indefinite integral is now: x 2F1(2, 1/m, (m + 1)/m, -0.5 x^m) which should converge to finite value as x → ∞
- Yes, limit is [5.31, 2.00, 1.16] for m = [0.7, 1.0, 1.8]
- Small x expansion is still
- Can we not modify the equation so that we get better behavior for l > l_0 ?
- S(x) = 2 [1 - C(x)] with x ≡ l/l_0
- Small x: 2[ 1 - (1 + x^m)^-1] = 2 [1 - (1 - x^m)] = 2 x^m
- So multiply that by something that is ≈ 1 for small x, but falls as x^-m at large x
- 1 / 1 + x^m would get us back to original
- e^-x + (1 - e^-x) x^-m might work
- Small x: 1 - x + (1 - 1 + x) x^-m = 1 - x + x1-m
- So that is OK if m < 1
- S(x) = 2 [1 - C(x)] with x ≡ l/l_0
- Alternative definition in terms of power spectrum P(k)
- L = ∫_0^∞ k^-1 P(k) dk / ∫_0^∞ P(k) dk
- [ ] or should I be integrating over d^3k ?
- Yes, I think so. Pope has a similar equation, but with a constant \[ L11 = \frac{π}{2 \langle u_1^2 \rangle} ∫_0^∞ \frac{E(κ)}{κ} \, dκ \]
- And the E(κ) that he uses is one where Kolmogorov spectrum is κ-5/3 so this is different from the k^-β version
- In the Lazarian terminology, the power spectrum is k^n with n = -11/3 for Kolmogorov
- In Pope terminology the energy spectrum is E(κ) ~ κ^-p with p = 5/3 for Kolmogorov
- So n = -(p + 2)
- In Lazarian terminology, “steep” is n < -3, which corresponds to p > 1 in Pope
- In Pope sec 6.1.3, there is a discussion of the possible limits on p if it is extended over the entire wavenumber range
- This shows that p should be in range 1 < p < 3 to avoid high-κ divergence of energy (on shallow side) or low-κ divergence of dissipation (on steep side)
- This implies that -n is in the range 3 to 5
- Of course, it is possible to have a shallower spectrum over a certain range of κ, so long as there is a high-κ cut-off
- In Pope sec 6.1.3, there is a discussion of the possible limits on p if it is extended over the entire wavenumber range
- Relation with turbustat terminology
- For our simulated fields we use P(k) = k^-β
- but only between k = 1 and k = N
- with k being spatial frequency in units (box size)^-1
- Therefore L = ∫_1^N k-(β+1) dk / ∫_1^N k^-β dk
- ∫_1^N k-(β+1) dk = [-1 / β] (N^-β - 1) = (1 - N^-β) / β
- ∫_1^N k^-β dk = [-1 / (β - 1)] (N-(β-1) - 1) = (1 - N-(β-1)) / (β - 1)
- L = [(β - 1) / β] [(1 - N^-β) / (1 - N-(β-1))]
- For large N, the second term ≈ 1
- So L ≈ (β - 1) / β ≈ (m + 1) / (m + 2)
- So for m = [0.7, 1.0, 1.8], we have L = [0.63, 0.67, 0.74]
- I always get confused by this
- Claims from Ossenkopf:2006a
- Intro
- Δ-variance slopes
- Δ-variance amplitudes
- Figure (note, I only show the lower plot - the upper plot is for a shallow spectrum, but it is just the same really)
- Text to go with Figure
Note especially the comment that they only show scales up to about 1/3 of the box size. - Further commentary by Will
- The figure looks pretty convincing and shows the effect of the projection smoothing, which steepens the power law by 1
- This is because they are projecting the entire cube
- I suspect that using a thin layer of a cube would produce different results
- Intro
- Structure functions of molecular gas
- Henshaw:2020t
- Ubiquitous velocity fluctuations throughout the molecular interstellar medium
- Carina
- Line widths typically have σ = 15 km/s, but this includes the thermal contribution.
- thermal σ is 9 km/s => non-thermal is 12 km/s
- For individual components, we have 5.6 km/s, so we have σ(LOS) = 2 σ(POS), as in Orion
- 30 Doradus
- We have σ = 55 km/s for H alpha at highest brightness
- For Orion with the same instrument we have 48.5 km/s
- Subtracting in quadrature gives 26 km/s for LOS σ
- Note however, that 15 km/s is what Melnick et al found for the same region
- Is this because they have decomposed the lines into components first
- I need to check this with the original data
- Note however, that 15 km/s is what Melnick et al found for the same region
- For POS σ we have sqrt(252) = 16 km/s
- So, we have LOS and POS σ being the same if we believe Melnick, although the MUSE widths suggest LOS σ is 2 x higher
- NGC 346
- Ha widths are 49 +/- 1 in the regions of moderate brightness
- Subtracting off the Orion value (48.5) in quadrature gives 7 (+5) (-7)
- Although this assumes the Orion width is zero, which is not entirely true
- If the Orion width is really 8, then this means that the instrumental width is actually sqrt(48.5**2 - 8**2) = 47.8
- So we get 11 +4 -7 km/s
- As compared with POS sigma = 5.65 km/s
- So this is similar to the Orion and Carina case, with σ(LOS) = 2 σ(POS)
- Note that there are higher σ at the positions of the globule heads:
- Highest is 64 km/s raw value => sqrt(64**2 - 47.8**2) = 42 km/s
- But these are a negligible fraction of the area
- Higher resolution spectroscopy from Danforth:2003m
- For Ha they find a corrected FWHM of 23.8 km/s => σ(LOS) = 10.1 km/s
- This consistent with the MUSE value, but is much more precise
- Orion
- M8
| Q_H | L(Ha) | SFR | n | R_S, pc | L, pc | Σ_SFR | \ell_0, pc | σ | m | D, kpc | |
|---|---|---|---|---|---|---|---|---|---|---|---|
| Orion | 1e49 | 1.2e37 | 5.3e-5 | 4000 | 0.187 | 0.6 | 46.862 | 0.05 | 3.1 | 1.1 | 0.4 |
| Orion large | 1e49 | 1.2e37 | 5.3e-5 | 100 | 2.184 | 4 | 1.054 | 0.4 | |||
| M8 small | 2e49 | 2.3e37 | 1.0e-4 | 600 | 0.833 | 20 | 0.080 | 1.3 | 3 | 0.8 | 1.3 |
| M8 large | 2e49 | 2.3e37 | 1.0e-4 | 60 | 3.868 | 20 | 0.080 | 6 | 4 | 1 | 1.3 |
| Carina | 2e51 | 2.3e39 | 1.0e-2 | 500 | 4.368 | 15 | 14.147 | 0.5 | 4 | 0.7 | 2.0 |
| NGC 346 | 4e50 | 4.7e38 | 2.1e-3 | 100 | 7.469 | 64 | 0.163 | 1.5 | 5.6 | 1.05 | 61.7 |
| 30 Dor | 5e51 | 5.8e39 | 2.6e-2 | 250 | 9.411 | 98.9 | 0.846 | 2.7 | 16 | 1.22 | 50 |
| Hubble V | 3e49 | 3.5e37 | 1.5e-4 | 90 | 3.379 | 100 | 0.005 | 3.6 | 2.8 | 1.8 | 500 |
| Hubble X | 6e49 | 7.0e37 | 3.1e-4 | 30 | 8.856 | 150 | 0.004 | 4.7 | 3.6 | 1.7 | 500 |
| NGC 604 T | 1e51 | 1.2e39 | 5.3e-3 | 50 | 16.092 | 400 | 0.011 | 11 | 7.3 | 1.7 | 840 |
| NGC 595 | 5e50 | 5.8e38 | 2.6e-3 | 50 | 12.772 | 300 | 0.009 | 11 | 6.6 | 1.7 | 840 |
More recent census of stars in 30 Doradus gives slightly larger QH - Bestenlehner:2020a give 2.75e+51 /s
L(Ha) = (hν) α_Ha VEM Q_H = α_B VEM = L(Ha) α_B/α_eff 1/hν
hν = 13.6 (1/4 - 1/9) eV = 3.026e-12 erg α_B / α_eff = 2.6
Q_H = 8.56e+11 L(Ha)
NGC 604: L(Ha) = 1e39 => Q_H = 1e51 erg NGC 595: L(Ha) =
Hubble V: L(Ha) = 1e49 en fotones => Q(H) = 2.6e49 Hubble X: L(Ha) = 2.4e49 en fotones => Q(H) = 6.24e49
- Our Galaxy:
- Orion
- M8
- Carina
- Local Group
- LMC
- 30 Dor
- NGC 6822
- Hubble V
- Hubble X
- M 33 (Triangulum Galaxy)
- NGC 595
- NGC 604
- LMC
- These are all regions with MUSE data
| Source | Distance, kpc | Lum | Size |
|---|---|---|---|
| M17 | |||
| RCW 49 (Wd 2) | 4.0 | >5e50 | |
| LMC-N44 | 50 | 2.5e50 | 45 |
| LMC-N180 | 50 | 1.3e50 | 40 |
- RCW 49 parameters
- Zeidler:2018a
- Zeidler:2021a has the velocity maps
- Shows a coherent pattern from blueshifted on the outside to redshifted on the inside
- Gradient of 15 to 20 km/s
- Seems correlated with the extinction
- Global expansion of H II region
- Tiwari:2021a has C I shell, like in Orion
- Also expanding at 15 km/s
- They say it is stellar wind driven but I am not convinced
- N44 and N80 parameters
- Table 2 of McLeod:2019a gives Ha fluxes:
- 7.2e-10 erg/s/cm2 for N44
- 4.5e-10 erg/s/cm2 for N180
- Not corrected for extinction
- D = 50 kpc => 4 π D^2 = 2.991e+47
- => L(Ha) = 2.0e38 and 1.3e+38
- So these are 15 and 20x Orion
- This plugs a useful gap in our samples
- This is confirmed by stellar census
- McLeod:2019a sections 3.1 and 3.2
- Q(H) = 2.5e50 s^-1 for N44
- Q(H) = 1.3e50 s^-1 for N180
- Note that there are clear sub-regions to each of these regions
- N44 A has 3.2e49 (O7III, O5V, O8V)
- N44 B has 2.6e48 (O9.5III)
- N44 B1 has 7.6e48 (O8III)
- N44 C has 3.3e49 (O5III, O8V, O9.5V)
- N44 D has 1.7e49 (O5V)
- N44 Main has the rest: 1.6e50 (WN4b + O8.5I + …)
- N180 similarly has 4 sub-regions
- Some may be too small to measure the structure function
- And some seem to be badly affected by extinction
- There is also a per-region L(Ha) analysis
- McLeod:2019a Table 7 that compares the Q(H) from stars and from recombinations, finding escape fraction of 0.3 to 0.8
- Table 2 of McLeod:2019a gives Ha fluxes:
- [2021-05-04 Tue] I am downloading 3 sample datasets:
- LMC-N44 strip: ADP.2016-08-11T18:26:36.424
- LMC-N180 strip: ADP.2016-08-12T09:14:55.932
- SMC-NGC346 square: ADP.2017-10-16T11:04:19.247
- The NGC 346 data is a restricted FOV of 1.5 arcmin in center of cluster, but it seems much higher data quality than the others
- Munoz-Tunon:1996a has TAURUS spectroscopy of NGC 595 and 588
- They present I-sigma maps
- Gives a high-intensity minimum corrected sigma(LOS) for NGC 604 of 18 km/s
- Compare our sigma(POS) of 7.3 km/s
- Compare mainly with cometary stirring model of Tenorio-Tagle:1993e
- Hodge:1989c calculate luminosity function of H II regions in NGC 6822
- For Hubble V and X they find
- log F(Ha) = 6.95, 6.93
- These are the brightest two regions in the galaxy
- In units of 1e-18 erg/s/cm2
- So about 8.9e-12 erg/s/cm2
- Which becomes 4e38 erg/s
- [ ] This is 10x larger than what we had in the table - need to fix!
- Isophotal diameter = 55, 81 pc
- Although up to twice as large if they take a lower brightness isophote
- log F(Ha) = 6.95, 6.93
- For Hubble V and X they find
- Bresolin:2020a do an I-sigma analysis of H II regions in M 101
- In principle they could have done structure functions but didn’t
- They also discuss the difference in width between [O III] and Hβ
- But they seem to be largely wrong about it
- Gaia-ESO spectroscopy
- Damiani:2016a Carina
- 10 arcmin at 2.3 kpc = 6 pc
- 10 arcsec = 0.1 pc
- Can analyze blue and red components separately and together
- Damiani:2017b Lagoon
- D = 1250 pc
- Total extent: 0.8 deg = 17 pc
- Damiani:2016a Carina
- MUSE spectroscopy
- McLeod:2016a Carina
- 2 arcmin at 2.3 kpc = 1.3 pc
- 1 arcsec = 0.01 pc
- Mc-Leod:2016a Orion
- 5 arcmin at 410 pc = 0.6 pc
- 1 arcsec = 0.002 pc
- Amplitude is about 10 km/s
- Castro:2018a Tarantula (30 Dor)
- These are excellent
- 2 arcmin at 50 kpc = 30 pc
- 1 arcsec = 0.242 pc
- Amplitude is about 40 km/s
- Looks like 5 arcsec is typical fluctuation scale => 1 pc
- Anti correlation of intensity with σ(LOS) - see Fig 7
- Could use density and de-extincted Hα to get LOS thickness
- McLeod:2015a Eagle pillars
- A bit too noisy to do much with
- McLeod:2016a Carina
- Longslit spectroscopy
- Arthur:2016a Orion
- Find scale of 0.05 pc decorrelation scale from struc func
- Similar scale from power spectrum of brightness
- Inner scale of 0.02 pc, but not clear what this is
- Arthur:2016a Orion
- Also see the images in introduction to my Greece talk
- Tiwari:2018a have CO, etc for Lagoon region
- Ionized gas is at negative velocities with respect to CO
- Lada:1976a have +11 km/s LSR for CO at biggest clump
- This is reminiscent of champagne flow, as in Orion
- Ionized gas is at negative velocities with respect to CO
- Esteban:1999a have multiple optical lines for M42 (Orion), M17 (Eagle), and M8 (Lagoon)
- For Lagoon
- Slit is 25 arcsec S of Hourglass
- V(HEL) approx -2 km/s for low-ionization lines
- Higher ionization lines are more like -10
- Discrepant results for [O I] and [N I] but may suffer from sky line contamination
- They find a velocity-IP correlation in all cases, indicating blue-shifted champagne flows
- For Lagoon
- Looks like main ionizing star (9 Sgr) is 1pc in front of cloud in Lagoon
- From Fig 20 of Damiani:2017a
- Chakraborty:1999a calculated the [O III] structure function of Hourglass region
- Separations of 2-30 arcsec (so very little dynamic range)
- As compared with up to 1800 arcsec in Damiani
- They did something very strange to eliminate large scale gradients
- And their absolute values are very large: saturates at 280 km^2/s^2
- Chakraborty:1997a give [N II] velocity image but do not calculate statistics
- And they don’t even have absolute velocities
- Haenel:1987a did FP spectroscopy of the entire nebula and has a grid at 50 arcsec resolution (0.014 deg)
- 34 x 26 pixels = 0.5 x 0.4 degrees
- Velocities are LSR
- Seems to agree more or less with Damiani, given the 10 km/s difference between heliocentric and LSR
- So this is very similar number of points to the Damiani data, but coverage is more uniform
- Damiani Carina
- Damiani Lagoon
- Jupyter notebook: mariano-test.ipynb
- Pure python version: mariano-test.py
- [ ] Install required packages
- [ ] Load data with astropy
- [ ] Convert to pandas
- [ ] Clean up data as necessary
- [ ] Look at correlations
- [ ] Make maps
- [ ] Calculate structure functions
- McLeod:2016a best region is around “defiant finger”, just to W of Keyhole.
- That is the brightest region, and only one that overlaps the Gaia-ESO observations
- Other regions are fainter and are in the periphery of the nebula
- As well as the published McLeod:2016a stuff, there are new observations of the Tr14 region, which are available from the data archive
- There is a python package for working with MUSE data:
mpdaf, which might or might not be useful- There is the option of working with pixel tables, which have not been resampled
- This might help avoid some of the artefacts seen in the velocity maps
- This is in data/Tarantula/MUSE_R136toWill
- Example of extracting coordinates for each pixel in the velocity maps is given in
- Haenel:1987a have maps at arcmin scale for whole nebula
- [ ] Could extend velocity statistics for Orion by combining this with the Arthur:2016a Garcia-Diaz:2008a data
- Look at data like in Moiseev:2015a
- Damiani 2017 Lagoon
- Damiani 2016 Carina
- Arthur 2016 Orion
- Medina 2014 Simulaciones
- García-Díaz 2008 Orion