Hopkins PDF Generalization

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

∞∫ − ∞PV(lnρ)dlnρ ≡ 1
∞∫ − ∞ρPV(lnρ)dlnρ ≡ ρ0 = 1

ρ0 = 1 is the definition used by Hopkins 2013. Part of the goal of this document is to generalize to ρ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

PV(lnρ ⁄ ρ0)dlnρ = [I1(2√(λu))e − λ − u√((λ)/(u)) + e − λδ(u)]du

The mass-weighted PDF is

PM(lnρ ⁄ ρ0)dlnρ = (ρ)/(ρ0)[I1(2√(λu))e − λ − u√((λ)/(u)) + e − λδ(u)]du

With definitions:

u ≡ (λ)/(1 + T) − (lnρ ⁄ ρ0)/(T);(u ≥ 0)
λ ≡ (Slnρ, V)/(2T2)

The ratio ρ ⁄ ρ0 is needed to ensure self-consistency: the mass PDF must also conserve probability, which means it must be scaled by 1 ⁄ ρ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 δ 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[ρ] ≡  < ρ > V, etc. The variance is computed as V[x] = E[x2] − E[x]2, which means we'll need to compute the moments of various combinations of variables.

Also, Slnρ, V ≡ V[lnρV]

 < ρ > V ≡ ρ0 ≡ ∞∫ − ∞ρPV(lnρ ⁄ ρ0)dlnρ
 < lnρ > V ≡ ∞∫ − ∞lnρPV(lnρ ⁄ ρ0)dlnρ
 < ρ > M ≡ ∞∫ − ∞ρ(ρPV(lnρ ⁄ ρ0))dlnρ
 < lnρ > M ≡ ∞∫ − ∞lnρ(ρPV(lnρ ⁄ ρ0))dlnρ
 < ρ3 > V ≡ ∞∫ − ∞ρ3PV(lnρ ⁄ ρ0)dlnρ
 < ln2ρ > V ≡ ∞∫ − ∞(lnρ)2PV(lnρ ⁄ ρ0)dlnρ
 < ln2ρ > M ≡ ∞∫ − ∞(lnρ)2ρPV(lnρ ⁄ ρ0)dlnρ

In the moment equations, only the input to the PDF is scaled by a mean density ρ0; the "weighting factors" for the expectation value are not (i.e., we are measuring E[ρ], not E[ρ ⁄ rho0]).

Delta terms

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

 < ρ > V ≡ ρ0 ≡ ∞∫ − ∞ρe − λδ((λ)/(1 + T) − (lnρ ⁄ ρ0)/(T))(dlnρ)/(T)
 < lnρ > V ≡ ∞∫ − ∞lnρe − λδ((λ)/(1 + T) − (lnρ ⁄ ρ0)/(T))(dlnρ)/(T)
 < ρ > M ≡ ∞∫ − ∞ρ((ρ)/(ρ0)e − λδ((λ)/(1 + T) − (lnρ ⁄ ρ0)/(T)))(dlnρ)/(T)
 < lnρ > M ≡ ∞∫ − ∞lnρ((ρ)/(ρ0)e − λδ((λ)/(1 + T) − (lnρ ⁄ ρ0)/(T)))(dlnρ)/(T)

Substitution: v = (lnρ ⁄ ρ0)/(T), dv = (1)/(T)dlnρ, ρ = ρ0ev*T, lnρ = vT + lnρ0

 < ρ > Vδ ≡ ρ0 ≡ ∞∫ − ∞ρ0evTe − λδ((λ)/(1 + T) − v)dv
 < lnρ > Vδ ≡ ∞∫ − ∞(vT + lnρ0)e − λδ((λ)/(1 + T) − v)dv
 < ρ > Mδ ≡ ∞∫ − ∞ρ0e2vT(e − λδ((λ)/(1 + T) − v))dv
 < lnρ > Mδ ≡ ∞∫ − ∞(vT + lnρ0)evT(e − λδ((λ)/(1 + T) − v))dv

