#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ neptune-aurora.py Science Journaling Club, own simplified calculation. Companion to: articles/neptune-aurora.html WHAT THIS SCRIPT IS ------------------- Two small models, both ours, both deliberately crude, written to sit beside Melin et al. (2025), "Discovery of H3+ and infrared aurorae at Neptune with JWST", Nature Astronomy 9, 666-671, doi:10.1038/s41550-025-02507-9. The paper's authors measured things. We did not. Everything below is arithmetic performed on published constants by a school club with a laptop, and it is here so that a reader can see the shape of the two claims for themselves. MODEL A -- where a tilted, offset dipole puts the auroral oval. MODEL B -- how fast H3+ infrared emission dies when the gas cools. MODEL A: THE OFFSET TILTED DIPOLE (OTD) --------------------------------------- Neptune's internal magnetic field is genuinely messy: the Voyager 2 fit (the "O8" model of Connerney, Acuna & Ness 1991, doi:10.1029/91JA01165) needs strong quadrupole and octupole terms. We throw all of that away and keep the two headline numbers from Ness et al. (1989, doi:10.1126/science.246.4936.1473): a dipole tilted about 47 degrees from the rotation axis and displaced about 0.55 Neptune radii from the planet's centre. Geometry, in planet-centred Cartesian coordinates with lengths in R_N (planet radius) and +z along the rotation axis: dipole axis a_hat = (sin t, 0, cos t) t = tilt angle dipole centre d = delta * u_hat delta = offset magnitude field direction B(r) proportional to 3 (a_hat . s_hat) s_hat - a_hat where s = r - d (position relative to the DIPOLE centre) The magnitude of the moment never enters, because we only ever follow the direction of B. Field lines are traced by RK4 integration of dr/ds = B_hat(r), started on the dipole's magnetic equator at a chosen L value (the equatorial crossing distance in R_N, measured from the dipole centre) and stopped when the path first reaches r = 1, which we are calling the cloud tops. The last step is refined by bisection so the footpoint sits on the surface to 1e-9 R_N. Seeding 360 degrees of azimuth around the dipole axis at fixed L traces out a closed curve of footpoints. That curve is the auroral oval in this model. ASSUMPTIONS, STATED PLAINLY A1. Dipole only. No quadrupole, no octupole. Neptune has both, strongly. A2. Internal field only. No magnetopause, no ring current, no solar-wind compression, so the oval is symmetric about the dipole axis in a way the real one cannot be. A3. The planet is a sphere of radius 1 R_N. Neptune's real flattening is about 1.7%, which moves footpoint latitudes by well under a degree at these colatitudes but is not zero. A4. "Cloud tops" means r = 1 R_N. The H3+ that JWST sees actually radiates from roughly 1000 km higher, about 0.04 R_N, which pushes footpoints very slightly poleward in this model. A5. The offset DIRECTION is not published in a form we could verify, only its magnitude. So the headline run uses a CENTRED tilted dipole, which spends only the one number everybody quotes, and the offset is then swept through every direction from axial to anti-axial as a separate study. This matters enormously: an offset pointed along the dipole axis leaves the whole system axisymmetric and cannot change the oval's shape at all, while an offset lying across the axis drags the footpoints tens of degrees. The sweep is section A3 of the output, and the observed latitudes turn out to constrain the direction. A6. The oval is taken to be the footprint of a single L shell. Real auroral ovals have width because precipitation happens over a range of L. A7. Planetocentric latitude, not planetographic. On a sphere they are the same thing; on the real oblate Neptune they differ by up to ~1 degree. LIMITS This model cannot tell you where the aurora is in longitude in any useful way, and neither can the paper: Neptune's rotation period is uncertain enough that the rotational phase at the time of the JWST observation is not known. So we compare LATITUDES only, and we say so in the article. It also cannot tell you how bright the aurora is, which requires a precipitation model we have no business writing. MODEL B: H3+ BRIGHTNESS AGAINST TEMPERATURE -------------------------------------------- H3+ is a triangular molecular ion, three protons and two electrons, and it is the thing that makes giant-planet upper atmospheres visible in the infrared. It radiates in the nu2 bending fundamental near 3.3 to 4.0 micrometres. The population of any emitting level follows a Boltzmann factor, so the emitted power per ion is eps(T) = (1 / Q(T)) * SUM over upper levels u of g_ns(K_u) (2J_u + 1) exp(-E_u / kT) * A_eff * h c nu_eff with Q(T) the total internal partition function. We build both sums explicitly from rigid oblate symmetric-top energy levels, E(J, K) = E_vib + B J(J+1) + (C - B) K^2 using B0 = 43.5685 cm-1 and C0 = 20.6170 cm-1 for the ground vibrational state and a nu2 band origin of 2521.31 cm-1. Nuclear spin statistics for three identical protons give a weight of 4 when K is divisible by 3 and 2 otherwise, and the (K = 0, J even) levels of the ground state do not exist at all. ASSUMPTIONS B1. Local thermodynamic equilibrium. At the pressures where Neptune's H3+ lives this is roughly true and is what everyone assumes. B2. Rigid rotor, no centrifugal distortion, no l-type doubling, no splitting of the nu2 = 1 vibrational angular momentum states. The real level structure is messier. B3. A single effective Einstein A coefficient and a single effective photon energy for the whole band. Both cancel exactly in the RATIO, which is the only quantity we quote, so this assumption costs us nothing in the headline number and everything in absolute radiance. We never quote an absolute radiance. B4. Constant H3+ column density. Cooling would in reality also change the ion chemistry and the column. The paper measures the column directly, so we do not need to guess it. LIMITS The ratio is the whole result. Do not read an absolute brightness out of this code, because it is not in there. MODEL C: A SENSITIVITY LADDER (the crudest thing here) ------------------------------------------------------ To ask why 36 years of searching failed we need a rough idea of how faint each generation of instrument could go. We use a photon-counting figure of merit, FOM = sqrt(A * t / B), with A the collecting area, t the on-source integration and B a penalty factor for the thermal sky and telescope background that ground-based 3 to 4 micrometre observing suffers and space-based observing does not. The detection floor is then k / FOM, and k is fixed by one calibration: JWST detected the emission, so we set JWST's floor a factor of two below the brightness the paper actually measured. Everything else follows. C1. Background-limited, Poisson noise, nothing else. C2. The background penalty B = 100 for the ground is an order-of-magnitude judgement, not a measurement. C3. Integration times for the historical attempts are our estimates. This ladder is an illustration of a scaling, not a replication of anybody's published upper limit. It gets one of its four cases wrong, and we say so. REPRODUCING python neptune-aurora.py > neptune-aurora-output.txt Requires numpy. All random draws are seeded (SEED below). """ import math import numpy as np SEED = 20250326 # the paper's publication date, 26 March 2025 np.random.seed(SEED) # ---------------------------------------------------------------------------- # Published constants we take from the literature. Nothing here is ours. # ---------------------------------------------------------------------------- TILT_DEG = 46.9 # dipole tilt from rotation axis, Ness et al. 1989 [2] OFFSET_RN = 0.55 # dipole displacement from planet centre, in R_N [2] R_NEPTUNE_KM = 24764.0 # equatorial 1-bar radius T_VOYAGER = 750.0 # exospheric temperature, Voyager 2 UVS 1989, +/-150 K T_VOYAGER_ERR = 150.0 T_JWST = 358.0 # Melin et al. 2025, globally averaged, +/- 8 K T_JWST_ERR = 8.0 OBS_LAT_N = -30.0 # equatorward edge of the JWST auroral patch OBS_LAT_S = -60.0 # poleward edge PAPER_RATIO = 0.008 # "at 358 K the intensity of H3+ is 0.8% of the # intensity at 750 K for the same density" [1] # H3+ spectroscopy, Lindsay & McCall 2001 [13] and Neale et al. 1996 [12] B0_GROUND = 43.5685 # cm-1 C0_GROUND = 20.6170 # cm-1 NU2_ORIGIN = 2521.31 # cm-1, the bending fundamental band origin B1_NU2 = 43.4 # cm-1, nu2 = 1 rotational constants, approximate C1_NU2 = 20.0 # cm-1 NU2_2_ORIGIN = 4998.0 # cm-1, the 2nu2 manifold, lumped, for Q(T) only K_CM = 0.69503476 # Boltzmann constant in cm-1 per kelvin JMAX = 22 # rotational levels summed to J = 22 # ============================================================================ # MODEL A -- offset tilted dipole and field-line tracing # ============================================================================ def dipole_frame(tilt_deg, offset, beta_deg): """Return (a_hat, d): dipole axis unit vector and dipole centre position. beta_deg is the angle between the offset vector and the dipole axis, measured inside the tilt plane. beta = 0 points the offset straight at the southern magnetic pole; beta = 90 puts it square across the axis; beta = 180 points it at the northern magnetic pole. offset = 0 makes beta irrelevant and gives a plain centred dipole. """ t = math.radians(tilt_deg) b = math.radians(beta_deg) a_hat = np.array([math.sin(t), 0.0, math.cos(t)]) e1 = np.array([math.cos(t), 0.0, -math.sin(t)]) # perp to a_hat, in x-z d = offset * (math.cos(b) * (-a_hat) + math.sin(b) * e1) return a_hat, d def bhat(r, a_hat, d): """Unit vector along B for a dipole of axis a_hat sitting at position d.""" s = r - d sn = math.sqrt(s[0] * s[0] + s[1] * s[1] + s[2] * s[2]) if sn < 1e-12: return np.zeros(3) sh = s / sn b = 3.0 * float(np.dot(a_hat, sh)) * sh - a_hat bn = math.sqrt(b[0] * b[0] + b[1] * b[1] + b[2] * b[2]) if bn < 1e-12: return np.zeros(3) return b / bn def trace_to_surface(start, a_hat, d, max_steps=40000): """RK4-trace a field line from `start` along +B until it reaches r = 1. Returns (footpoint_xyz, n_steps) or (None, n_steps) if it never arrives. Step length scales with distance from the dipole so the trace is fine near the planet and coarse out in the magnetosphere. """ r = np.array(start, dtype=float) for n in range(max_steps): s = r - d sn = math.sqrt(float(np.dot(s, s))) h = min(0.25, max(0.0015, 0.04 * sn)) k1 = bhat(r, a_hat, d) if not k1.any(): return None, n k2 = bhat(r + 0.5 * h * k1, a_hat, d) k3 = bhat(r + 0.5 * h * k2, a_hat, d) k4 = bhat(r + h * k3, a_hat, d) rn = r + (h / 6.0) * (k1 + 2.0 * k2 + 2.0 * k3 + k4) if float(np.dot(rn, rn)) <= 1.0: # bisect the final step onto the sphere r = 1 lo, hi = 0.0, h for _ in range(60): mid = 0.5 * (lo + hi) q1 = bhat(r, a_hat, d) q2 = bhat(r + 0.5 * mid * q1, a_hat, d) q3 = bhat(r + 0.5 * mid * q2, a_hat, d) q4 = bhat(r + mid * q3, a_hat, d) rm = r + (mid / 6.0) * (q1 + 2.0 * q2 + 2.0 * q3 + q4) if float(np.dot(rm, rm)) <= 1.0: hi = mid else: lo = mid q1 = bhat(r, a_hat, d) q2 = bhat(r + 0.5 * hi * q1, a_hat, d) q3 = bhat(r + 0.5 * hi * q2, a_hat, d) q4 = bhat(r + hi * q3, a_hat, d) rf = r + (hi / 6.0) * (q1 + 2.0 * q2 + 2.0 * q3 + q4) return rf / math.sqrt(float(np.dot(rf, rf))), n if float(np.dot(rn, rn)) > 4000.0: return None, n r = rn return None, max_steps def oval(tilt_deg, L, offset=0.0, beta_deg=90.0, n_az=72): """Footpoints of the southern oval: the L-shell traced from the magnetic equator at n_az azimuths. Returns (lats_deg, lons_deg) as arrays.""" a_hat, d = dipole_frame(tilt_deg, offset, beta_deg) t = math.radians(tilt_deg) e1 = np.array([math.cos(t), 0.0, -math.sin(t)]) e2 = np.array([0.0, 1.0, 0.0]) lats, lons = [], [] for i in range(n_az): phi = 2.0 * math.pi * i / n_az seed = d + L * (math.cos(phi) * e1 + math.sin(phi) * e2) fp, _ = trace_to_surface(seed, a_hat, d) if fp is None: continue lats.append(math.degrees(math.asin(max(-1.0, min(1.0, fp[2]))))) lons.append(math.degrees(math.atan2(fp[1], fp[0])) % 360.0) return np.array(lats), np.array(lons) # ============================================================================ # MODEL B -- H3+ emission against temperature # ============================================================================ def levels(vib_origin, B, C, forbid_k0_parity): """Enumerate (energy_cm1, statistical_weight) for one vibrational state. forbid_k0_parity: 0 forbids (K=0, J even), 1 forbids (K=0, J odd). Nuclear spin weight is 4 when K is a multiple of 3, 2 otherwise. K > 0 carries an extra factor of 2 for the degenerate +K/-K pair. """ out = [] for J in range(0, JMAX + 1): for K in range(0, J + 1): if K == 0 and (J % 2) == forbid_k0_parity: continue g_ns = 4 if (K % 3 == 0) else 2 g_k = 1 if K == 0 else 2 E = vib_origin + B * J * (J + 1) + (C - B) * K * K out.append((E, g_ns * g_k * (2 * J + 1))) return out GROUND = levels(0.0, B0_GROUND, C0_GROUND, forbid_k0_parity=0) NU2_1 = levels(NU2_ORIGIN, B1_NU2, C1_NU2, forbid_k0_parity=1) NU2_2 = levels(NU2_2_ORIGIN, B1_NU2, C1_NU2, forbid_k0_parity=0) def partition(T): """Total internal partition function, ground + nu2 + 2nu2.""" tot = 0.0 for state in (GROUND, NU2_1, NU2_2): for E, g in state: tot += g * math.exp(-E / (K_CM * T)) return tot def band_emission(T): """Relative power per ion radiated in the nu2 fundamental, LTE. Absolute scale is arbitrary and identical at every T, so ratios are exact within the model.""" s = 0.0 for E, g in NU2_1: s += g * math.exp(-E / (K_CM * T)) return s / partition(T) # ============================================================================ # MODEL C -- the sensitivity ladder # ============================================================================ INSTRUMENTS = [ # name, aperture diameter (m), on-source integration (s), background penalty ("Voyager 2 IRIS (1989 flyby)", 0.50, 600.0, 1.0), ("IRTF 3.0 m / iSHELL (2017)", 3.00, 25200.0, 100.0), ("Keck II 10 m / NIRSPEC (2009)", 10.00, 7200.0, 100.0), ("JWST 6.5 m / NIRSpec IFU (2023)", 6.50, 3420.0, 1.0), ] def fom(diam, t, bg): area = math.pi * (diam / 2.0) ** 2 return math.sqrt(area * t / bg) # ============================================================================ # RUN # ============================================================================ def line(ch='-', n=76): print(ch * n) def main(): print("=" * 76) print("NEPTUNE'S AURORA: THE CLUB'S OWN CALCULATION") print("Science Journaling Club | seed = %d | numpy %s" % (SEED, np.__version__)) print("Every number below is produced by this script. None of it is a") print("measurement, and none of it is the paper's own analysis.") print("=" * 76) # ---------------------------------------------------------------- MODEL A print() print("MODEL A -- WHERE A TILTED, OFFSET DIPOLE PUTS THE OVAL") line() print("Field lines traced by RK4 from the magnetic equator at fixed L down") print("to r = 1 R_N. Southern footpoints only. Latitudes are planetocentric") print("and negative in the southern hemisphere.") print() print("A0. Sanity check, centred dipole, zero tilt.") print(" Textbook: cos^2(lat) = 1/L, so lat = -arccos(sqrt(1/L)).") print(" %-6s %-14s %-14s %s" % ("L", "traced", "analytic", "diff")) for L in (5.0, 10.0, 20.0, 30.0): la, _ = oval(0.0, L, offset=0.0, n_az=8) ana = -math.degrees(math.acos(math.sqrt(1.0 / L))) print(" %-6.1f %-14.4f %-14.4f %+.4f deg" % (L, la.mean(), ana, la.mean() - ana)) print(" The tracer agrees with the closed-form dipole to better than") print(" 0.001 degrees, so the integrator is not the source of any error.") print() print("A1. The migration. Southern oval at L = 20 R_N, centred dipole.") print(" The oval is a ring of constant magnetic colatitude, so the tilt") print(" simply carries it down the planet. Nothing subtle happens here.") print() print(" %-7s %-12s %-12s %-13s %-9s" % ("tilt", "poleward", "centre", "equatorward", "width")) print(" %-7s %-12s %-12s %-13s %-9s" % ("(deg)", "edge (deg)", "(deg)", "edge (deg)", "(deg)")) tilt_grid = [0, 5, 10, 15, 20, 25, 30, 35, 40, 43, 46.9] mig = [] for t in tilt_grid: la, _ = oval(t, 20.0, offset=0.0, n_az=72) mig.append((t, la.min(), la.mean(), la.max(), la.max() - la.min())) print(" %-7.1f %-12.2f %-12.2f %-13.2f %-9.2f" % (t, la.min(), la.mean(), la.max(), la.max() - la.min())) print() print(" At zero tilt the ring sits at a single latitude, %.2f deg," % mig[0][3]) print(" %.2f deg from the pole, because a circle centred on the pole is" % (90.0 + mig[0][3])) print(" a line of latitude. Tilt it and the ring keeps its shape but no") print(" longer lies flat: its equatorward edge runs from %.2f deg to" % mig[0][3]) print(" %.2f deg, a walk of %.2f degrees, which is the tilt angle to the" % (mig[-1][3], abs(mig[-1][3] - mig[0][3]))) print(" second decimal place. In a centred dipole that is all a tilt can") print(" do, and it is enough.") print() print("A2. Add the offset and the ring stops lying flat against the lines") print(" of latitude. Same sweep, L = 20, offset 0.55 R_N held across") print(" the dipole axis.") print() print(" %-7s %-12s %-12s %-13s %-9s" % ("tilt", "poleward", "centre", "equatorward", "width")) mig_off = [] for t in tilt_grid: la, _ = oval(t, 20.0, offset=OFFSET_RN, beta_deg=90.0, n_az=72) mig_off.append((t, la.min(), la.mean(), la.max(), la.max() - la.min())) print(" %-7.1f %-12.2f %-12.2f %-13.2f %-9.2f" % (t, la.min(), la.mean(), la.max(), la.max() - la.min())) print() print(" The latitude width swings from %.1f deg down to %.1f deg and" % (max(m[4] for m in mig_off), min(m[4] for m in mig_off))) print(" back again as the ring is carried across the pole and down the") print(" far side. The ring itself never changes shape. What changes is") print(" how it lies against the lines of latitude, and that is the") print(" distortion a reader sees in Figure 3.") print() print("A3. The offset direction, which nobody has published in a form we") print(" could check. Tilt fixed at %.1f deg, L = 20, offset 0.55 R_N," % TILT_DEG) print(" swept from pointing at the southern magnetic pole (beta = 0)") print(" round to pointing at the northern one (beta = 180).") print() print(" %-7s %-10s %-10s %-12s %-8s %-8s %s" % ("beta", "poleward", "centre", "equatorward", "width", "miss", "overlap")) beta_scan = [] obs_mid = 0.5 * (OBS_LAT_N + OBS_LAT_S) for b in (0, 15, 30, 45, 60, 75, 90, 120, 150, 180): la, _ = oval(TILT_DEG, 20.0, offset=OFFSET_RN, beta_deg=float(b), n_az=72) lo_, hi_ = la.min(), la.max() ov = max(0.0, min(hi_, OBS_LAT_N) - max(lo_, OBS_LAT_S)) miss = abs(la.mean() - obs_mid) beta_scan.append((b, lo_, la.mean(), hi_, hi_ - lo_, ov, miss)) print(" %-7.0f %-10.2f %-10.2f %-12.2f %-8.2f %-8.2f %.0f%%" % (b, lo_, la.mean(), hi_, hi_ - lo_, miss, 100.0 * ov / 30.0)) print() print(" 'miss' is the distance from our band centre to the centre of the") print(" observed band. Prefer it to the overlap column, which rewards a") print(" model for being vague: a wide enough oval overlaps everything.") print() best = min(beta_scan, key=lambda r: r[6]) print(" Smallest miss: beta = %.0f deg, %.2f degrees out." % (best[0], best[6])) print(" A displacement lying across the axis, beta near 90, throws the") print(" oval %.0f degrees away from the emission and is ruled out on" % beta_scan[6][6]) print(" these latitudes alone. Displacements close to the axis, in") print(" either direction, keep the oval where JWST found it. So the") print(" observation does constrain the offset geometry, and it says the") print(" displacement runs near the dipole axis. We did not put that in.") print() print("A4. Against the observation, centred dipole, tilt %.1f deg." % TILT_DEG) print() la47, _ = oval(TILT_DEG, 20.0, offset=0.0, n_az=144) pred_lo, pred_hi = la47.min(), la47.max() print(" JWST enhanced H3+ column, southern hemisphere [1]:") print(" %.1f to %.1f degrees latitude" % (OBS_LAT_S, OBS_LAT_N)) print(" Our L = 20 oval:") print(" %.1f to %.1f degrees latitude" % (pred_lo, pred_hi)) ov_lo = max(pred_lo, OBS_LAT_S) ov_hi = min(pred_hi, OBS_LAT_N) overlap = max(0.0, ov_hi - ov_lo) print(" Overlap: %.1f degrees of the observed %.0f degree band," % (overlap, OBS_LAT_N - OBS_LAT_S)) print(" which is %.0f%% of it." % (100.0 * overlap / 30.0)) print(" Our band centre %.1f deg; observed band centre %.1f deg;" % (la47.mean(), 0.5 * (OBS_LAT_N + OBS_LAT_S))) print(" miss distance %.1f degrees." % abs(la47.mean() - 0.5 * (OBS_LAT_N + OBS_LAT_S))) print() print(" For contrast, the same model with the tilt set to zero, which is") print(" where every textbook picture of a giant-planet aurora puts it:") la0, _ = oval(0.0, 20.0, offset=0.0, n_az=72) print(" %.1f to %.1f degrees. That is %.0f degrees of latitude" % (la0.min(), la0.max(), abs(la0.mean() - la47.mean()))) print(" poleward of where the emission actually is.") print() print(" And the L shells. The paper maps the southern emission to field") print(" lines running from under 3 out to 20 R_N [1]. Here is what our") print(" model does with that span, tilt %.1f deg, centred:" % TILT_DEG) print() print(" %-8s %-16s %-12s %-12s" % ("L (R_N)", "mag. latitude", "poleward", "equatorward")) Lband = [] for L in (3, 5, 10, 15, 20, 30): la, _ = oval(TILT_DEG, float(L), offset=0.0, n_az=48) maglat = math.degrees(math.acos(math.sqrt(1.0 / L))) Lband.append((L, maglat, la.min(), la.max())) print(" %-8d %-16.2f %-12.2f %-12.2f" % (L, maglat, la.min(), la.max())) print() print("A5. Monte Carlo over the things we do not know exactly.") print(" tilt ~ N(46.9, 1.0) deg offset ~ U(0, 0.55) R_N") print(" beta ~ U(0, 180) deg L ~ U(15, 25) R_N") print(" azimuth around the oval ~ U(0, 360) deg") print(" 2000 seeded draws, one traced field line each. The offset") print(" direction is drawn flat because we genuinely do not know it.") print() n_mc = 2000 tilts = np.random.normal(TILT_DEG, 1.0, n_mc) offs = np.random.uniform(0.0, OFFSET_RN, n_mc) betas = np.random.uniform(0.0, 180.0, n_mc) Ls = np.random.uniform(15.0, 25.0, n_mc) azs = np.random.uniform(0.0, 2.0 * math.pi, n_mc) mc = [] for i in range(n_mc): a_hat, d = dipole_frame(float(tilts[i]), float(offs[i]), float(betas[i])) tt = math.radians(float(tilts[i])) e1 = np.array([math.cos(tt), 0.0, -math.sin(tt)]) e2 = np.array([0.0, 1.0, 0.0]) seed = d + Ls[i] * (math.cos(azs[i]) * e1 + math.sin(azs[i]) * e2) fp, _ = trace_to_surface(seed, a_hat, d) if fp is None: continue mc.append(math.degrees(math.asin(max(-1.0, min(1.0, fp[2]))))) mc = np.array(mc) q = np.percentile(mc, [2.5, 16, 50, 84, 97.5]) print(" draws that reached the surface: %d of %d" % (len(mc), n_mc)) print(" median footpoint latitude : %.2f deg" % q[2]) print(" 68%% interval : %.2f to %.2f deg" % (q[1], q[3])) print(" 95%% interval : %.2f to %.2f deg" % (q[0], q[4])) inside = float(np.mean((mc >= OBS_LAT_S) & (mc <= OBS_LAT_N))) print(" fraction landing inside the observed %.0f to %.0f deg band: %.1f%%" % (OBS_LAT_S, OBS_LAT_N, 100.0 * inside)) print() print(" Now the same run with the offset held at zero, which is the") print(" model the article headlines:") mc0 = [] for i in range(n_mc): a_hat, d = dipole_frame(float(tilts[i]), 0.0, 0.0) tt = math.radians(float(tilts[i])) e1 = np.array([math.cos(tt), 0.0, -math.sin(tt)]) e2 = np.array([0.0, 1.0, 0.0]) seed = d + Ls[i] * (math.cos(azs[i]) * e1 + math.sin(azs[i]) * e2) fp, _ = trace_to_surface(seed, a_hat, d) if fp is None: continue mc0.append(math.degrees(math.asin(max(-1.0, min(1.0, fp[2]))))) mc0 = np.array(mc0) q0 = np.percentile(mc0, [2.5, 16, 50, 84, 97.5]) inside0 = float(np.mean((mc0 >= OBS_LAT_S) & (mc0 <= OBS_LAT_N))) print(" median %.2f deg, 68%% %.2f to %.2f, 95%% %.2f to %.2f" % (q0[2], q0[1], q0[3], q0[0], q0[4])) print(" fraction inside the observed band: %.1f%%" % (100.0 * inside0)) print() print(" Read the two together. Knowing the tilt alone puts you in the") print(" right latitudes %.0f times out of a hundred. Letting the offset" % (100.0 * inside0)) print(" point anywhere drops that to %.0f. The honest statement is that" % (100.0 * inside)) print(" the 47 degree tilt is what predicts a mid-latitude aurora, and") print(" the offset is a source of error we cannot close from here.") print() # ---------------------------------------------------------------- MODEL B print() print("MODEL B -- WHY THE GLOW WENT OUT") line() print("LTE nu2 band emission per H3+ ion, rigid symmetric top, levels to") print("J = %d, nuclear spin weights 4 for K divisible by 3 and 2 otherwise." % JMAX) print() print(" %-8s %-14s %-16s %-14s" % ("T (K)", "Q(T)", "eps(T) rel.", "vs 750 K")) e750 = band_emission(T_VOYAGER) for T in (200, 250, 300, 350, T_JWST, 400, 450, 500, 550, 600, 650, 700, 750, 800, 900, 1000): e = band_emission(float(T)) print(" %-8.0f %-14.3f %-16.6e %-14s" % (T, partition(float(T)), e / e750, ("%.4f%%" % (100.0 * e / e750)) if e / e750 < 0.1 else ("%.3f" % (e / e750)))) print() e358 = band_emission(T_JWST) e350 = band_emission(350.0) r358 = e358 / e750 r350 = e350 / e750 print(" Cooling from %.0f K to %.0f K:" % (T_VOYAGER, T_JWST)) print(" emission per ion falls to %.4f%% of its old value," % (100.0 * r358)) print(" a factor of %.1f fainter." % (1.0 / r358)) print(" Cooling from %.0f K to 350 K:" % T_VOYAGER) print(" %.4f%%, a factor of %.1f." % (100.0 * r350, 1.0 / r350)) print() print(" The paper states 0.8% at 358 K for the same density [1].") print(" We get %.4f%%. Our answer is a factor of %.2f from theirs," % (100.0 * r358, (PAPER_RATIO / r358) if r358 < PAPER_RATIO else (r358 / PAPER_RATIO))) print(" which for a rigid-rotor model with a made-up band-average") print(" Einstein coefficient is about as well as we deserve to do.") print() print("B1. How much of that is the Boltzmann factor and how much is Q?") ratio_boltz = (sum(g * math.exp(-E / (K_CM * T_JWST)) for E, g in NU2_1) / sum(g * math.exp(-E / (K_CM * T_VOYAGER)) for E, g in NU2_1)) ratio_q = partition(T_VOYAGER) / partition(T_JWST) print(" population of the emitting levels, 358 K over 750 K : %.4e" % ratio_boltz) print(" which on its own is a fall by a factor of : %.0f" % (1.0 / ratio_boltz)) print(" partition function Q(750)/Q(358) : %.4f" % ratio_q) print(" net emission ratio, the product : %.4e" % (ratio_boltz * ratio_q)) print() print(" The two effects pull against each other. Cooling empties the") print(" emitting levels by a factor of %.0f. Cooling also crowds the ions" % (1.0 / ratio_boltz)) print(" into fewer states, which hands back a factor of %.2f. The fight" % ratio_q) print(" is not close. Net, a factor of %.0f." % (1.0 / (ratio_boltz * ratio_q))) print() print("B2. Sensitivity of the answer to the temperatures themselves.") print(" Voyager's number carries +/- %.0f K, which is the wide one." % T_VOYAGER_ERR) print() print(" %-24s %-14s %-12s" % ("case", "ratio", "factor")) for lab, Ta, Tb in (("750 -> 358 (nominal)", T_VOYAGER, T_JWST), ("600 -> 358 (cool 1989)", 600.0, T_JWST), ("900 -> 358 (hot 1989)", 900.0, T_JWST), ("750 -> 350", T_VOYAGER, 350.0), ("750 -> 366", T_VOYAGER, 366.0)): rr = band_emission(Tb) / band_emission(Ta) print(" %-24s %-14s %-12.1f" % (lab, "%.4f%%" % (100.0 * rr), 1.0 / rr)) print() # ---------------------------------------------------------------- MODEL C print() print("MODEL C -- THE SENSITIVITY LADDER (our crudest model, read the") print(" docstring before believing any of it)") line() f_jwst = fom(6.5, 3420.0, 1.0) k_cal = 0.5 * r358 * f_jwst # JWST floor = half the measured brightness print(" Calibration: JWST's floor is set to half the brightness the") print(" paper actually measured, %.4f%% of the 750 K expectation." % (100.0 * 0.5 * r358)) print(" That fixes k = %.4f. Everything else is then predicted." % k_cal) print() print(" %-32s %-9s %-10s %-12s" % ("instrument", "FOM", "floor", "floor as %")) floors = {} for name, dia, t, bg in INSTRUMENTS: f = fom(dia, t, bg) floor = k_cal / f floors[name] = floor print(" %-32s %-9.1f %-10.2e %-12s" % (name, f, floor, "%.3f%%" % (100.0 * floor))) print() print(" Now run each one against a hot Neptune and a cold one.") print(" Emission is in units of the 750 K expectation, so a hot Neptune") print(" sits at 1.000 and the Neptune JWST actually saw sits at %.5f." % r358) print() print(" %-32s %-16s %-16s" % ("instrument", "at 750 K", "at 358 K")) for name, dia, t, bg in INSTRUMENTS: fl = floors[name] a = "DETECT" if 1.0 > fl else "miss" b = "DETECT" if r358 > fl else "miss" print(" %-32s %-16s %-16s" % (name, "%s (%.0fx floor)" % (a, 1.0 / fl), "%s (%.2fx floor)" % (b, r358 / fl))) print() print(" Read that table honestly. Three of the four rows are the story:") print(" every ground-based attempt had the sensitivity to find H3+ on a") print(" 750 K Neptune and none of them had the sensitivity to find it on") print(" a 358 K one. The fourth row is where our ladder breaks. It says") print(" Voyager's IRIS should have seen the ion during the flyby, and no") print(" such detection was ever reported. IRIS was a low-resolution") print(" interferometer whose H3+ lines would have been blended into the") print(" continuum, which our figure of merit knows nothing about. We are") print(" leaving the wrong answer in the table rather than tuning it out.") print() print(" The exposure penalty. Background-limited signal to noise scales") print(" as brightness times sqrt(area x time), so recovering a factor of") print(" %.1f in brightness costs a factor of %.0f in area x time." % (1.0 / r358, (1.0 / r358) ** 2)) print(" Keck would have needed roughly %.0f nights of pure integration." % (((1.0 / r358) ** 2 * 7200.0) / (8.0 * 3600.0))) print() # --------------------------------------------------------- FIGURE SERIES print() print("FIGURE DATA (transcribed by hand into the article's SVGs)") line() print("F2. centred dipole: tilt, poleward, centre, equatorward, width [L=20]") for t, lo, mean, hi, w in mig: print(" %6.1f %8.2f %8.2f %8.2f %7.2f" % (t, lo, mean, hi, w)) print() print("F2b. with the 0.55 R_N offset across the axis, same columns") for t, lo, mean, hi, w in mig_off: print(" %6.1f %8.2f %8.2f %8.2f %7.2f" % (t, lo, mean, hi, w)) print() print("F2c. beta scan at full tilt: beta, poleward, centre, equatorward,") print(" width, overlap with the observed band, miss (all degrees)") for b, lo, mean, hi, w, ovl, miss in beta_scan: print(" %6.0f %8.2f %8.2f %8.2f %7.2f %7.2f %7.2f" % (b, lo, mean, hi, w, ovl, miss)) print() print("F2d. L shell, magnetic latitude, poleward, equatorward [tilt 46.9]") for L, ml, lo, hi in Lband: print(" %6d %8.2f %8.2f %8.2f" % (L, ml, lo, hi)) print() print("F3. oval shape at four tilts: azimuth-ordered (lon, lat) pairs,") print(" 24 points each, L = 20, centred dipole") for t in (0.0, 20.0, 35.0, TILT_DEG): la, lo = oval(t, 20.0, offset=0.0, n_az=24) print(" tilt %.1f" % t) print(" lon " + " ".join("%6.1f" % v for v in lo)) print(" lat " + " ".join("%6.1f" % v for v in la)) print(" and the same four with the 0.55 R_N offset held across the axis") for t in (0.0, 20.0, 35.0, TILT_DEG): la, lo = oval(t, 20.0, offset=OFFSET_RN, beta_deg=90.0, n_az=24) print(" tilt %.1f (offset, beta = 90)" % t) print(" lon " + " ".join("%6.1f" % v for v in lo)) print(" lat " + " ".join("%6.1f" % v for v in la)) print() print("F4. temperature, relative emission (750 K = 1)") for T in range(200, 1001, 25): print(" %5d %.6e" % (T, band_emission(float(T)) / e750)) print() print("F4b. Monte Carlo latitude histogram, 30 bins from -75 to -15 deg") hist, edges = np.histogram(mc, bins=30, range=(-75.0, -15.0)) for i in range(30): print(" %7.2f %7.2f %5d" % (edges[i], edges[i + 1], hist[i])) print() # ------------------------------------------------------------- SUMMARY print() print("=" * 76) print("THE FOUR NUMBERS THE ARTICLE QUOTES") print("=" * 76) print(" 1. Oval latitude at 46.9 deg tilt, L=20 : %.1f to %.1f deg" % (pred_lo, pred_hi)) print(" (centred; the offset sweep spans : %.1f to %.1f deg)" % (min(r[1] for r in beta_scan), max(r[3] for r in beta_scan))) print(" 2. Observed JWST auroral band : %.1f to %.1f deg" % (OBS_LAT_S, OBS_LAT_N)) print(" 3. H3+ emission, 750 K -> 358 K : %.4f%%, factor %.1f" % (100.0 * r358, 1.0 / r358)) print(" 4. Area x time penalty to recover it : %.0f" % ((1.0 / r358) ** 2)) print() print(" Every one of those is the club's own arithmetic on published") print(" constants. The measurements belong to Melin and colleagues [1].") print("=" * 76) if __name__ == '__main__': main()