Tearing - BUGFIX - Doppler-shift each GGJ surface into its own E×B frame in the coupled growth rates

src/Tearing/Dispersion/Coupled.jl

@@ -15,7 +15,7 @@
 #   det = mc(Q::ComplexF64)
 #
 # At each evaluation, for k = 1 .. msing_max, the inner-layer Δ is computed
-# at a Q rescaled by `tauk_k / tauk_ref` and offset by that surface's real
+# at a Q rescaled by `tauk_k / tauk_ref` and offset by that surface's
 # Doppler shift `q_shift`, then subtracted (with the dc offset) from the
 # diagonal of an `msing_max × msing_max` upper-left submatrix of `dp_matrix`.
 # The off-diagonal Δ' couplings are passed through unchanged.
@@ -38,7 +38,8 @@ A root of `mc` in the complex `Q` plane is a coupled tearing eigenvalue.
 `ratio_k = tauk_k/tauk_ref`: `Q` is defined as `tauk·ω`, so one shared physical
 ω reaches surface `k` as `Q_k = tauk_k·ω`, on the same scaling as its `Q_e`/`Q_i`.
 
-`q_shift_k` is the real offset carried on each `SurfaceCoupling` (zero by default).
+`q_shift_k` is the E×B Doppler offset carried on each `SurfaceCoupling` (zero by default),
+applied in the layer's own frequency variable (see `SurfaceCoupling`).
 """
 struct MultiSurfaceCoupling{V<:AbstractVector{<:SurfaceCoupling}}
     surfaces::V
@@ -96,7 +97,7 @@ function (mc::MultiSurfaceCoupling)(Q::Number)
     M = mc.dp_matrix[1:n, 1:n]
     @inbounds for k in 1:n
         sc = mc.surfaces[k]
-        Q_k = Qc * (sc.tauk / ref_tauk) + sc.q_shift
+        Q_k = layer_argument(sc, Qc * (sc.tauk / ref_tauk))
         # m×m scalar coupling: use only the tearing channel. The
         # interchange (Glasser-stabilization) channel is carried in the
         # full 4m×4m dispersion in `CoupledFullMatch.jl`; this reduced

src/Tearing/Dispersion/CoupledFullMatch.jl

@@ -89,9 +89,9 @@ end
 Construct the 4m × 4m dispersion matrix driver. `dp_raw` must be the
 2m × 2m matrix in side-major ordering (the `intr.delta_prime_raw`
 field populated by `ForceFreeStates.compute_delta_prime_matrix!` on the
-Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k + q_shift_k`, with
-`ratio_k` and the real E×B Doppler offset `q_shift_k` exactly as in the reduced
-`multi_surface_coupling`.
+Riccati path). Surface k's inner layer is evaluated at `Q·ratio_k` Doppler-shifted by its
+E×B offset `q_shift_k`, exactly as in the reduced `multi_surface_coupling`: added to Q for a
+SLAYER surface, and as `i·q_shift_k` to the growth rate s for a GGJ surface (see `SurfaceCoupling`).
 
 # Keyword arguments
 
@@ -151,7 +151,7 @@ function (mc::MultiSurfaceCouplingFull)(Q::Number)
 
         # Map the shared scanned Q onto this surface's normalization, then Doppler it into
         # the surface's E×B frame (GGJ carries tauk = 1, so the ratio is a no-op there).
-        Q_k = Qc * (sc.tauk / ref_tauk) + sc.q_shift
+        Q_k = layer_argument(sc, Qc * (sc.tauk / ref_tauk))
         resp = solve_inner(sc.model, sc.params, Q_k; mc.inner_kwargs...)
 
         # delta1 = interchange (parity +), delta2 = tearing (parity −); named

src/Tearing/Dispersion/GGJGrowthRates.jl

@@ -53,10 +53,12 @@ its own 2 × 2 block of the raw Δ'). The box `gamma_range` × `omega_range`
 winding number; see `find_growth_rates`. The result is in the units of the
 SLAYER extraction: `gamma_Hz` is γ [s⁻¹], `omega_Hz` is −ω [rad/s], and
 `Q_root` is the normalized `(−ω + iγ)/Q₀`. The reported root is the least
-stable zero in the box. Without drifts the GGJ roots come in ±ω pairs, so
-either member of a pair can be reported. GGJ surfaces carry no time
-normalization, so no E×B Doppler shift is applied and the roots are in the
-plasma frame.
+stable zero in the box. Each surface's layer is evaluated in its own E×B frame,
+`s_k = s + i·n·Ω_E,k` from the `omega_E` and `n` given to `surface_coupling`;
+the scanned s, and so ω, is then in the lab frame, and in the plasma frame when
+no surface rotates, so with rotation set `omega_range` must cover the lab-frame ω,
+shifted by about n·Ω_E from the plasma frame. Without drifts or rotation the GGJ roots come in ±ω pairs,
+so either member of a pair can be reported.
 
 `gamma_range` must lie at γ > 0: the GGJ layer is evaluated on the principal
 branch of a square root in s, with a cut on the negative real s-axis. The search

src/Tearing/Dispersion/SurfaceCoupling.jl

@@ -27,8 +27,10 @@ Per-surface dispersion data: `(model, params, dp_diag, dc, scale, tauk, q_shift)
 r(Q) = dp_diag - scale * solve_inner(model, params, Q + q_shift).tearing - dc
 ```
 
-`q_shift` is a real offset on the layer's Q argument, zero by default. The Tearing runner
-sets it to the surface's E×B Doppler shift.
+`q_shift` is the surface's E×B Doppler shift, a real frequency offset in the layer's own
+frequency variable, zero by default. SLAYER's argument is the normalized Q, so the shift is added
+to Q; GGJ's argument is the growth rate s = γ + iω, so the shift is added to ω and enters as
+`s + i·q_shift`. The Tearing runner sets it for SLAYER surfaces.
 
 A root of `sc` in the complex `Q` plane is a **tearing** eigenvalue at
 this surface in the *uncoupled* approximation (only the tearing channel
@@ -48,10 +50,14 @@ struct SurfaceCoupling{M<:InnerLayerModel,P}
 end
 
 function (sc::SurfaceCoupling)(Q::Number)
-    Δ = solve_inner(sc.model, sc.params, ComplexF64(Q) + sc.q_shift).tearing
+    Δ = solve_inner(sc.model, sc.params, layer_argument(sc, Q)).tearing
     return sc.dp_diag - sc.scale * Δ - sc.dc
 end
 
+# Doppler-shift a scanned value into the layer's own frequency variable: real for SLAYER's Q, imaginary for GGJ's s.
+layer_argument(sc::SurfaceCoupling, Q::Number) = ComplexF64(Q) + sc.q_shift
+layer_argument(sc::SurfaceCoupling{<:GGJModel}, s::Number) = ComplexF64(s) + im * sc.q_shift
+
 """
     surface_coupling(model::SLAYERModel, params::SLAYERParameters,
                      dp_diag::Number; dc::Real=0.0, q_shift::Real=0.0)