Solutions:

 < ρ > Vδ = ρ0exp⎡⎣(Tλ)/(1 + T) − λ⎤⎦ = ρ0exp⎡⎣ − λ(1)/(1 + T)⎤⎦
 < lnρ > Vδ = e − λ(λT)/(1 + T) + e − λlnρ0
 < ρ > Mδ = ρ0exp⎡⎣(2Tλ)/(1 + T) − λ⎤⎦ = ρ0exp⎡⎣λ(T − 1)/(T + 1)⎤⎦
 < lnρ > Mδ = ⎛⎝(λT)/(1 + T) + lnρ0⎞⎠exp⎡⎣(Tλ)/(1 + T) − λ⎤⎦
 = ⎛⎝(λT)/(1 + T) + lnρ0⎞⎠exp⎡⎣( − λ)/(T + 1)⎤⎦

(for these next 3, I skipped intermediate steps)

 < ρ3 > Vδ = ρ20exp⎡⎣(3Tλ)/(1 + T) − λ⎤⎦ = ρ20exp⎡⎣λ(2T − 1)/(T + 1)⎤⎦
 < ln2ρ > Mδ = ⎛⎝(λT)/(1 + T) + lnρ0⎞⎠2exp⎡⎣( − λ)/(T + 1)⎤⎦
 < ln2ρ > Vδ = ⎛⎝(λT)/(1 + T) + lnρ0⎞⎠2e − λ

Using ρ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:

∞∫0I1(x)e − x2 ⁄ (4L)dx = eL − 1
∞∫0x2I1(x)e − x2 ⁄ (4L)dx = 4L2*eL
∞∫0x4I1(x)e − x2 ⁄ (4L)dx = 16L3(L + 2)*eL

We use L instead of λ in these equations because it is often substituted in later equations.

Expectation Value of the Volume-Weighted Density E[ρ]

E[ρ] ≡ ∫ρPv(lnρ ⁄ ρ0)dlnρ = ∫ρ[I1(2√(λu))e − λ − u√((λ)/(u)) + e − λδ(u)]du

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

x = 2√(λu)

which gives us ρ in terms of x:

ρ = ρ0exp⎡⎣T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎤⎦

and leads to the rearrangement:

E[ρ] = ∫ρ0exp⎡⎣T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎤⎦[I1(x)e − x2 ⁄ (4λ) − λ]dx + ρ0exp⎛⎝ − (λ)/(1 + T)⎞⎠

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.

L → (λ)/(1 + T)
E[ρ] = ρ0⎡⎣exp⎛⎝ − λ + (Tλ)/(1 + T)⎞⎠∫[I1(x)e − x2 ⁄ (4L)]dx + exp⎛⎝ − (λ)/(1 + T)⎞⎠⎤⎦
 = ρ0⎡⎣exp⎛⎝ − λ + (Tλ)/(1 + T)⎞⎠(eL − 1) + exp⎛⎝ − (λ)/(1 + T)⎞⎠⎤⎦
 = ρ0⎡⎣exp⎛⎝ − λ + (Tλ)/(1 + T)⎞⎠(eλ ⁄ 1 + T − 1) + exp⎛⎝ − (λ)/(1 + T)⎞⎠⎤⎦
 = ρ0⎡⎣e − λ ⁄ (1 + T)(eλ ⁄ 1 + T − 1) + exp⎛⎝ − (λ)/(1 + T)⎞⎠⎤⎦
 = ρ0

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[ρ] = Slnρ, V

V[ρ] = E[ρ2] − E[ρ]2 = ρ20⎡⎣exp⎛⎝λ(2T2)/(1 + 3T + 2T2)⎞⎠ − 1⎤⎦

However, the "correction factor" is still important:

Vδ[ρ] = ρ20⎡⎣exp⎛⎝λ(T − 1)/(T + 1)⎞⎠ − exp⎛⎝ − 2(λ)/(1 + T)⎞⎠⎤⎦

Expectation Value of the Mass-Weighted Density EM[ρ]

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

