Adam Ginsburg's blog

PPVI Talk 3: Paul Clark

2013/07/15
THE ORIGIN AND UNIVERSALITY OF THE INITIAL MASS FUNCTION

Authors:

S. Offner (Yale University, United States), P. Clark (Heidelberg University, ITA, Germany), P. Hennebelle (Observatoire de Paris, LERMA, Paris, France), N. Bastian (Liverpool John Moores University, United Kingdom), M. Bate (University of Exeter, United Kingdom), P. Hopkins (California Institute of Technology, Department of Astronomy, United States), E. Moraux (Institut de Planetologie et d'Astrophysique de Grenoble (IPAG), France), A. Whitworth (Cardiff University, Physics & Astronomy, United Kingdom)

Abstract

We review current theories for the origin of the Stellar Initial Mass Function (IMF) with particular focus on the extent to which the IMF can be considered universal across various environments. To place the issue in an observational context, we summarize the techniques used to determine the IMF for different stellar populations, the uncertainties affecting the results, and the evidence for systematic departures from universality under extreme circumstances. We next consider theories for the formation of prestellar cores by turbulent fragmentation and the possible impact of various thermal, hydrodynamic and magneto-hydrodynamic instabilities. We address the conversion of prestellar cores into stars, and evaluate the roles played by different processes: competitive accretion, dynamical fragmentation, ejection and starvation, filament fragmentation and filamentary accretion flows, disc formation and fragmentation, critical scales imposed by thermodynamics, and magnetic braking. We present explanations for the characteristic shapes of the Present-Day Prestellar Core Mass Function and the IMF and consider what significance can be attached to their apparent similarity. Substantial computational advances have occurred in recent years, and we review the numerical simulations that have been performed to predict the IMF directly and discuss the influence of dynamics, time-dependent phenomena, and initial conditions.

Last time, IMF discussion was very heated. (technical difficulties)
Intro to IMF

Salpeter, Kroupa, Chabrier, de Marchi ("Tapered" power law) (aside) Thies & Kroupa: BDs not part of IMF

SFR indicators only see >10 Msun sources. Must extrapolate downward. This is why extragalactic people care.
Are there variations in the IMF?

Locally, no.

More evolved clusters, yes on the low-mass end:

Pleiades, NGC 2516, etc.

Lose low-mass stars with age dynamically

fraction of BDs remaining drops to zero over time

Difficulties in measurement:

Patchy extinction

young clusters not done sampling IMF

don't understand hayashi tracks that well

Modal plot from van Dokkum

Hints at change in IMF

Extragalactic - big camps argue top, bottom heavy and universal

rely on mass to light ratios

M/L models depend on population synthesis models
What determines the IMF?

Hoyle, Larson: the cloud fragmentation process. Clouds fragment into cores, cores make stars.

Does supersonic turbulence regulate? (movie from Federrath)

Padoan:

CMF driven only by turbulence.

# of BDs controlled by Mach number.

Cloud-in-cloud problem: small cores may be in big cores

Press-Schechter & Excursion Set formalism

Gravitational decoupling sets "critical" scale

Hopkins: Density as a function of scale

large scale: shear support

intermediate: turbulence

small scale: hrm...

last-crossing => IMF

Hennebelle & Chabrier 2008 - account for non-thermal support

GMCs are a separate population

Useful: turbulence -> more BDs, so extragalactic maybe more turbulent?
Does the CMF map onto the IMF?

Cores must:

neither accrete nor merge (or you break the CMF my moving cores)

fragment self-similarly: same number of fragments per core

Same efficiency regardless of mass

formation must happen on similar timescale (can't have high-mass first)

Cores must be able to form stars, i.e. must be grav. bound

Can cores form stars?

Dust cores depend on column & temperature & line of sight effects

number of Jeans masses goes as T^-3

Velocity structure? B-fields?

Fragmentation within cores

Used to think that everything fragmented...

strong B-fields reduce fragmentation

Still haven't studied parameter space well

SFE within cores

Binaries or triples form, requires 30%-50% efficiency

can work with outflows

but not clear... may depend on well-aligned B-field

Problems with CMF:

Pipe cores follow CMF shape, but are pressure bound

CMF in Aquila is too low

need CMF of Class 0 sources
Clusters

What effects have cluster simulations had on IMF studies...

Initial Conditions imprints thermal jeans mass in box, which imprints onto IMF

thermodynamics matter

Fraction of cloud doesn't matter

IMF is same in small periodic box and isolated cloud

colliding flows can build clouds, be more consistent

How we drive turbulence doesn't seem to affect IMF (solenoidal, compressive same)

power spectrum of turbulence doesn't affect IMF

Metallicity?

Myers 2011: 20x metallicity change -> same IMF

Multiplicity may change (Bate)

even 10^-5 Z_sun -> same IMF

Things that DO affect IMF:

radiative transfer: reduces fragmentation

need wind feedback, otherwise radiation -> flat IMF by preventing fragmentation but keeping gas

B-fields reduce characteristic mass by creating more filaments with different jeans masses

State of Simulations:

Top of IMF about right

bottom is wrong
Secondary Metrics

Multiplicity

Krumholz & Bate both reproduce observations fairly well

Require radiative transfer
Summary

IMF studies focused on low-mass end instead of high

Cores better understood

Future: Put simulations into observational plane
Questions
Q: Properties of turbulence in simulations give same IMF...

A: Outflows drive turbulence, but turbulence doesn't change IMF. Outflows remove gas.

Q Zinnecker: Binaries.

A: Cathy will talk about this. Fragmentation knobs are the issue Unclear that binary properties will stay constant (does field match cluster?) N-body evolution matters

Q (Galactic center star studier...): not done with discussion of high-mass issues. Radiatitve transfer still matters.

A: Many effects are actually from mass segregation

Q: Alyssa Goodman: Are the CMF and IMF both lognormals?

A: Sometimes fitting distributions they expect. Maybe different processes in low mass clouds.

PPVI Talk 2: Philippe Andre

2013/07/15
FROM FILAMENTARY NETWORKS TO DENSE CORES IN MOLECULAR CLOUDS: TOWARD A NEW PARADIGM FOR STAR FORMATION

Authors:

P. Andre (CEA Saclay, Laboratoire AIM - Service d'Astrophysique, Gif-sur-Yvette, France), J. Di Francesco (National Research Council Canada, Herzberg Institute of Astrophysics, Canada), S.-I. Inutsuka (Nagoya University, Japan), R. Pudritz (McMaster University, Origins Institute, Canada), D. Ward-Thompson (University of Central Lancashire, Jeremiah Horrocks Institute, UK), J. Pineda (University of Manchester, Jodrell Bank Centre for Astrophysics, United Kingdom)

Abstract

We review recent progress in our understanding of the physics controlling the earliest evolutionary phases of star formation. Since PPV seven years ago, one area that has seen the most dramatic advances has been the characterization of the link between star formation and the structure of the cold interstellar medium (ISM). In particular, extensive studies of the nearest star-forming clouds of our Galaxy with the Herschel Space Observatory have provided us with unprecedented images of the initial and boundary conditions of the star formation process. The Herschel images reveal an intricate network of filamentary structures in every interstellar cloud. The observed filaments share common properties such as their central widths - but only the densest ones contain prestellar cores, the seeds of future stars. Overall, the Herschel submillimeter data, as well as other observations from, e.g., near- IR extinction studies, favor a scenario in which interstellar filaments and prestellar cores represent two key steps in the star formation process: first supersonic turbulence stirs up the gas, giving rise to a universal web-like structure in the ISM, then gravity takes over and controls the further fragmentation of filaments into prestellar cores and ultimately protostars. The new observational results connect remarkably well with nearly a decades worth of numerical simulations and theory which have consistently shown that the ISM should be highly filamentary on all scales and that star formation is intimately connected with self-gravitating filaments. We thus attempt to synthesize a comprehensive physical picture that arises from the confrontation of recent observations and simulations. We also emphasize how the apparent complexity of cloud structure and star formation may be governed by relatively simple universal processes - from filamentary clumps to galactic scales.
Content

Goal: convince that filamentary structures are essential

(Talk was relatively disorganized; this is reflected in the disorganization below)

Filaments well-known from both observations and theory.
IRAS image of Taurus

Filamentary IRDCs

Turbulence, collapse, etc. all show filamentary structures in simulations
Herschel is great!
G11.11-0.12 one of the densest, darkest clouds

Filamentary networks in Vela C

Sousbie's DisPerSE shows "networks of filaments"
Universals in filaments:
Striations perpendicular to large-scale filaments (e.g., Taurus)

Serpens S shows H-band polarization perpendicular to filaments

similar column density profiles: fixed "scales"

QUESTION: If these are genuinely common size scales, can they be used as a distance tracer? [implied answer: No, they are not universal, they are observationally selected]

Plummer-like density profile: r^-2, not the expected cylindrical Ostriker solution r^-4

colder inside (non-isothermal)

size of central inner plateau "roughly the same"

same over wide range of column densities
Possible explanation
0.1 pc is the scale of turbulence

in Larson's law, 0.1pc is where velocity dispersion meets sound speed

Sharp transition (Pineda 2011; Goodman 1998?)

Turbulent fragmentation may cause this?

width of compressed filaments is the sonic scale

Hennebelle: Compression + shear -> filaments...
Equilibrium with ISM?
unbound filaments below N ~ 5e21 (Inutsuka in prep)

Gravitating filaments above that sacel

widths centered at 0.1 +/- 0.05 pc?
More observations (Helen Kirk et al)
Radial infall motions in HNC along minor axis of filament

gradient along filament "suggesting infall" (no, that doesn't suggest infall, at least no more than it suggests rotation)

Overall picture is collapse in all dimensions

Accretion of "background material" required to maintain "fixed widths"
Arzoumanian 2013:
velocity dispersion increases with Column density

velocity dispersion gives mass per unit length...
Heitsch 2013:
Accretion onto dense filaments

Most likely location in parameter space is around 0.1 pc....
Hacar, Tafalla et al 2013 (1S036)
CO -> velocity-coherent "fibers" branching off of filament
Aquila:
Wavelet/Curvelet decomposition

wavelets extract cores well

curvelets highlight filaments
Konyves cores
core extraction comparison underway

Menschikov getsources, Kirk csar, Molinari cutex (csar: see this thesis? http://orca.cf.ac.uk/14483/ just used,didn't develop, csar)
Mass vs Size diagram for starless vs starry cores
Simpson: model evoluationary tracks

Prestellar CMF resembles IMF...

one-to-one mapping....

supports fragmentation models...

75% of cores form "on" filaments

preferentially form above A_V ~ 8

Instability proportional to mass per length for filaments...

16 msun/pc corresponds to 160 msun/pc^2 or 1600 msun/pc^3 (2e4 cm^3)

threshold explained by filaments...

Bonnor-Ebert mass corresponding to 16 msun/pc is ~0.5 pc, corresponding to CMF peak
Conclusions

First you form filaments Then filaments fragment Perhaps massive stars come from filament convergence into hubs "This scenario may possibly account for the global rate of star formation on Galactic scales"
Q: Chris: A_V = 8. Arises naturally?

Q: Is there any relationship between the critical column density for SF and the column density of the parent cloud

A: Not really. Average column density in clouds is very low. No correlation between average column density and presence/absence of overdensities

Q: 75% of dense cores in filaments. Must be a large number of isolated ones. Any difference between on-core and off-core?

A: Yes, see Polychroni's poster. More massive on filaments.

Q: Tan - Universal widths. Seem to be a factor of 2 from resolution for all widths. How sure?

A: The physical width is constant to within a factor of 2. Angular width varies.

Only the average width is constant

PPVI Talk 1: Mark Krumholz

2013/07/15
FORMATION OF MOLECULAR CLOUDS AND GLOBAL CONDITIONS FOR STAR FORMATION

Authors:

C. Dobbs (Exeter University, United Kingdom), M. Krumholz (University of California, Santa Cruz, United States), J. Ballesteros-Paredes (CRyA-UNAM, Mexico), A. Bolatto (University of Maryland, United States), Y. Fukui (Nagoya University, Japan), M. Heyer (University of Massachusetts, United States), M. Mac Low (American Museum of Natural History, United States), E. Ostriker (Princeton University, United States), E. Vazquez-Semadeni (UNAM, Morelia, Mexico)

Abstract

Giant molecular clouds (GMCs) are the primary reservoirs of cold, star-forming molecular gas in the Milky Way and similar galaxies, and thus any understanding of star formation must encompass a model for GMC formation, evolution, and destruction. These models are necessarily constrained by measurements of interstellar molecular and atomic gas, and the emergent, newborn stars. Both observations and theory have undergone great advances in recent years, the latter driven largely by improved numerical simulations, and the former by the advent of large-scale surveys with new telescopes and instruments. This chapter reviews the current state of the field.
Content
Galactic Scale

HI column density "hits a wall"

molecular gas correlated with SF

"depletion time" ~ 2 Gyr

power law index of molecular clouds: -1.5 to -2 -> most mass in big clouds

only in gas-rich regions

steeper in outer regions of galaxies: -2 to -2.5

typical mean surface density 100 Msun/pc^2

PDF of surface densities "lognormal" with powerlaw tail
Theory

Virial Parameter: $\sigma_v = (\alpha_G \pi G \Sigma R / 5)^{1/2}$

$\alpha_G = -2 T/ W$

$\alpha_G = 2$ -> free-fall collapse

3 options: pressure confinement, virial balance, or free-fall collapse

Most velocity dispersion on largest scale

median $\alpha_G \sim 0.1$

$M/M_{crit} = \frac{M}{\phi/(4\pi G)^{1/2}} \approx 2$

B-fields not strong enough to support against gravity in mol. clds (Crutcher)
Lifetime of Molecular Clouds

Spatial distribution + galactic dynamics.

Method 1: compare dynamics to distribution

Formed in the arms and survive passage between arms? (implies lifetime ~100 Myr)

Method 2: compare star cluster ages to spatial cloud distribution

lifetimes closer to 30 Myr

Method 3: look at T Tau vs post T Tau (Hartmann 2009)

3 Myr for Taurus

Why difference?

size matters

should compare lifetimes to cloud dynamical times

$t_{ff} = \left(\frac{3\pi}{32 G \rho}\right)^{3/2}$
Efficiency

The rate at which gas converts into stars.

$\epsilon_{ff} = SFR / (M_{gas}/t_{ff}) \sim 0.01$

if pure free-fall, $\epsilon_{ff} = 1$

Low star formation rate

Lifetime << 100 $t_{ff}$

Lifetime * efficiency < 1, so low overall star formation efficiency
Metallicity

Break between HI and CO:

10 Msun in Galaxy, local spirals

50 Msun in SMC
How do Molecular Clouds form?

"Stories", not conclusions
Converging Flows

Collision -> cooling

vazquez-semadeni et al

"local turbulence" triggers collision of warm HI

Can't make big clouds: limited to mean surface density * scale height? [needs more]
Cloud Collision in Spiral Arms

Spiral arms are convergence points

Collisions unlikely in most of galaxy

orbits crowd in arms: collisions more likely

mechanism probably can't work in flocculent galaxies

explains counter-orbiting GMCs [are these really a common feature?]
Gravitational Magneto-Jeans (Kim & Ostriker 2006)

Toomre Q instability at Q < 1.5

makes 10^7-10^8 Msun clouds; smaller ones must fragment

produces low spin clouds. Mol. Clds apparently have much less angular momentum than expected

easily explains "beads on a string" morphology [why?]
Parker + Thermal instability (Mouschovias et al 2005)

Buoyancy draws B-fields out of plane

may not work with turbulence
Chemistry: Forming H2 and CO

H2 dominance begins at ~10 Msun/pc^2

Glover and Clark: Is the molecular transition a necessary prerequisite to star formation?

Extincting FUV allows cooling, H2 formation
Evolution of Molecular Clouds

Morphology?
Highest intensity offset from highest column (DR21; Schneider 2001)

MHD collapse (Hennebelle) vs. pure turbulence (Federrath)
Non-thermal motions in GMCs:
Molecular clouds are invariably turbulent

Reynolds number ~10^9

Turbulence decays: Why sigma so large?

Global grav collapse?

Internal energy injection?

External driving?
External driving:
Accretion energy

Large-scale galactic shear

Requires feedback to get star formation right
Internal driving:
Gritschneder: ejected clumps stirring up

Someone else: B-fields distorted by ionization, store energy
Mol Cld Disruption

Dont' survive, must be disrupted....
Why SF rate low?
feedback

or turbulence
Disruption Scenarios

ionization can disrupt up to 10^5, may not work above that mass

turbulence slows down SF

turbulence must be maintained
Scales

"Normal" from 10-100 msun/pc^2 (all clouds are the same; cloud counting)

extreme above 100 msun/pc^2 (densities in mol clds very high; mean surface density of galaxy matches individual clouds)

"all hell break loose" < 10 msun/pc^2
Future
Obs: big surveys, better data

Thry: Need to combine large and small scales
Citations

Leroy 2008

Roman-Duval 2010

Gratier 2012

Narayanan 2008: Taurus cloud CO 3-color by velocity

Crutcher 201

Koda 2009 M51

Kawamura 2009 LMC

Rosolowsky & Blitz 2005
Questions
Q: Correlation with CO, H2... have there been correlations with others?

A: Yes, very well correlated with HCN, but no maps

Q Tan: Tasker showed that GMC collisions can happen on short timescales...????

Q Tan: B-fields relatively weak by "factor of 2, but don't believe to factor of

2". Could B-fields still be important?

A: I know about your work with Tasker. Requires high surface density. Can't

work in, e.g., M33.

A: B-fields. If only a few measurements, fine. But with so much data...

(systematic??)

Comment: Reconciling ages in LMC vs local. Function of measurement method: LMC, use OB stars outside. Inside clouds, locally, don't use OB. Ages may be biased...

A: What do you call a cloud? Clump?

Q: Should we stop talking about virial equilibrium, start talking about

equipartition?

A: They may be in equilibrium. No "true" equilibrium, but a "statistical

equilibrium". Equipartition is uninteresting since the energies are automatically equal.

GPS Data Analysis

2013/06/09

I've been running up the same mountain about a dozen times with a GPS on. With a little hacking, I was able to take the GPS data and normalize it to the same start & end points to see how well I did on different runs.

Hopkins PDF Generalization

2013/06/04
Hopkins 2013 presents a non-lognormal form of the probability distribution function of density in turbulence, depending primarily on a parameter $T$. I examine some of its properties here. All of the equations shown on this page are implemented in https://github.com/keflavich/turbulence_pdfs, in particular the functions moments_theoretical_hopkins and moments in hopkins_pdf.py.

The conservation equations, conservation of probability and mass

\begin{equation*} \int_{-\infty}^\infty P_V(\ln \rho) d \ln \rho \equiv 1 \end{equation*}

\begin{equation*} \int_{-\infty}^\infty\rho P_V(\ln \rho) d \ln \rho \equiv \rho_0 = 1 \end{equation*}

$\rho_0=1$ is the definition used by Hopkins 2013. Part of the goal of this document is to generalize to $\rho_0 <> 1$

The most important thing to do is determine the values of the moments of the distribution, including the integral corrections for the Dirac Delta function component of the PDF.
The PDF

The full form of the Hopkins PDF is

\begin{equation*} P_V(\ln \rho/\rho_0) d \ln \rho = \left[I_1(2\sqrt{\lambda u}) e^{-\lambda-u} \sqrt{\frac{\lambda}{u}} + e^{-\lambda} \delta(u)\right]du \end{equation*}

The mass-weighted PDF is

\begin{equation*} P_M(\ln \rho/\rho_0) d \ln \rho = \frac{\rho}{\rho_0} \left[I_1(2\sqrt{\lambda u}) e^{-\lambda-u} \sqrt{\frac{\lambda}{u}} + e^{-\lambda} \delta(u)\right]du \end{equation*}

With definitions:

\begin{equation*} u\equiv \frac{\lambda}{1+T} - \frac{\ln \rho/\rho_0}{T} ; (u \geq 0) \end{equation*}

\begin{equation*} \lambda \equiv \frac{S_{\ln \rho,V}}{2 T^2} \end{equation*}

The ratio $\rho/\rho_0$ is needed to ensure self-consistency: the mass PDF must also conserve probability, which means it must be scaled by $1/\rho_0$.

These distributions are slightly different than those shown in Hopkins 2013. The $m=0$ term in the summation version of the PDF yields a Dirac $\delta$ function in the integral form of the PDF. To quote Phil:
```
If you check, the summation expression is exact, and does exactly integrate
to unity (and properly conserve mass) for any value of T or lambda (with
u>0). But this is with the implicit understanding that the gamma-function
distribution with exponent "m" should be taken to go to a delta-function
about u=0 when m=0 (since otherwise it is ill defined; and physically this
is the "no event/multiplication" case, so that's the only allowed case for
the term in u).
```
Moment equations

The important moments are the mean and variance of the volume and mass density probability distribution and the volume and mass log-density probability distribution.

The definitions used here are $E[\rho] \equiv <\rho>_V$, etc. The variance is computed as $V[x]=E[x^2]-E[x]^2$, which means we'll need to compute the moments of various combinations of variables.

Also, $S_{\ln \rho,V} \equiv V[\ln \rho_V]$

\begin{equation*} <\rho>_V \equiv \rho_0 \equiv \int_{-\infty}^{\infty} \rho P_V(\ln \rho/\rho_0) d \ln \rho \end{equation*}

\begin{equation*} <\ln \rho>_V \equiv \int_{-\infty}^{\infty} \ln \rho P_V(\ln \rho/\rho_0) d \ln \rho \end{equation*}

\begin{equation*} <\rho>_M \equiv \int_{-\infty}^{\infty} \rho (\rho P_V(\ln \rho/\rho_0)) d \ln \rho \end{equation*}

\begin{equation*} <\ln \rho>_M \equiv \int_{-\infty}^{\infty} \ln \rho (\rho P_V(\ln \rho/\rho_0)) d \ln \rho \end{equation*}

\begin{equation*} <\rho^3>_V \equiv \int_{-\infty}^{\infty} \rho^3 P_V(\ln \rho/\rho_0) d \ln \rho \end{equation*}

\begin{equation*} <\ln^2 \rho>_V \equiv \int_{-\infty}^{\infty} (\ln \rho)^2 P_V(\ln \rho/\rho_0)d \ln \rho \end{equation*}

\begin{equation*} <\ln^2 \rho>_M \equiv \int_{-\infty}^{\infty} (\ln \rho)^2 \rho P_V(\ln \rho/\rho_0)d \ln \rho \end{equation*}

In the moment equations, only the input to the PDF is scaled by a mean density $\rho_0$; the "weighting factors" for the expectation value are not (i.e., we are measuring $E[\rho]$, not $E[\rho/rho_0]$).

Delta terms

These can be integrated analytically. In this section, I am just computing the means; the variances are more complicated.

\begin{equation*} <\rho>_V \equiv \rho_0 \equiv \int_{-\infty}^{\infty} \rho e^{-\lambda} \delta(\frac{\lambda}{1+T} - \frac{\ln\rho/\rho_0}{T}) \frac{d \ln \rho}{T} \end{equation*}

\begin{equation*} <\ln \rho>_V \equiv \int_{-\infty}^{\infty} \ln \rho e^{-\lambda} \delta(\frac{\lambda}{1+T} - \frac{\ln\rho/\rho_0}{T}) \frac{d \ln \rho}{T} \end{equation*}

\begin{equation*} <\rho>_M \equiv \int_{-\infty}^{\infty} \rho (\frac{\rho}{\rho_0} e^{-\lambda} \delta(\frac{\lambda}{1+T} - \frac{\ln\rho/\rho_0}{T})) \frac{d \ln \rho}{T} \end{equation*}

\begin{equation*} <\ln \rho>_M \equiv \int_{-\infty}^{\infty} \ln \rho (\frac{\rho}{\rho_0} e^{-\lambda} \delta(\frac{\lambda}{1+T} - \frac{\ln\rho/\rho_0}{T})) \frac{d \ln \rho}{T} \end{equation*}

Substitution: $v=\frac{\ln \rho/\rho_0}{T}$, $dv = \frac{1}{T} d \ln \rho$, $\rho=\rho_0 e^{v*T}$, $\ln \rho = v T + \ln \rho_0$

\begin{equation*} <\rho>_{V\delta} \equiv \rho_0 \equiv \int_{-\infty}^{\infty} \rho_0 e^{vT} e^{-\lambda} \delta(\frac{\lambda}{1+T} - v) d v \end{equation*}

\begin{equation*} <\ln \rho>_{V\delta} \equiv \int_{-\infty}^{\infty} (vT + \ln \rho_0) e^{-\lambda} \delta(\frac{\lambda}{1+T} - v) d v \end{equation*}

\begin{equation*} <\rho>_{M\delta} \equiv \int_{-\infty}^{\infty} \rho_0 e^{2vT} ( e^{-\lambda} \delta(\frac{\lambda}{1+T} - v)) d v \end{equation*}

\begin{equation*} <\ln \rho>_{M\delta} \equiv \int_{-\infty}^{\infty} (vT + \ln \rho_0) e^{vT} ( e^{-\lambda} \delta(\frac{\lambda}{1+T} - v)) d v \end{equation*}

Solutions:

\begin{equation*} <\rho>_{V\delta} = \rho_0 \exp\left[\frac{T \lambda }{1+T} - \lambda\right] = \rho_0 \exp\left[-\lambda \frac{1}{1+T}\right] \end{equation*}

\begin{equation*} <\ln \rho>_{V\delta} = e^{-\lambda} \frac{\lambda T}{1+T} + e^{-\lambda} \ln \rho_0 \end{equation*}

\begin{equation*} <\rho>_{M\delta} = \rho_0 \exp\left[\frac{2 T \lambda }{1+T} - \lambda\right] = \rho_0 \exp\left[\lambda\frac{T-1}{T+1}\right] \end{equation*}

\begin{equation*} <\ln \rho>_{M\delta} = \left( \frac{\lambda T}{1+T} + \ln \rho_0 \right) \exp\left[\frac{T \lambda }{1+T} - \lambda\right] \end{equation*}

\begin{equation*} = \left( \frac{\lambda T}{1+T} + \ln \rho_0 \right) \exp\left[\frac{ -\lambda }{T+1}\right] \end{equation*}

(for these next 3, I skipped intermediate steps)

\begin{equation*} <\rho^3>_{V\delta} = \rho_0^2 \exp\left[\frac{3 T \lambda }{1+T} - \lambda\right] = \rho_0^2 \exp\left[\lambda\frac{2T-1}{T+1}\right] \end{equation*}

\begin{equation*} <\ln^2 \rho>_{M\delta} = \left( \frac{\lambda T}{1+T} + \ln \rho_0 \right)^2 \exp\left[\frac{ -\lambda }{T+1}\right] \end{equation*}

\begin{equation*} <\ln^2 \rho>_{V\delta} = \left( \frac{\lambda T}{1+T} + \ln \rho_0 \right)^2 e^{-\lambda} \end{equation*}

Using $\rho_0=1$ as defined in Hopkins 2013 simplifies all of these a great deal.

PDF Integrals

These cannot be integrated analytically.

However, we can work from a few simple mathematica/sympy results:

\begin{equation*} \int_0^\infty I_1(x) e^{-x^2/(4L)} dx = e^L - 1 \end{equation*}

\begin{equation*} \int_0^\infty x^2 I_1(x) e^{-x^2/(4L)} dx = 4 L^2 * e^L \end{equation*}

\begin{equation*} \int_0^\infty x^4 I_1(x) e^{-x^2/(4L)} dx = 16 L^3 (L+2) * e^L \end{equation*}

We use $L$ instead of $\lambda$ in these equations because it is often substituted in later equations.

Expectation Value of the Volume-Weighted Density $E[\rho]$

\begin{equation*} E[\rho] \equiv \int \rho P_v(\ln \rho/\rho_0) d \ln \rho = \int \rho \left[I_1(2\sqrt{\lambda u}) e^{-\lambda-u} \sqrt{\frac{\lambda}{u}} + e^{-\lambda} \delta(u)\right]du \end{equation*}

To get to the form of the above equations, we use the substitution

\begin{equation*} x = 2\sqrt{\lambda u} \end{equation*}

which gives us $\rho$ in terms of $x$:

\begin{equation*} \rho = \rho_0 \exp\left[T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right)\right] \end{equation*}

and leads to the rearrangement:

\begin{equation*} E[\rho] = \int \rho_0 \exp\left[T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right)\right] \left[I_1(x) e^{-x^2/(4\lambda)-\lambda} \right]dx + \rho_0 \exp\left(- \frac{\lambda}{1+T}\right) \end{equation*}

where the rightmost term is kept from the first moment above. The integral term can straightforwardly be broken apart into equations of the form shown above.

\begin{equation*} L \rightarrow \frac{\lambda}{1+T} \end{equation*}

\begin{equation*} E[\rho] = \rho_0 \left[ \exp \left(-\lambda+\frac{T\lambda}{1+T}\right) \int \left[I_1(x) e^{-x^2/(4L)} \right]dx +\exp\left(- \frac{\lambda}{1+T}\right) \right] \end{equation*}

\begin{equation*} = \rho_0 \left[ \exp \left(-\lambda+\frac{T\lambda}{1+T}\right)(e^L-1) +\exp\left(- \frac{\lambda}{1+T}\right) \right] \end{equation*}

\begin{equation*} = \rho_0 \left[ \exp \left(-\lambda+\frac{T\lambda}{1+T}\right)(e^{\lambda/1+T}-1) +\exp\left(- \frac{\lambda}{1+T}\right) \right] \end{equation*}

\begin{equation*} = \rho_0 \left[ e^{-\lambda/(1+T)}(e^{\lambda/1+T}-1) +\exp\left(- \frac{\lambda}{1+T}\right) \right] \end{equation*}

\begin{equation*} = \rho_0 \end{equation*}

The same general approach can be followed for all expectation values, but we'll skip the detailed algebra.

Variance of the Volume-Weighted Density $V[\rho]=S_{\ln \rho,V}$

\begin{equation*} V[\rho] = E[\rho^2] - E[\rho]^2 = \rho_0^2 \left[ \exp\left(\lambda\frac{2 T^2}{1+3T+2T^2}\right) - 1 \right] \end{equation*}

However, the "correction factor" is still important:

\begin{equation*} V_\delta[\rho] = \rho_0^2 \left[ \exp\left(\lambda\frac{T-1}{T+1}\right) - \exp\left(-2\frac{\lambda}{1+T}\right) \right] \end{equation*}

Expectation Value of the Mass-Weighted Density $E_M[\rho]$

Start from halfway through $E[\rho]$, simply adding a factor of 2 in the exponent:

\begin{equation*} E_M[\rho] = \int \rho_0 \exp\left[2T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right)\right] \left[I_1(x) e^{-x^2/(4\lambda)-\lambda} \right]dx + \rho_0 \exp\left(- \frac{\lambda(T-1)}{1+T}\right) \end{equation*}

Follow the same math, with $L=\frac{\lambda}{1+2T}$

\begin{equation*} = \rho_0 \left[ \exp \left(-\lambda+\frac{2T\lambda}{1+T}\right)(e^L-1) +\exp\left(- \frac{\lambda(T-1)}{1+T}\right) \right] \end{equation*}

\begin{equation*} = \rho_0 \left[ \exp \left(\frac{(T-1)\lambda}{1+T}\right)(e^{\lambda/(1+2T)}-1) +\exp\left(- \frac{\lambda}{1+T}\right) \right] \end{equation*}

\begin{equation*} E_M[\rho] = \rho_0 \left[ \exp\left(\lambda\frac{2 T^2}{1+3T+2T^2}\right) - \exp\left(\lambda\frac{T-1}{T+1}\right) + \exp\left(\lambda\frac{T-1}{T+1}\right) \right] \end{equation*}

The right 2 terms cancel, yielding the value shown in Equation 7 of Hopkins 2013 scaled by $\rho_0^2$. However, the right-most term is the correction factor from the Dirac Delta term needed to correct any numerical computation of the mass-weighted density.

\begin{equation*} E_{\delta,M}[\rho] = \exp\left(\lambda\frac{T-1}{T+1}\right) \end{equation*}

\begin{equation*} E_M[\rho] = \rho_0 \exp\left(\lambda\frac{2 T^2}{1+3T+2T^2}\right) \end{equation*}

Expectation Value of the Mass-Weighted Density Squared $E_M[\rho^2]$

\begin{equation*} E_M[\rho^2] = \int \rho^2 \frac{\rho}{\rho_0} e^{-\lambda} \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx + \int \rho^2 \frac{\rho}{\rho_0} e^{-\lambda} \delta(u) du \end{equation*}

\begin{equation*} E_{\delta,M}[\rho^2] = \rho_0^2 \exp\left[\lambda\frac{2T-1}{T+1}\right] \end{equation*}

\begin{equation*} E_{left}[\rho^2] = e^{-\lambda} \int \rho_0^2 \exp\left[3T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right)\right] \left[I_1(x) e^{-x^2/(4\lambda)} \right] dx \end{equation*}

\begin{equation*} = \rho_0^2 \exp\left[\frac{(2T-1)\lambda}{1+T}\right] \int I_1(x) e^{-(3T+1)x^2/(4\lambda)} dx \end{equation*}

\begin{equation*} = \rho_0^2 \exp\left[\frac{(2T-1)\lambda}{1+T}\right] \left( \exp\left[\frac{\lambda}{3T+1}\right] - 1\right) \end{equation*}

\begin{equation*} = \rho_0^2 \left(\exp\left[\frac{6\lambda T^2}{3T^2+4T+1}\right] - \exp\left[\frac{(2T-1)\lambda}{1+T}\right] \right) \end{equation*}

\begin{equation*} E_{M}[\rho^2] = \rho_0^2 \exp\left[\frac{6\lambda T^2}{3T^2+4T+1}\right] \end{equation*}

Variance of the Mass-Weighted Density $V_M[\rho] = E_M[\rho^2] - E_M[\rho]^2$

Since correction factors are given for $E_M[\rho^2]$ and $E_M[\rho]$, they are not included separately here:

\begin{equation*} V_M[\rho] = E_M[\rho^2] - E_M[\rho]^2 = \rho_0^2 \left( \exp\left[\frac{6\lambda T^2}{3T^2+4T+1}\right] -\exp\left[\lambda\frac{4 T^2}{1+3T+2T^2}\right] \right) \end{equation*}

Expectation Value of the Volume-Weighted Log Density $E[\ln \rho]$

\begin{equation*} E[\ln \rho] = \int \ln \rho e^{-\lambda} \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx + \int \ln \rho e^{-\lambda} \delta(u) du \end{equation*}

\begin{equation*} E_\delta[\ln \rho] = e^{-\lambda} \left[ \frac{\lambda T}{1+T} + \ln \rho_0 \right] \end{equation*}

\begin{equation*} E_{left}[\ln \rho] = \int \left[\ln \rho_0 + T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right) \right] e^{-\lambda} \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx \end{equation*}

