[matt-pharr]

Release note

• Audience: users
• Numerical impact: none for existing models; adds a new layer model (harness @ b84b559 vs InnerLayer - REFACTOR - Move the GGJ solvers into a shared ToroidalLayer library behind a layer-model interface #653, pending)
• Migration: none

Adds CCCModel, the drift-resistive linear MHD inner layer of Cho, Choi & Choi (1987) in the form of Kim, Cho & Choi (1988): the GGJ layer with electron diamagnetic drift (Q*) and perpendicular transport (M⊥/M), solved with the rotated-ray backend.

Wasow analysis was done in mathematica to obtain the asymptotic series using a wolframscript routine I have been using for my other derivation which essentially just re-implemented Zhirui's mathematica notebook from GW2020. This is essentially the same as the GGJ ray solver but with coefficients and asymptotic series for the Cho-Choi-Choi layer, so this is a small, drop-in change.

Regression report

Validation and benchmarks

These runs used the prototype of this model (the same equations threaded through GGJ on an earlier branch); CCCModel reproduces it to 1.5e-8 at the pinned point (q4 surface, Q* = 0.5, M⊥/M = 1.4, Q = 5i).

Kim, Cho & Choi (1988) Figs. 1b and 2b. Their surface (E = −0.01, F = H = 0, G = K = 100, M⊥/M = 1), tearing (odd) parity, Q* = 0 (Fig. 1b) and Q* = 0.05 (Fig. 2b). Roots of Δ₋(Q) = Δ₀ found by secant continuation along Δ₀, overlaid on the published curves with the scan axes calibrated from their tick marks. Kim et al. plot their Fourier-space tearing quantity Δ_o, which their Eq. (10b) relates to the real-space Δ′_o we compute; we convert with that equation and fit nothing.

• The Δ₀ = 0 roots, which do not depend on how Δ is normalized, land on the published curves:

  	ours	read from Kim et al.
  Fig. 1b, Q* = 0	−0.0187 ± 0.0324i	−0.018 ± 0.033i
  Fig. 2b, Q* = 0.05, branch A	−0.0184 − 0.0626i	−0.019 − 0.062i
  Fig. 2b, Q* = 0.05, branch B	−0.0152 + 0.0097i	−0.014 + 0.010i

• Away from Δ₀ = 0, Kim's axis is Δ_o = C_o·Δ′_o with C_o = |Γ(1+p₊)/Γ(−p₊)·cot(πp₊/2)| and p₊ = −1/2 + √(−D_I) (Eq. 10b); for this surface C_o = 0.6293327. With that conversion our curves overlay all of Fig. 1b and all four Q* = 0.05 curves of Fig. 2b. C_o depends only on the Mercier exponent, so the same value serves both Q*. Our Δ uses the GGJ normalization cross-checked against Fortran rmatch.

Convergence to GGJ. On the q4 benchmark surface, Q* from 1e-7 to 1 and M⊥/M − 1 from 1e-7 to 1, at Q = 0.01, 0.1 + 0.1i, i, 3 + 3i, 10i and 100i:

• The departure from GGJ is linear in Q* (fitted slope 1.00) at every Q down to the solver's accuracy (1e-7 to 1e-8), with no singular behaviour where the two exponential families merge.
• The size of the drift effect falls with |Q|: at Q* = 1e-2, about 3e-2 relative at |Q| ≈ 0.1 and 4e-5 at Q = 100i.
• M⊥/M enters linearly at small |Q| and is below 1e-5 for |Q| ≳ 3 even at M⊥/M = 2.

Accuracy and cost.

• delta_convergence spread (contour angle, mesh, series order and radius, march tolerance, BVP edge, decaying-pair separation) on the q4 surface with Q* = 0.5, M⊥/M = 1.4: ≤ 3.4e-7 at Q = 0.5 + 0.5i and ≤ 7.8e-8 at Q = 5i.

• Single-solve times on one laptop core:

  case	time
  q4, Q* = 0.5, M⊥/M = 1.4, Q = 5i	0.3 s
  q4, Q* = 2, Q = 5i	5 s
  q4, Q* = 2, Q = 20 + 20i	18 s (GGJ at the same Q: 0.6 s)

  The cost grows with K·|Q|·Q*², because the march resolves the faster exponential family across the slower one's whole band.

Tearing roots on GPEC examples (dispersion roots by continuation from the SLAYER root, Q* = −ω*ₑ/Q₀ from the kinetic profiles, M⊥/M = 1; not yet wired into TearingProblem):

• DIII-D-like SLAYER example (examples/DIIID-like_SLAYER_example, Q* ≈ 1.2–1.5 at 2/1–4/1). The CCC roots rotate at SLAYER's lab frequency on 2/1 and 3/1 (within 3%) and near it on 4/1 (−3.6e3 against −2.8e3 rad/s). The damping on 3/1 and 4/1 comes out at −310 and −420 s⁻¹, against SLAYER's −145 and −450 s⁻¹, while plain GGJ gives −1230 and −1756 s⁻¹. 2/1 differs: SLAYER is unstable at +164 s⁻¹ (without its curvature offset), CCC damped at about −330 s⁻¹ (not converged to tolerance), GGJ −1004 s⁻¹.
• TokaMaker aspect-ratio scan (ε = 0.05–0.75, β_N ≈ 0.1, κ and δ rising with ε; synthetic kinetic profiles with T_e = T_i = 2 keV and n from the equilibrium pressure, Ω_E = 0). For the 2/1:
  • CCC tracks SLAYER with its toroidal curvature offset to 10–20% in γ and a few % in ω for 0.11 < ε < 0.4, and both stabilize near ε ≈ 0.45–0.5;
  • SLAYER without the offset stays unstable across the scan;
  • plain GGJ has no unstable 2/1 root above ε ≈ 0.25, while CCC stays unstable to ε ≈ 0.44. This is Kim et al.'s finding that diamagnetic drift weakens the favourable-curvature stabilization.

Notes for reviewers

• Stacked on InnerLayer - REFACTOR - Move the GGJ solvers into a shared ToroidalLayer library behind a layer-model interface #653 (the ToroidalLayer library); review only the commits above 552e97504.
• CCCModel / CCCParameters (src/InnerLayer/CCC/) implement the ToroidalLayerParameters interface: Kim (1988) Eqs. (1a)–(1c) as rows and weak form, a Wasow basis from Kim's A₀–A₂ with a block-diagonalizing T that equals GW2020's at Q* = 0, M⊥/M = 1, a general Sylvester splitting solve, and the two exponential families of Kim Eq. (2) setting the ray angle and mesh. CCCParameters holds a GGJParameters plus Qstar, mperp, and takes GGJ's scales and Δ normalization. Only :ray is available.
• The library gains two optional hooks, parabolic_resolution_rate and ray_angle_admissible, whose defaults leave GGJ unchanged.
• Validation (unit tests in runtests_innerlayer.jl): the Wasow matrices reproduce the ODE matrix to 1e-16; the first-order Ξ series coefficient matches a closed form from a separate Mathematica derivation to 1e-8; Δ is independent of the contour angle to 1e-7; CCCModel at Q* = 0, M⊥/M = 1 reproduces GGJ to 3e-10; Kim Figs. 1b and 2b (tearing dispersion at Q* = 0 and 0.05) are reproduced with no fitted parameter once our real-space Δ′_o is converted to their Fourier-space Δ_o (Kim Eq. 10b; factor 0.6293 for this surface).
• Still to come on this branch: TearingProblem, dispersion, HDF5 output, deck (inner_model = "CCC") and MatchProblem support through ToroidalLayerModel, with Q* = −ω*ₑ/Q₀ from the kinetic profiles, a regression case and deck/API parity tests.
• Known limitation: solves are slow when K·Q·Q² is large (the march resolves the faster exponential family across the slower one's band); at DIII-D-like Q ≈ 1 an adaptive dispersion scan can take hours.

---

🤖 (Partially :P) generated with Claude Code

[github-actions]

This pull request is missing an assignee.

If you are not ready to name them, mark this pull request as a draft.
docs/development/contributors.md suggests lead developers to ask.
Merging is not blocked here, but no pull request may be merged without human review.

[logan-nc]

Awesome!