EM[ρ] = ∫ρ0exp⎡⎣2T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎤⎦[I1(x)e − x2 ⁄ (4λ) − λ]dx + ρ0exp⎛⎝ − (λ(T − 1))/(1 + T)⎞⎠

Follow the same math, with L = (λ)/(1 + 2T)

 = ρ0⎡⎣exp⎛⎝ − λ + (2Tλ)/(1 + T)⎞⎠(eL − 1) + exp⎛⎝ − (λ(T − 1))/(1 + T)⎞⎠⎤⎦
 = ρ0⎡⎣exp⎛⎝((T − 1)λ)/(1 + T)⎞⎠(eλ ⁄ (1 + 2T) − 1) + exp⎛⎝ − (λ)/(1 + T)⎞⎠⎤⎦
EM[ρ] = ρ0⎡⎣exp⎛⎝λ(2T2)/(1 + 3T + 2T2)⎞⎠ − exp⎛⎝λ(T − 1)/(T + 1)⎞⎠ + exp⎛⎝λ(T − 1)/(T + 1)⎞⎠⎤⎦

The right 2 terms cancel, yielding the value shown in Equation 7 of Hopkins 2013 scaled by ρ20. 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.

Eδ, M[ρ] = exp⎛⎝λ(T − 1)/(T + 1)⎞⎠
EM[ρ] = ρ0exp⎛⎝λ(2T2)/(1 + 3T + 2T2)⎞⎠

Expectation Value of the Mass-Weighted Density Squared EM[ρ2]

EM[ρ2] = ∫ρ2(ρ)/(ρ0)e − λ[I1(x)e − x2 ⁄ (4λ)]dx + ∫ρ2(ρ)/(ρ0)e − λδ(u)du
Eδ, M[ρ2] = ρ20exp⎡⎣λ(2T − 1)/(T + 1)⎤⎦
Eleft[ρ2] = e − λ∫ρ20exp⎡⎣3T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎤⎦[I1(x)e − x2 ⁄ (4λ)]dx
 = ρ20exp⎡⎣((2T − 1)λ)/(1 + T)⎤⎦∫I1(x)e − (3T + 1)x2 ⁄ (4λ)dx
 = ρ20exp⎡⎣((2T − 1)λ)/(1 + T)⎤⎦⎛⎝exp⎡⎣(λ)/(3T + 1)⎤⎦ − 1⎞⎠
 = ρ20⎛⎝exp⎡⎣(6λT2)/(3T2 + 4T + 1)⎤⎦ − exp⎡⎣((2T − 1)λ)/(1 + T)⎤⎦⎞⎠
EM[ρ2] = ρ20exp⎡⎣(6λT2)/(3T2 + 4T + 1)⎤⎦

Variance of the Mass-Weighted Density VM[ρ] = EM[ρ2] − EM[ρ]2

Since correction factors are given for EM[ρ2] and EM[ρ], they are not included separately here:

VM[ρ] = EM[ρ2] − EM[ρ]2 = ρ20⎛⎝exp⎡⎣(6λT2)/(3T2 + 4T + 1)⎤⎦ − exp⎡⎣λ(4T2)/(1 + 3T + 2T2)⎤⎦⎞⎠

Expectation Value of the Volume-Weighted Log Density E[lnρ]

E[lnρ] = ∫lnρe − λ[I1(x)e − x2 ⁄ (4λ)]dx + ∫lnρe − λδ(u)du
Eδ[lnρ] = e − λ⎡⎣(λT)/(1 + T) + lnρ0⎤⎦
Eleft[lnρ] = ∫⎡⎣lnρ0 + T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎤⎦e − λ[I1(x)e − x2 ⁄ (4λ)]dx
 = e − λ⎛⎝∫⎡⎣lnρ0 + (Tλ)/(1 + T)⎤⎦[I1(x)e − x2 ⁄ (4λ)]dx − ∫(Tx2)/(4λ)[I1(x)e − x2 ⁄ (4λ)]dx⎞⎠
 = e − λ⎛⎝∫⎡⎣lnρ0 + (Tλ)/(1 + T)⎤⎦(eλ − 1) − (4Tλ2eλ)/(4λ)⎞⎠
 = ⎛⎝⎡⎣lnρ0 + (Tλ)/(1 + T)⎤⎦(1 − e − λ) − Tλ⎞⎠