\begin{equation*} = e^{-\lambda} \left( \int \left[\ln \rho_0 + \frac{T\lambda}{1+T}\right] \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx - \int \frac{T x^2}{4\lambda} \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx \right) \end{equation*}

\begin{equation*} = e^{-\lambda} \left( \int \left[\ln \rho_0 + \frac{T\lambda}{1+T}\right](e^{\lambda}-1) - \frac{4 T \lambda^2 e^{\lambda}}{4\lambda} \right) \end{equation*}

\begin{equation*} = \left( \left[\ln \rho_0 + \frac{T\lambda}{1+T}\right](1-e^{-\lambda}) - T \lambda \right) \end{equation*}

\begin{equation*} E[\ln \rho] = \ln \rho_0 + \frac{T\lambda}{1+T} - T \lambda \end{equation*}

\begin{equation*} = \ln \rho_0 - \frac{T^2\lambda}{1+T} \end{equation*}

Expectation Value of the Mass-Weighted Log Density $E_M[\ln \rho]$

\begin{equation*} E_M[\ln \rho] = \int \ln \rho \frac{\rho}{\rho_0} e^{-\lambda} \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx + \int \ln \rho \frac{\rho}{\rho_0} e^{-\lambda} \delta(u) du \end{equation*}

\begin{equation*} E_\delta[\ln \rho] = \left( \frac{\lambda T}{1+T} + \ln \rho_0 \right) \exp\left[\frac{ -\lambda }{T+1}\right] \end{equation*}