@@ -74,24 +80,34 @@ end
 
 """
     surface_coupling(model::GGJModel, params::GGJParameters,
-                     dp_diag::Number) -> SurfaceCoupling
+                     dp_diag::Number; omega_E::Real=0.0, n::Integer=0) -> SurfaceCoupling
 
 GGJ convenience constructor. `scale` is `1.0` because GGJ's `solve_inner`
 applies its own `rescale_delta` (S^(2p₁/3)·v1^(2p₁)) internally, so the
 returned Δ is already in outer units. `tauk` defaults to `1.0` (GGJ has no
 direct analogue of SLAYER's per-surface time normalization, so multi-surface
 Q rescaling is a no-op for GGJ surfaces unless overridden).
 
+`omega_E` is the surface's E×B rotation per unit toroidal mode number [rad/s] and `n` the
+toroidal mode number. GGJ is a plasma-frame layer: rotation enters only beyond the order its
+equations retain (Johnson & Greene 1967, Plasma Phys. 9, 611, Appendix B), so a rotating
+surface sees the lab-frame growth rate Doppler-shifted, `s_k = s + i·n·omega_E`, and
+`q_shift = n·omega_E`. With rotation set on every surface, the scanned `s` is in the lab
+frame; differential rotation between coupled surfaces changes the coupled γ, while a common
+rotation only shifts ω.
+
 **No `dc` kwarg**: GGJ's 4m×4m Pletzer-Dewar residual already includes the
 interchange channel, which provides Glasser (Mercier) stabilization
 natively. A Δ_crit proxy (χ_parallel-matching offset on the diagonal) is
 meaningful only for tearing-only slab-layer approximations like SLAYER;
 for GGJ it would double-count the interchange physics. The `SurfaceCoupling`
-struct's `dc` field is hard-wired to 0 here, and so is `q_shift`.
+struct's `dc` field is hard-wired to 0 here.
 """