E[lnρ] = lnρ0 + (Tλ)/(1 + T) − Tλ
 = lnρ0 − (T2λ)/(1 + T)

Expectation Value of the Mass-Weighted Log Density EM[lnρ]

EM[lnρ] = ∫lnρ(ρ)/(ρ0)e − λ[I1(x)e − x2 ⁄ (4λ)]dx + ∫lnρ(ρ)/(ρ0)e − λδ(u)du
Eδ[lnρ] = ⎛⎝(λT)/(1 + T) + lnρ0⎞⎠exp⎡⎣( − λ)/(T + 1)⎤⎦
Eleft[lnρ] = e − λ∫⎡⎣⎛⎝lnρ0 + T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎞⎠exp⎛⎝T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎞⎠⎤⎦[I1(x)e − x2 ⁄ (4λ)]dx

Again, separate into integrable terms:

 = exp⎛⎝(Tλ)/(1 + T) − λ⎞⎠⎡⎣⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠[I1(x)e − (1 + T)x2 ⁄ (4λ)] + ⎛⎝ − (Tx2)/(4λ)⎞⎠[I1(x)e − (1 + T)x2 ⁄ (4λ)]⎤⎦
L = (λ)/(1 + T)
Eleft[lnρ] = exp⎛⎝(Tλ)/(1 + T) − λ⎞⎠⎡⎣⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠⎛⎝exp⎡⎣(λ)/(1 + T)⎤⎦ − 1⎞⎠ + ⎛⎝ − (T)/(4λ)⎞⎠⎛⎝(4λ2)/((1 + T)2)exp⎡⎣(λ)/(1 + T)⎤⎦⎞⎠⎤⎦
Eleft[lnρ] = exp⎛⎝( − λ)/(1 + T)⎞⎠⎡⎣⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠⎛⎝exp⎡⎣(λ)/(1 + T)⎤⎦ − 1⎞⎠ + ⎛⎝ − (T)/(4λ)⎞⎠⎛⎝(4λ2)/((1 + T)2)exp⎡⎣(λ)/(1 + T)⎤⎦⎞⎠⎤⎦
Eleft[lnρ] = ⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠⎛⎝1 − exp⎡⎣( − λ)/(1 + T)⎤⎦⎞⎠ − ⎛⎝(Tλ)/((1 + T)2)⎞⎠
EM[lnρ] = ⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠ − ⎛⎝(Tλ)/((1 + T)2)⎞⎠
 = lnρ0 + (T2λ)/((1 + T)2)

Expectation Value of the Mass-Weighted Log Density Squared EM[ln2ρ]

EM[ln2ρ] = ∫(lnρ)2(ρ)/(ρ0)e − λ[I1(x)e − x2 ⁄ (4λ)]dx + ∫(lnρ)2(ρ)/(ρ0)e − λδ(u)du
Eδ[ln2ρ] = ⎛⎝(λT)/(1 + T) + lnρ0⎞⎠2exp⎡⎣( − λ)/(T + 1)⎤⎦
Eleft[ln2ρ] = e − λ∫⎡⎣⎛⎝lnρ0 + T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎞⎠2exp⎛⎝T⎛⎝ − (x2)/(4λ) + (λ)/(1 + T)⎞⎠⎞⎠⎤⎦[I1(x)e − x2 ⁄ (4λ)]dx

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