\begin{equation*} E_{left}[\ln \rho] = e^{-\lambda} \int \left[ \left(\ln \rho_0 + T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right) \right) \exp\left(T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right)\right) \right] \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx \end{equation*}

Again, separate into integrable terms:

\begin{equation*} = \exp\left(\frac{T\lambda}{1+T} -\lambda\right) \left[ \left(\ln \rho_0 + \frac{T\lambda}{1+T} \right) \left[I_1(x) e^{-(1+T)x^2/(4\lambda)} \right] + \left(-\frac{Tx^2}{4\lambda} \right) \left[I_1(x) e^{-(1+T)x^2/(4\lambda)} \right] \right] \end{equation*}

\begin{equation*} L = \frac{\lambda}{1+T} \end{equation*}

\begin{equation*} E_{left}[\ln \rho] = \exp\left(\frac{T\lambda}{1+T} -\lambda\right) \left[ \left(\ln \rho_0 + \frac{T\lambda}{1+T} \right) \left(\exp\left[\frac{\lambda}{1+T}\right]-1\right) + \left(-\frac{T}{4\lambda} \right) \left(\frac{4 \lambda^2}{(1+T)^2} \exp\left[\frac{\lambda}{1+T}\right]\right) \right] \end{equation*}

\begin{equation*} E_{left}[\ln \rho] = \exp\left(\frac{-\lambda}{1+T}\right) \left[ \left(\ln \rho_0 + \frac{T\lambda}{1+T} \right) \left(\exp\left[\frac{\lambda}{1+T}\right]-1\right) + \left(-\frac{T}{4\lambda} \right) \left(\frac{4 \lambda^2}{(1+T)^2} \exp\left[\frac{\lambda}{1+T}\right]\right) \right] \end{equation*}