-function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number)
+function surface_coupling(model::GGJModel, params::GGJParameters, dp_diag::Number; omega_E::Real=0.0, n::Integer=0)
+    iszero(omega_E) || n >= 1 ||
+        throw(ArgumentError("surface_coupling: E×B rotation omega_E=$(omega_E) rad/s needs the toroidal mode number n ≥ 1; got n=$(n)"))
     return SurfaceCoupling(model, params, ComplexF64(dp_diag),
-        0.0, 1.0, 1.0, 0.0)
+        0.0, 1.0, 1.0, Float64(n * omega_E))
 end
 
 """

src/Tearing/Runner/run_slayer.jl

@@ -227,8 +227,6 @@ function _q_shift(p::SLAYERParameters, Ω_E::Real)
         throw(ArgumentError("run_slayer: E×B rotation Ω_E=$(Ω_E) rad/s needs the toroidal mode number, but SLAYERParameters has n=$(p.n)"))
     return -p.tauk * p.n * Ω_E
 end
-# GGJ carries no time normalization (tauk = 1), so a physical rotation has no Q-space image.
-_q_shift(::InnerLayerParameters, ::Real) = 0.0
 
 # ---------------------------------------------------------------------
 # Reference-length conversion of the outer Δ' for the slab layer

test/runtests_dispersion_ggj.jl

@@ -109,17 +109,35 @@ GeneralizedPerturbedEquilibrium.Dispersion._residual_domain_radius(f::GGJDomainT
         mc2 = multi_surface_coupling_full([surface_coupling(model, p, r0.tearing), surface_coupling(model, p2, r02.tearing)], dp2)
         g2 = GGJDispersion(mc2)
         # Independent reference: a complex secant iteration on the determinant, started at the uncoupled root.
-        a, b = 0.1234im, 0.1234im * 1.01
-        fa, fb = g2(a), g2(b)
-        for _ in 1:40
-            a, fa, b = b, fb, b - fb * (b - a) / (fb - fa)
-            fb = g2(b)
-            abs(b - a) < 1e-14 && break
+        function secant_root(g, x0)
+            x1 = 1.01 * x0
+            f0, f1 = g(x0), g(x1)
+            for _ in 1:60
+                x0, f0, x1 = x1, f1, x1 - f1 * (x1 - x0) / (f1 - f0)
+                f1 = g(x1)
+                abs(x1 - x0) < 1e-14 && break
+            end
+            return x1
         end
+        b = secant_root(g2, 0.1234im)
         @test abs(b - 0.1234im) > 1e-5                 # the coupling moved the root
         gr2, _ = ggj_growth_rates(mc2; gamma_range=(0.08Q0, 0.17Q0), omega_range=(-0.031Q0, 0.047Q0),
             nre0=6, nim0=6, passes=3, parallel=false)
         @test length(gr2.valid_roots) == 1 && all(w -> w ≈ 1, gr2.valid_roots_winding)
         @test abs(gr2.Q_root - b) < 1e-8 * abs(b)
+
+        @testset "E×B rotation Doppler-shifts each surface's s" begin
+            # Rotation enters the plasma-frame GGJ layer as s_k = s + i·n·Ω_E,k on each surface.
+            mc_rot(Ω1, Ω2) = multi_surface_coupling_full([surface_coupling(model, p, r0.tearing; omega_E=Ω1, n=1),
+                    surface_coupling(model, p2, r02.tearing; omega_E=Ω2, n=1)], dp2)
+            Ω = 0.02 * Q0
+            @test_throws ArgumentError surface_coupling(model, p, r0.tearing; omega_E=Ω)
+            # A common rotation is a change of frame: γ is unchanged and the lab ω shifts by n·Ω (Re Q by Ω/Q₀, as in SLAYER).
+            b_common = secant_root(GGJDispersion(mc_rot(Ω, Ω)), b + Ω / Q0)
+            @test abs(b_common - (b + Ω / Q0)) < 1e-8 * abs(b)
+            # Differential rotation moves the coupled growth rate.
+            b_diff = secant_root(GGJDispersion(mc_rot(0.0, 0.5 * Q0)), b)
+            @test abs(imag(b_diff) - imag(b)) > 1e-4 * imag(b)
+        end
     end
 end