Gamma vs Log-Normal Particle Size Distribution
Author: Elspeth K.H. Lee (https://github.com/ELeeAstro)
In this notebook we compare the two particle size distributions available in virga:
Log-normal (default) — parameterised by geometric standard deviation
sigGamma (Christie et al. 2022, Appendix B) — shape parameter
Ais derived from the samesigvia the trigamma function, so the two distributions have matched variance in log-space
We run both on the same Hot Jupiter PT profile and compare particle radii, optical depth, and the dn/dr distribution.
[1]:
from bokeh.io import output_notebook
from bokeh.plotting import show
output_notebook()
import numpy as np
import astropy.units as u
import virga.justdoit as jdi
import virga.justplotit as jpi
We use the built-in Hot Jupiter PT profile. Both runs use identical atmospheric parameters — only the particle size distribution (dist) differs.
[2]:
mieff_dir = '/Users/nbatalh1/Documents/data/virga'
metallicity = 1 # atmospheric metallicity relative to Solar
mean_molecular_weight = 2.2 # atmospheric mean molecular weight
# --- Log-normal run (default) ---
a_ln = jdi.Atmosphere(['MnS','Cr','MgSiO3'],
fsed=1, mh=metallicity, mmw=mean_molecular_weight)
a_ln.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
a_ln.ptk(df=jdi.hot_jupiter())
out_ln = jdi.compute(a_ln, as_dict=True, directory=mieff_dir)
# --- Gamma run (Christie et al. 2022, Appendix B) ---
a_gm = jdi.Atmosphere(['MnS','Cr','MgSiO3'],
fsed=1, mh=metallicity, mmw=mean_molecular_weight,
dist='gamma')
a_gm.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
a_gm.ptk(df=jdi.hot_jupiter())
out_gm = jdi.compute(a_gm, as_dict=True, directory=mieff_dir)
[3]:
print(f"lognormal dist: {a_ln.dist}")
print(f"gamma dist: {a_gm.dist}")
print(f"gamma A: {a_gm.gamma_A:.5f}") # should be ~2.54278 for sig=2
lognormal dist: lognormal
gamma dist: gamma
gamma A: 2.54278
PT Profile with Condensation Curves
Same atmosphere for both runs — this is just a sanity check that the PT profile and condensation curves are identical.
[4]:
show(jpi.pt(out_ln, plot_height=450))
Particle Radii: Log-normal vs Gamma
The left panel shows the mean particle radius profile with pressure. The right panel shows the full dn/dr distribution at a selected pressure level.
For sig=2 (the default), Christie et al. (2022) find the effective radii are similar between the two distributions — this is a good first check.
[5]:
fig, dndr = jpi.radii([out_ln, out_gm], at_pressure=1e-3,
compare=True, legend=['lognormal', 'gamma'])
show(fig)
[6]:
print("lognormal rg (top layer):", out_ln['mean_particle_r'][0])
print("gamma rg (top layer):", out_gm['mean_particle_r'][0])
print("lognormal reff (top layer):", out_ln['droplet_eff_r'][0])
print("gamma reff (top layer):", out_gm['droplet_eff_r'][0])
print("gamma_A:", out_gm['scalar_inputs']['gamma_A'])
lognormal rg (top layer): [0.01972209 0.01103366 0.02471404]
gamma rg (top layer): [0.04862254 0.02720222 0.06092963]
lognormal reff (top layer): [0.06555384 0.03667455 0.08214648]
gamma reff (top layer): [0.0868661 0.04859784 0.1088532 ]
gamma_A: 2.5427835527445835
Optical Depth, Single Scattering, Asymmetry
Compare the full wavelength-resolved optics between the two distributions.
[7]:
show(jpi.all_optics(out_ln))
show(jpi.all_optics(out_gm))
1D Optics (wavelength-averaged)
Averaged over 1–2 microns for a quick visual comparison.
[8]:
show(jpi.all_optics_1d(out_ln, wave_range=[1,2]))
show(jpi.all_optics_1d(out_gm, wave_range=[1,2]))
Cumulative Optical Depth by Gas
[9]:
show(jpi.opd_by_gas(out_ln))
show(jpi.opd_by_gas(out_gm))
Mean Mass Mixing Ratio of the Condensates
[10]:
show(jpi.condensate_mmr(out_ln))
show(jpi.condensate_mmr(out_gm))
Extended Comparison of Different Distribtuions over Sigma and Fsed
The sections above used the default sig=2. Christie et al. (2022) show that the two distributions agree well at sig=2 but diverge at larger values — the gamma distribution peak stays close to r_w while the lognormal peak shifts significantly. The sections below explore this systematically.
Sweep over sigma
We run both distributions for sig = [1.5, 2.0, 2.5, 3.0] with all other parameters fixed. This is the key test from Christie et al. (2022): the distributions should agree at sig=2 and diverge as sig increases, with the gamma distribution producing fewer very large particles than the lognormal.
Particle radii and dn/dr at 1e-3 bar for each sigma — lognormal vs gamma overlaid:
[11]:
sig_values = [1.1, 1.5, 2.0, 2.5, 3.0]
outputs_ln_sig = []
outputs_gm_sig = []
for sig in sig_values:
# lognormal
a = jdi.Atmosphere(['MnS','Cr','MgSiO3'],
fsed=1, mh=metallicity, mmw=mean_molecular_weight, sig=sig)
a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
a.ptk(df=jdi.hot_jupiter())
outputs_ln_sig.append(jdi.compute(a, as_dict=True, directory=mieff_dir))
# gamma
a = jdi.Atmosphere(['MnS','Cr','MgSiO3'],
fsed=1, mh=metallicity, mmw=mean_molecular_weight, sig=sig,
dist='gamma')
a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
a.ptk(df=jdi.hot_jupiter())
outputs_gm_sig.append(jdi.compute(a, as_dict=True, directory=mieff_dir))
print("Done. sig values:", sig_values)
Done. sig values: [1.1, 1.5, 2.0, 2.5, 3.0]
[12]:
from bokeh.layouts import column
all_figs = []
for i, sig in enumerate(sig_values):
fig, _ = jpi.radii([outputs_ln_sig[i], outputs_gm_sig[i]], at_pressure=1e-3,
compare=True, legend=['lognormal', 'gamma'])
all_figs.append(fig)
# Link x/y ranges across rows so all panels share the same scale
p1_ref = all_figs[0].children[0] # mean radius vs pressure
p2_ref = all_figs[0].children[1] # dn/dr distribution
for fig in all_figs[1:]:
fig.children[0].x_range = p1_ref.x_range
fig.children[1].x_range = p2_ref.x_range
fig.children[1].y_range = p2_ref.y_range
show(column(*all_figs))
Sweep over fsed
We fix sig=2 and sweep fsed = [0.1, 0.5, 1.0, 3.0]. Lower fsed means more lofted, larger particles — this tests whether the distribution choice matters more in the thick cloud regime.
[13]:
fsed_values = [0.1, 0.5, 1.0, 3.0]
outputs_ln_fsed = []
outputs_gm_fsed = []
for fsed in fsed_values:
# lognormal
a = jdi.Atmosphere(['MnS','Cr','MgSiO3'],
fsed=fsed, mh=metallicity, mmw=mean_molecular_weight)
a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
a.ptk(df=jdi.hot_jupiter())
outputs_ln_fsed.append(jdi.compute(a, as_dict=True, directory=mieff_dir))
# gamma
a = jdi.Atmosphere(['MnS','Cr','MgSiO3'],
fsed=fsed, mh=metallicity, mmw=mean_molecular_weight,
dist='gamma')
a.gravity(gravity=7.460, gravity_unit=u.Unit('m/(s**2)'))
a.ptk(df=jdi.hot_jupiter())
outputs_gm_fsed.append(jdi.compute(a, as_dict=True, directory=mieff_dir))
print("Done. fsed values:", fsed_values)
Done. fsed values: [0.1, 0.5, 1.0, 3.0]
[14]:
from bokeh.layouts import column
all_figs_r = []
all_figs_m = []
for i, fsed in enumerate(fsed_values):
fig_r, _ = jpi.radii([outputs_ln_fsed[i], outputs_gm_fsed[i]], at_pressure=1e-3,
compare=True, legend=['lognormal', 'gamma'])
fig_m = jpi.condensate_mmr([outputs_ln_fsed[i], outputs_gm_fsed[i]])
all_figs_r.append(fig_r)
all_figs_m.append(fig_m)
# Link ranges for radii panels so all rows share the same scale
p1_ref = all_figs_r[0].children[0] # mean radius vs pressure
p2_ref = all_figs_r[0].children[1] # dn/dr distribution
for fig in all_figs_r[1:]:
fig.children[0].x_range = p1_ref.x_range
fig.children[1].x_range = p2_ref.x_range
fig.children[1].y_range = p2_ref.y_range
show(column(*[f for pair in zip(all_figs_r, all_figs_m) for f in pair]))
Optics comparison
For sig=2 (where distributions agree) and sig=3 (where they diverge), compare 1D wavelength-averaged optics. Christie et al. (2022) find negligible differences in transmission spectra at sig=2 but larger differences at wider distributions.
[15]:
from bokeh.layouts import column
all_figs = []
for label, out in [('lognormal sig=2', out_ln), ('gamma sig=2', out_gm),
('lognormal sig=3', outputs_ln_sig[3]), ('gamma sig=3', outputs_gm_sig[3])]:
print(f"── {label} ──")
all_figs.append(jpi.all_optics_1d(out, wave_range=[0.5, 5]))
show(column(*all_figs))
── lognormal sig=2 ──
── gamma sig=2 ──
── lognormal sig=3 ──
── gamma sig=3 ──