\begin{equation*} E_{left}[\ln \rho] = \left(\ln \rho_0 + \frac{T\lambda}{1+T} \right) \left(1-\exp\left[\frac{-\lambda}{1+T}\right]\right) - \left(\frac{ T \lambda}{(1+T)^2} \right) \end{equation*}

\begin{equation*} E_M[\ln \rho] = \left(\ln \rho_0 + \frac{T\lambda}{1+T} \right) - \left(\frac{ T \lambda}{(1+T)^2} \right) \end{equation*}

\begin{equation*} = \ln \rho_0 + \frac{T^2\lambda}{(1+T)^2} \end{equation*}

Expectation Value of the Mass-Weighted Log Density Squared $E_M[\ln^2 \rho]$

\begin{equation*} E_M[\ln^2 \rho] = \int (\ln \rho)^2 \frac{\rho}{\rho_0} e^{-\lambda} \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx + \int (\ln \rho)^2 \frac{\rho}{\rho_0} e^{-\lambda} \delta(u) du \end{equation*}

\begin{equation*} E_\delta[\ln^2 \rho] = \left( \frac{\lambda T}{1+T} + \ln \rho_0 \right)^2 \exp\left[\frac{ -\lambda }{T+1}\right] \end{equation*}

\begin{equation*} E_{left}[\ln^2 \rho] = e^{-\lambda} \int \left[ \left(\ln \rho_0 + T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right) \right)^2 \exp\left(T\left(-\frac{x^2}{4\lambda} + \frac{\lambda}{1+T}\right)\right) \right] \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx \end{equation*}

