Add magnon-dispersion-bec task

.github/scripts/validate_repository.py

@@ -6,6 +6,7 @@
 ROOT = Path(__file__).resolve().parents[2]
 EXPECTED_TASKS = (
     "heterodyne-shot-noise-analysis",
+    "magnon-dispersion-bec",
     "surface-ion-trap-shuttling",
     "trapped-ions-heating-rate",
 )

registry.json

@@ -12,6 +12,11 @@
         "path": "tasks/heterodyne-shot-noise-analysis",
         "digest": "sha256:497468dab7ed79dca0f5b5768eb5b854b0fdd420247952f49c741f63d6ae27f2"
       },
+      {
+        "name": "magnon-dispersion-bec",
+        "path": "tasks/magnon-dispersion-bec",
+        "digest": "sha256:9ff077258443093a1c471d49abfd3685b53494133400b5b4ff4f0c579f661938"
+      },
       {
         "name": "surface-ion-trap-shuttling",
         "path": "tasks/surface-ion-trap-shuttling",

tasks/magnon-dispersion-bec/environment/Dockerfile

@@ -0,0 +1,25 @@
+FROM python:3.12-slim
+
+ENV DEBIAN_FRONTEND=noninteractive \
+    OMP_NUM_THREADS=1 \
+    OPENBLAS_NUM_THREADS=1 \
+    MKL_NUM_THREADS=1 \
+    MPLBACKEND=Agg
+
+RUN pip install --no-cache-dir \
+    numpy==2.2.6 \
+    scipy==1.15.3 \
+    matplotlib==3.10.3 \
+    pytest==8.4.1 \
+    pytest-json-ctrf==0.3.5
+
+WORKDIR /app
+
+RUN mkdir -p \
+    /app \
+    /home/agent/.claude \
+    /home/agent/.codex \
+    /home/agent/.opencode \
+    /home/agent/.agents && \
+    chmod -R 777 /home/agent && \
+    chmod 777 /root

tasks/magnon-dispersion-bec/environment/skills/magnon-gas-thermodynamics/SKILL.md

@@ -0,0 +1,66 @@
+---
+name: magnon-gas-thermodynamics
+description: Reusable guidance for YIG-film spin-wave calculations and equilibrium or pumped magnon-gas thermodynamics.
+---
+
+# YIG spin waves and magnon gases
+
+Use this guide when a problem combines finite-thickness ferromagnetic-film
+spectra with a thermal or pumped magnon gas. Keep the field, wavevector,
+thickness, frequency, and energy units consistent before doing any fitting or
+integration.
+
+## Film spectra
+
+- Use the full finite-thickness dipole-exchange dispersion for an
+  in-plane-magnetized film. With `theta` measured from the magnetization,
+  `F(x) = (1 - exp(-x))/x` and
+
+  `f(k, theta) = gamma * sqrt((H + D*k^2 + 4piM*(1 - F(k*d))*sin(theta)^2) * (H + D*k^2 + 4piM*F(k*d)))`.
+
+  Take `F(0) = 1`, keep the propagation angle, and use the units supplied by
+  the problem. An exchange-only or bulk approximation can give the wrong
+  branch and minimum.
+- At zero wavevector, use the Kittel relation to infer the saturation
+  magnetization when an FMR frequency is supplied.
+- Infer exchange stiffness from an observed finite-wavevector minimum by
+  minimizing the backward-volume branch while fitting the stiffness. Recompute
+  the minimum after changing only the thickness when comparing films of the
+  same material.
+- Evaluate surface or Damon-Eshbach predictions with the same material
+  parameters and full form factor. Use a bounded minimizer and check that the
+  result is stable when the wavevector grid or bracket changes.
+
+## Thermal magnons
+
+- Decide whether the thermal cloud is a 3D bulk gas or a lower-dimensional
+  film gas from the physical wavelength and geometry. For an exchange-dominated
+  3D cloud, use the isotropic dispersion
+  `epsilon(k) = h * gamma * (H + D*k^2)`.
+- Compute number and energy from the Bose-Einstein occupation with a
+  three-dimensional density-of-states factor. Convert frequency to joules with
+  `h`, and keep the integration wavevector in the units used by `D`.
+- Check the spectral gap at the lowest energy, the density scaling with
+  temperature, quadrature convergence, and the high-wavevector cutoff.
+
+## Pumped quasi-equilibrium
+
+- Rapid thermalization with conserved magnon number and energy gives two
+  conservation equations for the pumped temperature and chemical potential.
+  Include the injected energy at the pump frequency.
+- Below condensation, solve for `(T, mu)` while enforcing `mu` below the
+  spectral minimum. The critical injected density is the increase needed to
+  reach that minimum at the reference temperature.
+- Above threshold, pin `mu` at the minimum and solve for the thermal
+  temperature plus condensate density. Count condensate magnons at the minimum
+  energy when enforcing energy conservation, then divide by total density for
+  the condensate fraction.
+- If one configuration feeds a fitted material parameter into another, pass
+  the computed value through rather than substituting a nominal material value.
+
+## Numerical checks
+
+Use SI constants for energies, document any unit conversions, bracket roots
+where possible, and verify both number and energy conservation from the final
+state. Keep scripts and intermediate calculations general; do not hardcode
+expected outputs or task-specific field names.

tasks/magnon-dispersion-bec/environment/skills/magnon-gas-thermodynamics/references/model.md

@@ -0,0 +1,40 @@
+# Units and worked checks
+
+## Constants
+- h = 6.62607015e-34 J s, kB = 1.380649e-23 J/K
+- gamma = 2.8e-3 GHz/Gauss = 2.8e6 Hz/Gauss (YIG)
+
+## Part A units
+- H and 4piM in Gauss (Oe); D in Oe cm^2 with k in cm^-1, so `D k^2` is in Gauss.
+- `omega(k,theta)` above is in GHz. f_FMR, f_min in GHz.
+- The dipolar form factor F(k) is dimensionless; `k d` uses d in cm.
+- Sanity: with H around a few hundred Gauss and 4piM ~ 1.8e3 Gauss, f_FMR is a
+  few GHz; the backward-volume minimum sits a few hundred MHz below f_FMR at
+  k0 ~ 1e5 cm^-1.
+
+## Part B units
+- Bulk magnon energy `eps(k) = h * gamma * (H + D k^2) * 1e9` in Joules
+  (the `1e9` converts GHz -> Hz before multiplying by h).
+- Number density `n = (1/(2 pi^2)) int_0^inf k^2 nBE dk` with k in cm^-1 gives
+  cm^-3 directly (the 3D density of states `d^3k/(2 pi)^3` integrated over solid
+  angle gives `k^2/(2 pi^2) dk`).
+- Energy density e in J/cm^3.
+- Bose factor uses `expm1` for numerical stability: `1/np.expm1(x)`.
+
+## Worked check (arbitrary parameters, not the task values)
+For H = 600 Gauss, D = 2e-9, T = 300 K, mu = 0:
+- eps(0)/kB = h*gamma*H*1e9/kB should be ~ 81 mK.
+- the number integral is dominated by k ~ 1e7 cm^-1 (where eps ~ kB T) and the
+  density comes out of order 1e21 cm^-3, following n(T) proportional to T^1.5
+  (so the density-vs-temperature exponent is ~1.5, not ~1).
+
+## Supercritical / condensed branch (units)
+- With `mu` pinned at `eps_min`, the thermal-cloud integrals `n(T, eps_min)` and
+  `e(T, eps_min)` are evaluated exactly as above but at `mu = eps_min` (use
+  `mu = eps_min (1 - tiny)` to avoid the k=0 singularity of the Bose factor).
+- The condensate density `n_c` (cm^-3) is a number, sitting at energy `eps_min`
+  per magnon; it carries energy `n_c eps_min` (J/cm^3).
+- Condensate fraction is dimensionless: `n_c / (n(T0,0) + delta_N)`.
+
+Use these to validate your dispersion and integration before computing the task
+values. Do NOT hardcode any result; derive everything from the stated inputs.

tasks/magnon-dispersion-bec/oracle/magnon_physics.py

@@ -0,0 +1,207 @@
+#!/usr/bin/env python3
+"""Physics core for the YIG magnon dispersion + pumped-magnon-gas task (oracle copy).
+
+Part A -- dipole-exchange spin-wave dispersion of a magnetized YIG film
+(Rezende, PRB 79, 174411 (2009)); extract the saturation magnetization 4piM
+(Kittel) and the exchange stiffness D (from the spectral minimum f_min), then
+predict the surface (Damon-Eshbach) branch and a second film's spectral minimum.
+
+Part B -- thermodynamics of the parametrically pumped magnon gas following
+Demokritov et al., Nature 443, 430 (2006): the exchange-dominated 3D bulk
+thermal magnon gas (Bloch T^3/2 law), the chemical potential / temperature of
+the pumped quasi-equilibrium gas from magnon-number and energy conservation,
+the critical pumped density for Bose-Einstein condensation, and the condensed
+branch (condensate fraction) of a supercritically pumped gas.
+"""
+
+from __future__ import annotations
+
+import numpy as np
+from scipy.integrate import quad
+from scipy.optimize import brentq, fsolve
+
+H_PLANCK = 6.62607015e-34  # J s
+KB = 1.380649e-23  # J / K
+GAMMA = 2.8e-3  # GHz / Gauss  (gyromagnetic ratio for YIG)
+
+# ============================ PART A ============================
+# NV3 experimental inputs (film 1)
+NV3 = {"H": 330.0, "f_FMR": 2.395, "f_min": 1.913, "d": 0.206e-4}  # d in cm (0.206 um)
+D2_FILM = 0.5e-4  # second film thickness (cm), same material, same field H
+K_S = 8.0e4  # surface-branch (theta = pi/2) probe wavevector (cm^-1) for film 1
+
+
+def _F(k, d):
+    k = np.asarray(k, float)
+    return np.where(k != 0, (1 - np.exp(-k * d)) / (k * d), 1.0)
+
+
+def film_dispersion_GHz(k, theta, H, fourpiM, D, d):
+    """Rezende dipole-exchange film magnon frequency (GHz). k in cm^-1."""
+    a = H + D * k**2 + fourpiM * (1 - _F(k, d)) * np.sin(theta) ** 2
+    b = H + D * k**2 + fourpiM * _F(k, d)
+    return GAMMA * np.sqrt(a * b)
+
+
+def extract_fourpiM(f_FMR, H):
+    """Kittel: f_FMR = gamma sqrt(H (H + 4piM))  ->  4piM."""
+    return f_FMR**2 / (GAMMA**2 * H) - H
+
+
+def film_fmin(H, fourpiM, D, d, n=400000):
+    ks = np.linspace(1.0, 3e5, n)
+    fk = film_dispersion_GHz(ks, 0.0, H, fourpiM, D, d)
+    i = int(np.argmin(fk))
+    return float(fk[i]), float(ks[i])
+
+
+def extract_D(f_min_meas, H, fourpiM, d):
+    """Fit D so the backward-volume (theta=0) spectral minimum equals f_min."""
+    return brentq(lambda D: film_fmin(H, fourpiM, D, d)[0] - f_min_meas, 1e-10, 3e-8)
+
+
+def part_A():
+    p = NV3
+    fourpiM = extract_fourpiM(p["f_FMR"], p["H"])
+    D = extract_D(p["f_min"], p["H"], fourpiM, p["d"])
+    fmin, kmin = film_fmin(p["H"], fourpiM, D, p["d"])
+    # surface (Damon-Eshbach, theta = pi/2) branch of film 1 at K_S
+    f_DE = float(film_dispersion_GHz(K_S, np.pi / 2, p["H"], fourpiM, D, p["d"]))
+    # second film (same material -> same 4piM, D) spectral minimum at thickness D2_FILM
+    fmin2, kmin2 = film_fmin(p["H"], fourpiM, D, D2_FILM)
+    return {
+        "fourpiM_G": fourpiM,
+        "D_Oecm2": D,
+        "f_min_GHz": fmin,
+        "k_min_percm": kmin,
+        "f_DE_GHz": f_DE,
+        "fmin2_GHz": fmin2,
+        "kmin2_percm": kmin2,
+    }
+
+
+# ============================ PART B ============================
+# Config B1: Nature-06 YIG system + sub-critical pumping
+NB = {"H": 700.0, "D": 2e-9}  # Gauss, Oe cm^2
+T0 = 300.0  # K
+NU_P = 4.0e9  # pumped magnon frequency (Hz)
+DELTA_N = 5e18  # pumped magnon density (cm^-3)
+
+# Config B2: a second YIG system, supercritically pumped (condenses)
+NB2 = {"H": 500.0, "D": 2e-9}
+T0_2 = 300.0
+NU_P2 = 4.0e9
+DELTA_N2 = 1.5e20  # cm^-3 (chosen well above the B2 critical density)
+
+# Config B3: interlocked system that uses the exchange stiffness extracted in Part A
+NB3 = {"H": 400.0}  # D supplied from Part A at runtime
+T0_3 = 300.0
+DELTA_N3 = 2e18  # cm^-3
+
+
+def _eps(k, H, D):
+    """3D bulk exchange magnon energy (J). k in cm^-1."""
+    return H_PLANCK * GAMMA * (H + D * k**2) * 1e9
+
+
+def bulk_density(T, mu, H, D):
+    """3D bulk thermal magnon density (cm^-3): n = (1/2pi^2) int k^2 nBE dk."""
+    val, _ = quad(lambda k: k**2 / np.expm1((_eps(k, H, D) - mu) / (KB * T)), 0, 3e8, limit=600)
+    return val / (2 * np.pi**2)
+
+
+def bulk_energy(T, mu, H, D):
+    val, _ = quad(lambda k: k**2 * _eps(k, H, D) / np.expm1((_eps(k, H, D) - mu) / (KB * T)), 0, 3e8, limit=600)
+    return val / (2 * np.pi**2)
+
+
+def _subcritical_solve(H, D, T0_cfg, nu_p, delta_n, n0, e0):
+    """Pumped sub-critical quasi-equilibrium: conserve number and energy -> (T, mu)."""
+    nT = n0 + delta_n
+    eT = e0 + H_PLANCK * nu_p * delta_n
+    eps_min = _eps(0.0, H, D)
+    sol = fsolve(
+        lambda v: [(bulk_density(v[0], v[1], H, D) - nT) / nT, (bulk_energy(v[0], v[1], H, D) - eT) / eT],
+        [T0_cfg, eps_min * 0.3],
+        full_output=True,
+    )
+    (T_p, mu_p), ier = sol[0], sol[2]
+    return T_p, mu_p, bool(ier == 1)
+
+
+def _condensed_solve(H, D, T0_cfg, nu_p, delta_n, n0, e0):
+    """Pumped supercritical branch: mu pinned at eps_min, conserve number+energy
+    for (T, condensate density n_c). Returns (T_p, condensate_fraction)."""
+    eps_min = _eps(0.0, H, D)
+    mu_pin = eps_min * (1 - 1e-12)
+    n_tot = n0 + delta_n
+    e_tot = e0 + H_PLANCK * nu_p * delta_n
+
+    def f_T(T):
+        n_th = bulk_density(T, mu_pin, H, D)
+        return bulk_energy(T, mu_pin, H, D) + (n_tot - n_th) * eps_min - e_tot
+
+    T_p = brentq(f_T, 20.0, 3.0 * T0_cfg)
+    n_c = n_tot - bulk_density(T_p, mu_pin, H, D)
+    return T_p, n_c / n_tot
+
+
+def part_B(D_from_A=None):
+    # ---- Config B1 (sub-critical) ----
+    H, D = NB["H"], NB["D"]
+    eps_min = _eps(0.0, H, D)  # k=0 gap = h gamma H
+    n0 = bulk_density(T0, 0.0, H, D)
+    e0 = bulk_energy(T0, 0.0, H, D)
+    n_crit = bulk_density(T0, eps_min * (1 - 1e-9), H, D)  # density at mu -> eps_min
+    dN_c = n_crit - n0
+    T_p, mu_p, conv1 = _subcritical_solve(H, D, T0, NU_P, DELTA_N, n0, e0)
+    # Bloch exponent of the gapless 3D gas: n0(T) ~ T^p near T0
+    p_bloch = float(np.log(bulk_density(2 * T0, 0.0, H, D) / n0) / np.log(2.0))
+
+    # ---- Config B2 (supercritical -> condensate) ----
+    H2, D2 = NB2["H"], NB2["D"]
+    eps_min2 = _eps(0.0, H2, D2)
+    n0_2 = bulk_density(T0_2, 0.0, H2, D2)
+    e0_2 = bulk_energy(T0_2, 0.0, H2, D2)
+    dN_c2 = bulk_density(T0_2, eps_min2 * (1 - 1e-9), H2, D2) - n0_2
+    bec2 = bool(DELTA_N2 > dN_c2)
+    T_p2, ncfrac2 = _condensed_solve(H2, D2, T0_2, NU_P2, DELTA_N2, n0_2, e0_2)
+
+    # ---- Config B3 (A->B interlock: uses D extracted in Part A) ----
+    D3 = D_from_A if D_from_A is not None else part_A()["D_Oecm2"]
+    H3 = NB3["H"]
+    eps_min3 = _eps(0.0, H3, D3)
+    n0_3 = bulk_density(T0_3, 0.0, H3, D3)
+    dN_c3 = bulk_density(T0_3, eps_min3 * (1 - 1e-9), H3, D3) - n0_3
+    bec3 = bool(DELTA_N3 > dN_c3)
+
+    return {
+        # config B1
+        "eps_min_mK": eps_min / KB * 1e3,
+        "n0_percm3": n0,
+        "e0_Jpercm3": e0,
+        "p_bloch": p_bloch,
+        "dN_c_percm3": dN_c,
+        "mu_p_mK": mu_p / KB * 1e3,
+        "T_p_K": T_p,
+        "bec": bool(DELTA_N > dN_c),
+        "converged": conv1,
+        # config B2 (supercritical)
+        "eps_min2_mK": eps_min2 / KB * 1e3,
+        "dN_c2_percm3": dN_c2,
+        "bec2": bec2,
+        "ncfrac2": ncfrac2,
+        "T_p2_K": T_p2,
+        # config B3 (interlock)
+        "n0_3_percm3": n0_3,
+        "dN_c3_percm3": dN_c3,
+        "bec3": bec3,
+    }
+
+
+if __name__ == "__main__":
+    import json
+
+    A = part_A()
+    B = part_B(D_from_A=A["D_Oecm2"])
+    print(json.dumps({"A": A, "B": B}, indent=2))