Q = lnρ0 + (Tλ)/(1 + T)
Eleft[ln2ρ] = e − λ ⁄ (1 + T)∫⎡⎣⎛⎝Q − (Tx2)/(4λ)⎞⎠2exp⎛⎝ − (Tx2)/(4λ)⎞⎠⎤⎦[I1(x)e − x2 ⁄ (4λ)]dx
Eleft[ln2ρ] = e − λ ⁄ (1 + T)∫⎡⎣⎛⎝Q2 − 2Q(Tx2)/(4λ) + (T2x4)/(16λ2)⎞⎠[I1(x)e − (1 + T)x2 ⁄ (4λ)]⎤⎦dx
Eleft[ln2ρ] = e − λ ⁄ (1 + T)⎡⎣Q2(eλ ⁄ (1 + T) − 1) − 2Q(T)/(4λ)⎛⎝(4λ2)/((1 + T)2)eλ ⁄ (1 + T)⎞⎠ + (T2)/(16λ2)⎛⎝16(λ3)/((1 + T)3)((λ)/(1 + T) + 2)eλ ⁄ (1 + T)⎞⎠⎤⎦
Eleft[ln2ρ] = Q2(1 − e − λ ⁄ (1 + T)) − 2Q(Tλ)/((1 + T)2) + (λT2)/((1 + T)3)⎛⎝(λ)/(1 + T) + 2⎞⎠

Add back the δ termG

EM[ln2ρ] = Q2 − 2Q(Tλ)/((1 + T)2) + (λT2)/((1 + T)3)⎛⎝(λ)/(1 + T) + 2⎞⎠

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

EM[ln2ρ] = ⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠2 − 2⎛⎝lnρ0 + (Tλ)/(1 + T)⎞⎠(Tλ)/((1 + T)2) + (λT2)/((1 + T)3)⎛⎝(λ)/(1 + T) + 2⎞⎠
EM[ln2ρ] = ln2ρ0 + 2lnρ0(Tλ)/(1 + T) + (T2λ2)/((1 + T)2) − 2lnρ0(Tλ)/((1 + T)2) − 2(Tλ)/(1 + T)(Tλ)/((1 + T)2) + 2(λT2)/((1 + T)3) + (λT2)/((1 + T)3)(λ)/(1 + T)
EM[ln2ρ] = ln2ρ0 + 2lnρ0⎛⎝(Tλ)/(1 + T) − (Tλ)/((1 + T)2)⎞⎠ + (T2λ2(1 + T))/((1 + T)3) − 2(T2λ2)/((1 + T)3) + 2(λT2)/((1 + T)3) + (λ2T2)/((1 + T)4)
EM[ln2ρ] = ln2ρ0 + 2lnρ0⎛⎝(T2λ)/(1 + T)⎞⎠ + (T2λ2(1 + T))/((1 + T)3) − 2(T2λ2)/((1 + T)3) + 2(λT2)/((1 + T)3) + (λ2T2)/((1 + T)4)
 = ln2ρ0 + 2lnρ0(T2λ)/((1 + T)2) + ⎛⎝(λT2)/((1 + T)2)⎞⎠2 + (2λT2)/((1 + T)3)
 = ⎛⎝lnρ0 + (T2λ)/((1 + T)2)⎞⎠2 + (2λT2)/((1 + T)3)

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

Slnρ, M = EM[lnρ2] − EM[lnρ]2 = ln2ρ0 + 2lnρ0(T2λ)/((1 + T)2) + ⎛⎝(λT2)/((1 + T)2)⎞⎠2 + (2λT2)/((1 + T)3) − ⎛⎝lnρ0 + (T2λ)/((1 + T)2)⎞⎠2
 = (2λT2)/((1 + T)3)

independent of ρ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 ∼ ℳC. If this is taken on face value, without recognizing the relation between T and σ, it means that for fixed σ,  < ρ > M ∝ e − T2.

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 δ terms; the agreement is essentially perfect excepting some numerical noise (the only visible discrepancy is for a perfect lognormal at high σV):

However, Hopkins 2013 also found that simulations generally produced T ∼ 0.1σ2s, M (though a relation of the form T ∼ 0.25ln(1 + 0.25σ4s, M) 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 σs. In these plots, half are labeled with σ and the others are labeled with T, though the two are equivalent. The highest σs observed in any of the simulations shown in Hopkins 2013 was about σs = 4 (marked with a black square below), so the maximum ρM ⁄ ρV ∼ 102.