This time it's just too ugly. Define a new variable:

\begin{equation*} Q = \ln \rho_0 + \frac{T\lambda}{1+T} \end{equation*}

\begin{equation*} E_{left}[\ln^2 \rho] = e^{-\lambda/(1+T)} \int \left[ \left(Q -\frac{T x^2}{4\lambda} \right)^2 \exp\left(-\frac{Tx^2}{4\lambda} \right) \right] \left[I_1(x) e^{-x^2/(4\lambda)} \right]dx \end{equation*}

\begin{equation*} E_{left}[\ln^2 \rho] = e^{-\lambda/(1+T)} \int \left[ \left(Q^2 - 2 Q \frac{T x^2}{4\lambda} + \frac{T^2 x^4}{16\lambda^2} \right) \left[I_1(x) e^{-(1+T)x^2/(4\lambda)} \right] \right]dx \end{equation*}

\begin{equation*} E_{left}[\ln^2 \rho] = e^{-\lambda/(1+T)} \left[ Q^2 (e^{\lambda/(1+T)}-1) - 2 Q \frac{T}{4\lambda} \left(\frac{4\lambda^2}{(1+T)^2} e^{\lambda/(1+T)}\right) + \frac{T^2}{16\lambda^2} \left(16 \frac{\lambda^3}{(1+T)^3} (\frac{\lambda}{1+T}+2) e^{\lambda/(1+T)} \right) \right] \end{equation*}

\begin{equation*} E_{left}[\ln^2 \rho] = Q^2 (1-e^{-\lambda/(1+T)}) - 2 Q \frac{T\lambda}{(1+T)^2} + \frac{\lambda T^2}{(1+T)^3} \left(\frac{\lambda}{1+T}+2\right) \end{equation*}

Add back the $\delta$ termG

\begin{equation*} E_M[\ln^2 \rho] = Q^2 - 2 Q \frac{T\lambda}{(1+T)^2} + \frac{\lambda T^2}{(1+T)^3} \left(\frac{\lambda}{1+T}+2\right) \end{equation*}

Re-expand to see if it's simplifiable....

\begin{equation*} E_M[\ln^2 \rho] = \left(\ln \rho_0 + \frac{T\lambda}{1+T}\right)^2 - 2 \left(\ln \rho_0 + \frac{T\lambda}{1+T}\right) \frac{T\lambda}{(1+T)^2} + \frac{\lambda T^2}{(1+T)^3} \left(\frac{\lambda}{1+T}+2\right) \end{equation*}

\begin{equation*} E_M[\ln^2 \rho] = \ln^2 \rho_0 + 2 \ln \rho_0 \frac{T\lambda}{1+T} + \frac{T^2\lambda^2}{(1+T)^2} - 2 \ln \rho_0 \frac{T\lambda}{(1+T)^2} - 2 \frac{T\lambda}{1+T} \frac{T\lambda}{(1+T)^2} + 2 \frac{\lambda T^2}{(1+T)^3} + \frac{\lambda T^2}{(1+T)^3} \frac{\lambda}{1+T} \end{equation*}

\begin{equation*} E_M[\ln^2 \rho] = \ln^2 \rho_0 + 2 \ln \rho_0 \left(\frac{T\lambda}{1+T} - \frac{T\lambda}{(1+T)^2}\right) + \frac{T^2\lambda^2(1+T)}{(1+T)^3} - 2 \frac{T^2\lambda^2}{(1+T)^3} + 2 \frac{\lambda T^2}{(1+T)^3} + \frac{\lambda^2 T^2}{(1+T)^4} \end{equation*}

\begin{equation*} E_M[\ln^2 \rho] = \ln^2 \rho_0 + 2 \ln \rho_0 \left(\frac{T^2\lambda}{1+T}\right) + \frac{T^2\lambda^2(1+T)}{(1+T)^3} - 2 \frac{T^2\lambda^2}{(1+T)^3} + 2 \frac{\lambda T^2}{(1+T)^3} + \frac{\lambda^2 T^2}{(1+T)^4} \end{equation*}

\begin{equation*} = \ln^2 \rho_0 + 2 \ln \rho_0 \frac{T^2 \lambda}{(1+T)^2} + \left(\frac{\lambda T^2}{(1+T)^2}\right)^2 + \frac{2\lambda T^2}{(1+T)^3} \end{equation*}

\begin{equation*} = \left(\ln \rho_0 + \frac{T^2\lambda}{(1+T)^2}\right)^2 + \frac{2\lambda T^2}{(1+T)^3} \end{equation*}

As expected, we recover the correct relation from Hopkins 2013:

\begin{equation*} S_{\ln \rho,M} = E_M[\ln \rho^2] - E_M[\ln \rho]^2 = \ln^2 \rho_0 + 2 \ln \rho_0 \frac{T^2 \lambda}{(1+T)^2} + \left(\frac{\lambda T^2}{(1+T)^2}\right)^2 + \frac{2\lambda T^2}{(1+T)^3} - \left(\ln \rho_0 + \frac{T^2\lambda}{(1+T)^2}\right)^2 \end{equation*}

\begin{equation*} = \frac{2\lambda T^2}{(1+T)^3} \end{equation*}

independent of $\rho_0$.

Properties of the Hopkins PDF

I began this investigation in order to find out whether a different parameter, i.e. something related to the compressiveness of the turbulent driving, could be responsible for the "discrepancy" between the Formaldehyde-derived density and the volume-averaged density of some clouds.

I naively expected that, in compressive turbulence, more of the mass will be concentrated at higher density, which should drive up the mass-weighted mean density.

Hopkins 2013 showed that $T\sim\mathcal{M}_C$. If this is taken on face value, without recognizing the relation between $T$ and $\sigma$, it means that for fixed $\sigma$, $<\rho>_M\propto e^{-T^2}$.

This figure shows the relation of the various moments and $T$. The plots show both the "theoretical" result (i.e., doing the integral by hand) and the numerical result with 50,000 data points plus a correction for the $\delta$ terms; the agreement is essentially perfect excepting some numerical noise (the only visible discrepancy is for a perfect lognormal at high $\sigma_V$):

However, Hopkins 2013 also found that simulations generally produced $T\sim0.1\sigma_{s,M}^2$ (though a relation of the form $T\sim 0.25 \ln(1+0.25\sigma_{s,M}^4)$ is also a good fit).

If you use the simpler of these relations, the expected relationship (mean mass increasing with increasing "compressiveness") is recovered. However, it comes with an increasing $\sigma_s$. In these plots, half are labeled with $\sigma$ and the others are labeled with $T$, though the two are equivalent. The highest $\sigma_s$ observed in any of the simulations shown in Hopkins 2013 was about $\sigma_s=4$ (marked with a black square below), so the maximum $\rho_M/\rho_V \sim 10^2$.

BGPS Point Source Recovery - Experiment 23 results

2013/05/20

Experiment 23, run in mid-May 2013, created power-law synthetic astrophysical sky maps with added point sources in order to examine the pipeline + bolocat's ability to recover point source flux densities accurately.

For the purpose of the BGPS v2 paper, this data is incorporated into Section 5.3: the angular transfer function and comparison with other data sets.

The parameter space explored includes 3 steps in power-spectrum background, $\alpha_{ps}=1,1.5,2$, which is the major parameter being explored. A background with $\alpha_{ps}=2$ is approximately what we measure in both HiGal and the reproduced regions of the BGPS power-spectrum, so it's probably the most realistic, but it is also very bad for the pipeline: most of the power is on the largest angular scales. $\alpha_{ps}=1$, on the other hand, beautifully creates point-source-like structures.

These two images show the results of $\alpha_{ps}=2$ and $\alpha_{ps}=1$ skies respectively with faint point sources on top. The point sources are more pointy and stick out more in the $\alpha_{ps}=1$ situation, but the $\alpha_{ps}=1$ map also looks more generally point-like, so there is greater confusion. If you believe that there are genuine point sources in the 1.1 mm maps, the $\alpha_{ps}=1$ would actually be very difficult to deal with since many extracted sources would actually be part of the background.

We also explored two different sets of source brightnesses, 0.1-1 and 1-10 Jy/beam, uniformly distributed. We picked a high source density, 500 in a 512x512 pixel map, to get decent source extraction statistics - the crowding is still fairly low and does not significantly affect the extracted source properties (only a few source overlap with others). The bright sources are obviously better-recovered than the faint, with good recovery for all of the backgrounds explored. The faint sources, on the other hand, were fairly sensitive to the background.

The background levels used were $\sim1,2,10$ Jy/beam at the peak for the $\alpha=2$ maps (they were lower for the other power-laws, so those will be ignored from now on). As described in the paper, the faint source recovery was good for the low peak backgrounds, but recovery was essentially nonexistent for the 10 Jy/beam background - point sources were not detected at all.

Despite the relative simplicity of this experiment, the data took 51 GB of storage and took about half a day to run.

In principle, one would like to examine a range of different source distributions (power-law flux distribution, upper/lower limits, sizes) on a range of different power-spectrum backgrounds - for the purpose of the v2 paper, this approach would be thoroughly excessive. However, I expect Tim will be taking this sort of approach for the next paper.

BGPS Point Source Recovery

2013/05/16
In the last couple drafts of the BGPS paper (e.g. the May 15, 2013 draft), we've included a single simulation demonstrating the flux recovery of the Bolocat:

These figures are regarded as very important, as they demonstrate the capability of recovering accurate flux density measurements from a realistic sky using bolocat.

However, the above figures show a simulation only for one set of parameters, with $\alpha=1$ in the power-law flux distribution, which unfortunately is not realistic according to "Figure 8":

So I've started up a new experiment, experiment #23, to examine this problem.

The problem has a few layers:
1. The reason I used $\alpha=1$ is that it looks much like a realistic BGPS map after processing, in the sense that most of the field is empty but there are a few hundred sources in the map. However, $\alpha=1$ is not a realistic representation of the measured power spectra.
2. $\alpha=2$ maps with the previous normalization had a peak value of 18 Jy, which resulted in heavily signal-dominated output maps that did not resemble BGPS maps.
3. The normalization is tricky. One of the key goals of the simulations was to test the effect of different atmospheric to astrophysical signal ratios on the angular transfer function; in order to accomplish this, it was necessary to scale the atmospheric power based on the astrophysical power at its peak in fourier space. i.e., in real timestreams, we can measure the astrophysical to atmospheric power ratio, but we have to perform that measurement somewhere that the angular transfer function is known to be reliable. This is done at about 1 Hz.
4. The normalization is important because of the noise level. In the simulations, we use a fixed noise level of about 30 mJy in the timestreams to match our best observations (though it is not difficult to scale this to other levels). This fixed noise level means that, for some normalizations, all pixels are statistically significant. Also, even though the noise level is fixed, it will be higher because of intrinsic noise in a power-law distributed map.
5. The normalizations used in experiment #21, the angular transfer function measurement, were selected such that there would be high signal-to-noise at all angular scales. This means that white noise would not be dominant on any angular scale, since white noise is equivalent to $\alpha=0$. So it wasn't crazy to use these ridiculously high-flux maps, but it is not feasible to use the same maps for analysis of point sources. In maps for which we're interested in small-angular-scale features (<100"), we want the maps to be primarily noise-dominated with a handful of bright features either caused by adding point sources directly or from the local peaks in the power-law distributed flux.
Some notes along the way:
- Using a power-law background, the point-source sensitivity is much worse than without a power-law background. This is intuitive: a 100 mJy source on a 200 mJy background (which may easily include power fluctuations on the smallest scales of comparable magnitude) is not going to be recovered.
- Doubly important: a 1 Jy source should have peak amplitude 1 Jy, but the current method of adding point sources adds them as delta functions that are later convolved, conserving the total flux rather than the peak flux. This needs to be changed! (has been now)
Here are some examples of what the before/after look like with point sources added. The first has bright sources, the second faint sources:

With these new figures, the 40" apertures work fine, but the 120" apertures are still utterly junk. This does not make sense.

A careful analysis of a single source shows that something is wrong. Here are some annular extractions followed by the image:
```
Input Map:
reg sum npix    mean    median  min         max     var         stddev      rms
--- --- ----    ----    ------  ---         ---     ---         ------      ---
1   26.3323 22  1.19692 1.18083 1.10787     1.28825 0.00208826  0.0456975   1.1978
2   77.8553 74  1.0521  1.047   1.00426     1.123   0.000929507 0.0304878   1.05254
3   124.868 122 1.02351 1.02869 0.996566    1.04427 0.000260295 0.0161337   1.02363

Output Map:
reg sum         npix    mean        median      min         max         var         stddev      rms
--- ---         ----    ----        ------      ---         ---         ---         ------      ---
1   3.89157     23      0.169199    0.175204    0.086872    0.255151    0.00206254  0.0454152   0.175188
2   2.06843     74      0.0279517   0.0275116   -0.0695484  0.155258    0.00210834  0.0459167   0.0537554
3   0.502601    123     0.00408619  0.00629906  -0.121023   0.0834974   0.00143969  0.0379432   0.0381626

Backgrounds:
Input Map:
reg  sum      npix  mean      median    min       max      var          stddev     rms
---  ---      ----  ----      ------    ---       ---      ---          ------     ---
1    297.054  291   1.0208    1.02293   0.98094   1.05234  0.000313419  0.0177037  1.02096
3    2538.6   2618  0.969671  0.972413  0.859551  1.12495  0.0015557    0.0394423  0.970473

Output Map:
reg  sum      npix  mean        median      min        max       var         stddev     rms
---  ---      ----  ----        ------      ---        ---       ---         ------     ---
1    1.49461  291   0.00513613  0.00729431  -0.121023  0.133141  0.001586    0.0398247  0.0401545
3    5.83372  2618  0.00222831  0.00194597  -0.195075  0.181155  0.00211747  0.046016   0.0460699

Bolocat for this source:
In [200]: fields
Out[200]: ['FLUX_40', 'FLUX_40_NOBG', 'BG_40', 'FLUX_120', 'FLUX_120_NOBG', 'BG_120']

In [198]: [inp[61][f] for f in fields]
Out[198]: [0.1896538, 1.3536235, 0.38071653, 0.83893013, 10.77737, 0.42352486]

In [199]: [m20[61][f] for f in fields]
Out[199]: [0.18015364, 0.18674377, -0.016305592, 0.2868295, 0.29759517, -0.00027596406]
```
With background apertures:

However, note that the background are computed using the mmm.pro sky-background estimation procedure over a range $2r$ to $4r$ (i.e., 40-80" and 120-240" radius for the 40" and 120" diameter apertures).

The numbers shown by ds9 disagree fairly severely with those from bolocat. In particular, it appears that the background estimate returned by mmm.pro is off by a factor of 2, in this case giving 0.42 instead of 0.97. Turns out this was due to an indexing error that did not affect the pipeline results in any way.

Out of date analysis: Bolocat's flux total in the 120" aperture is 10.77 Jy/beam, background is 0.42 Jy/beam. There are 218 pixels. The resulting flux should be (10.77-0.42*218/23.8), but this gives 6.9 instead of the expected 0.84. Why?

If we do the same with the ds9 numbers, we get a total of 9.62 Jy/beam, background 0.97 Jy/beam average, so: (9.62 - 0.97*218/23.8) = 0.74. This is consistent with bolocat, and very very wrong.

If we take our background to be 1.02 instead of 0.97, we get 0.28 Jy/beam, which is exactly the right answer according to the pipeline. 1.02 comes from taking a much more local background subtraction from r=40 to r=80 arcsec, which isn't really acceptable. If we go from r=60 to r=120, the disagreement remains fairly bad, with f=0.38 Jy/beam, but certainly a lot better.
```
Bolocat after the correction:
In [297]: [m20[61][f] for f in fields]
Out[297]: [0.18015364, 0.18674377, 0.00578439, 0.25477698, 0.29759517, 0.0041758944]

In [298]: [inp[61][f] for f in fields]
Out[298]: [0.1896538, 1.3536235, 1.0216572, 0.40137589, 10.77737, 1.0119308]
```
This may mean that we'll need to re-do aperture extraction with a tighter background region everywhere. I made the change to object_photometry.

However, even with the change, even in the best case, FLUX_120 appears to be totally unreliable. FLUX_80 is acceptable with the change, but only for bright sources (for faint sources, <1 Jy, there is no recovery at all - I think this must be an issue of the source brightness still not being calculated correctly):

I'll have to continue this analysis tomorrow once the full suite of simulations has completed, but I strongly suspect that we'll have to recommend strictly against using FLUX_120 if the background is expected to be $\alpha=2$ distributed.
Day 2

The simulations have partly completed. After correcting the error with convolved point sources vs. delta functions, I reset the flux distributions to be 0.05 to 1.0 Jy for the "faint" distribution and 0.1 to 50.0 Jy for the "bright" distribution (both with power-law distributions $\alpha=2$).

The results can be summarized:
In the presence of complex background, i.e power-law distributed flux density with amplitude comparable to the detected source brightness, only the smallest aperture (40") is reliable (i.e., recovers a flux in the input and processed map).

The recovered flux density has a very high dispersion in the presence of high-amplitude power-law flux. "Very high" = $\sigma \gtrsim 1$ around a mean of 1.

There is no evident dependence of the maps on the atmospheric properties. Therefore, there's no sense in varying the atmospheric properties in the simulations.
I decided the simulations needed changing again. First, there was excessive sampling in astro/atmo parameter space; this is not needed (see point 3). More important is the power-law distribution map's peak flux. Also, it is more important to get decent sampling of bright-ish sources than to have a "physically accurate" point source distribution; the distribution of sources does not affect the recovery, but it does affect the signal-to-noise of the recovery measurement. Please don't ask me to defend this statement; it would require another 10 hours of computer time + me time that I really don't want to allocate, but I'm confident that it is true.

The essential conclusion is that, for $\alpha=2$, point source recovery is only possible if the point source is brighter than the background, which is a very intuitive result. Background annulus subtraction isn't very effective at pulling out sources.

Conclusions:
These experiments show that source recovery is very poor in the presence of a bright power-law background: it is not possible to reliably extract point sources from a map filled with power-law distributed emission brighter than or comparable to the point sources.

The 120" apertures aren't really good for anything.

There is so much source-extraction parameter space out there that any further study would really deserve its own paper.
Additional plots of interest:

40" aperture in the presence of a bright background with faint sources:

Versus the same with a faint background:

Compare these to the 120" equivalents (bright then faint background):

It's fairly easy to see why there are issues with the bright background and the 120" apertures. In this image, bright background on the left, faint background on the right, with faint sources (0.1-1 Jy).

It's more helpful to look at that previous image with the source contours superposed. These images really give a nice feel for what it means to have point sources subsumed in $\alpha=2$ background.

Recovering from "Startup drive out of memory"

2013/05/13
To continue frozen processes, you need to find all sleeping procs and restart them with kill -CONT.

In principle, this command should do that:

ps -A -o tpgid,uid,pid,state,flags | awk '$1!="0" && $2==501 && $4=="S" && $5!="4000410c" {print "kill -CONT " $3}' | bash

Without the flags exclusion, this failed for me on apple's VNC server: ps -A -o tpgid,uid,pid,state,flags | awk '$2==501 && $4=="S" && $5!="4000410c" {print "kill -CONT " $3}' | bash

This is necessary for recovering mac processes: -A Display information about other users' processes, including those without controlling terminals.

The outputs from the ps process are:
- tpgid control terminal process group ID
- pid process ID
- uid effective user ID

Bolocam Mass Calculations

2013/03/25

Adopt an effective Bolocam center frequency, f = 271.1 GHz

Use Ossenkopf, V. & Henning, Th. 1994, A&A, 291, 943 Table 1. [http://adsabs.harvard.edu/abs/1994A%26A...291..943O] Interpolate last two entries (lambda = 1.0 mm and 1.3 mm to f)

This gives $\kappa = 0.01166 cm^2 g^{-1}$ (or round to 0.0117) for gas-to-dust = 100. But the slope between these 2 points is only 1.61!

The points in OH94 table 1 jump around and locally give a variety of slopes around 1.8. A fit to kappa longward of 100 microns gives a slope of 1.8.

Interpolating kappa from 1.00 mm where OH5 gives $\kappa(1.00) = 0.0137$ (for a gas-to-dust ratio = 100) to 271.1 GHz (= 1.1058 mm) using a slope of 1.8 gives:

$\kappa = 0.01143 cm^{2} g^{-1}$

This is the same as the value used by Enoch, so let's adopt it (round to 0.0114), but note that this was calculated for 271.1 GHz.

For the v1.0 data, the BGPS has the following properties:

effective beam FWHM = 33"

gaussian beam $\sigma = 14.01$" (= FWHM / $\sqrt{(8*\ln(2))}$)

effective beam radius = 19.82" (= FWHM / $\sqrt{(4*\ln(2))}$)

This is the radius of a circular top-hat function that has the same solid angle as a gaussian. At a distance of D = 1 kpc, this radius corresponds to a length = $2.97 x 10^{17}$ cm. The effective beam solid angle is $\Omega = 2.902 x 10^{-8}$ sr.

Adopt a fiducial dust temperature for mass calculations, $T_d = 20K$

The total mass per Jy at $D_{kpc}$ = 1 kpc is given by

\begin{equation*} M = 1.0 x 10^{-23} * S(Jy) * D^2 / \kappa * B_{\nu}(T) \end{equation*}

\begin{equation*} = 14.30 * [e^{13.01/T(K) - 1}] * S(Jy) * D^2_{kpc} M_{\odot} \end{equation*}

For $T = 20 K$:

\begin{equation*} M = 13.07 * S(Jy) * D^2_{kpc} M_{\odot} \end{equation*}

This corresponds to a beam-filling column density of $H_2$,

\begin{equation*} N(H2) = 2.19e+22 * [e^{13.01/T(K) - 1}] * S(Jy) (cm^{-2}) \end{equation*}

For $T = 20 K$:

\begin{equation*} N(H2) = 2.01 x 10^{22} * S(Jy) (cm^{-2}) A_V = 21.1 \mathrm{magnitudes} Jy^{-1} \end{equation*}

where we use the Bohlin & Savage calibration

\begin{equation*} A_V = 1 mag. <=> N(H) = 1.9 x 10^{20} cm^{-2} \mathrm{and} N(H_2) = N(H) / 2 \end{equation*}

Arecibo Mapping: Total Power Measurements

2013/03/22
Worked on reducing the Arecibo continuum data. There were some serious problems with the off positions previously, in particular generating weird slopes.

I suspect the problem was with these lines
```
sorted1 = sort(totals[0:-1:2])
sorted2 = sort(totals[1:-1:2])
```
since the indexing is wrong, but I'm not sure yet. I re-wrote the whole off-position measurement to average over all scans. Still in progress, but the fitting works now. Here's an example: