diff --git a/CHANGELOG.md b/CHANGELOG.md
index 3dccea7e0c..95ba234182 100644
--- a/CHANGELOG.md
+++ b/CHANGELOG.md
@@ -3,6 +3,7 @@ Changelog
 
 New Features
 
+- Adds Available Energy for nonlinear measure of turbelent transport [#2215](https://github.com/PlasmaControl/DESC/pull/2215) and local potential well plotting utitlity for analysis.
 - Adds ``desc.objectives.DeflationOperator``, a new objective class which can be used to apply deflation techniques to equilibrium and optimization problems to find multiple local minima or multiple solutions from a single initial point, either by wrapping an existing ``desc.objectives._Objective`` object or by including as an additional penalty or constraint. Also adds a tutorial showing this functionality.
 - Sub-objectives of an `ObjectiveFunction` can now have different `use_jit` values than the `ObjectiveFunction`. These objectives have to be built before building the `ObjectiveFunction`.
 - Adds ``num_neighbors`` parameter to ``CoilSetMinDistance`` that limits the pairwise distance computation to the nearest neighbors per coil, reducing memory useage for large coilsets.
@@ -30,6 +31,7 @@ Performance Improvements
 - Reduces import time of `desc` modules.
     - Now, `desc.compute._build_data_index` uses depth-first search algorithm to construct the dependency tree.
     - Some of the default value computations at import time are removed (i.e. `desc.integrals.bounce_integral.default_quad`)
+- Significantly improves convergence of bounce point computations. Now lower ``Y_B`` can be used.
 - [Significantly improves convergence of inverse stream maps](https://github.com/PlasmaControl/DESC/pull/1919).
 - Resolves a JAX memory regression in bounce integrals by avoiding materialization of a large tensor in memory. Previously, we had closed the issue by adding nuffts as a workaround. This update actually solves the issue for the case when a user specifies to not use nuffts as well.
 - ``ObjectiveFunction.print_value`` can now use the previously computed ``compute_scaled_error`` values to print. For bounded objectives, we fall back to computing ``compute_unscaled``. Additionally, ``compute_scaled_error`` and array splitting are used in other parts of the code to prevent recompilation for one-time tasks, which makes initialization faster.
diff --git a/desc/compute/__init__.py b/desc/compute/__init__.py
index 343ca29002..916ecea52e 100644
--- a/desc/compute/__init__.py
+++ b/desc/compute/__init__.py
@@ -42,6 +42,7 @@
     _profiles,
     _stability,
     _surface,
+    _turbulence,
 )
 from .data_index import all_kwargs, allowed_kwargs, data_index
 from .utils import (
diff --git a/desc/compute/_fast_ion.py b/desc/compute/_fast_ion.py
index e6abd74e19..8186efd3c9 100644
--- a/desc/compute/_fast_ion.py
+++ b/desc/compute/_fast_ion.py
@@ -33,7 +33,7 @@ def _gamma_c_data(data):
     # units Tesla meters. Smoothness is determined by positive lower bound of
     # log argument, and hence behaves as ∂log(|B|/B₀)/∂ρ |B| = ∂|B|/∂ρ.
     return {
-        "|grad(psi)|*kappa_g": data["|grad(psi)|"] * data["kappa_g"],
+        "|grad(psi)|*kappa_g": data["|grad(psi)|*kappa_g"],
         "|grad(rho)|*|e_alpha|r,p|": data["|grad(rho)|"] * data["|e_alpha|r,p|"],
         "|B|_r|v,p": data["|B|_r|v,p"],
         "K": data["iota_r"]
@@ -48,7 +48,7 @@ def _v_tau(data, B, pitch):
     return safediv(2.0, jnp.sqrt(jnp.abs(1 - pitch * B)))
 
 
-def _radial_drift_1(data, B, pitch):
+def _radial_drift(data, B, pitch):
     return (
         safediv(1 - 0.5 * pitch * B, jnp.sqrt(jnp.abs(1 - pitch * B)))
         * data["|grad(psi)|*kappa_g"]
@@ -57,21 +57,18 @@ def _radial_drift_1(data, B, pitch):
 
 
 def _poloidal_drift_periodic(data, B, pitch):
-    return (
-        safediv(1 - 0.5 * pitch * B, jnp.sqrt(jnp.abs(1 - pitch * B)))
-        * data["|B|_r|v,p"]
-        + jnp.sqrt(jnp.abs(1 - pitch * B)) * data["K"]
-    ) / B
+    g = jnp.sqrt(jnp.abs(1 - pitch * B))
+    return (safediv(1 - 0.5 * pitch * B, g) * data["|B|_r|v,p"] + g * data["K"]) / B
 
 
-def _radial_drift_2(data, B, pitch):
+def _radial_drift_wb_inverse(data, B, pitch):
     return safediv(
         data["cvdrift0"] * (1 - 0.5 * pitch * B),
         jnp.sqrt(jnp.abs(1 - pitch * B)),
     )
 
 
-def _poloidal_drift_secular(data, B, pitch):
+def _poloidal_drift_secular_wb_inverse(data, B, pitch):
     # TODO (#465), multiply by (omega + zeta) instead of zeta
     return safediv(
         (data["gbdrift (periodic)"] + data["gbdrift (secular)/phi"] * data["zeta"])
@@ -104,10 +101,9 @@ def _poloidal_drift_secular(data, B, pitch):
         "b",
         "grad(phi)",
         "grad(psi)",
-        "|grad(psi)|",
         "|grad(rho)|",
         "|e_alpha|r,p|",
-        "kappa_g",
+        "|grad(psi)|*kappa_g",
         "iota_r",
         "V_psi",
     ]
@@ -153,7 +149,7 @@ def Gamma_c(data):
         def fun(pitch_inv):
             points = bounce.points(pitch_inv, opts.num_well)
             v_tau, radial_drift, poloidal_drift = bounce.integrate(
-                [_v_tau, _radial_drift_1, _poloidal_drift_periodic],
+                [_v_tau, _radial_drift, _poloidal_drift_periodic],
                 pitch_inv,
                 data,
                 ["|grad(psi)|*kappa_g", "|B|_r|v,p", "K"],
@@ -182,9 +178,8 @@ def fun(pitch_inv):
         Gamma_c, _gamma_c_data(data), data, angle, grid, opts.surf_batch_size
     )
     assert out.ndim == 1
-    data["Gamma_c"] = (
-        grid.expand(out) / data["V_psi"] / (opts.num_field_periods / grid.NFP * 2**0.5)
-    )
+    scalar = opts.num_field_periods / grid.NFP * 2**0.5
+    data["Gamma_c"] = grid.expand(out / scalar) / data["V_psi"]
     return data
 
 
@@ -208,10 +203,9 @@ def fun(pitch_inv):
         "b",
         "grad(phi)",
         "grad(psi)",
-        "|grad(psi)|",
         "|grad(rho)|",
         "|e_alpha|r,p|",
-        "kappa_g",
+        "|grad(psi)|*kappa_g",
         "iota_r",
     ]
     + Bounce2D.required_names,
@@ -243,7 +237,7 @@ def gamma_c0(data):
         bounce = Bounce2D(grid, data, data["angle"], **opts)
         points = bounce.points(pitch_inv, opts.num_well)
         radial_drift, poloidal_drift = bounce.integrate(
-            [_radial_drift_1, _poloidal_drift_periodic],
+            [_radial_drift, _poloidal_drift_periodic],
             pitch_inv,
             data,
             ["|grad(psi)|*kappa_g", "|B|_r|v,p", "K"],
@@ -324,7 +318,7 @@ def Gamma_c(data):
 
         def fun(pitch_inv):
             v_tau, radial_drift, poloidal_drift = bounce.integrate(
-                [_v_tau, _radial_drift_2, _poloidal_drift_secular],
+                [_v_tau, _radial_drift_wb_inverse, _poloidal_drift_secular_wb_inverse],
                 pitch_inv,
                 data,
                 ["cvdrift0", "gbdrift (periodic)", "gbdrift (secular)/phi"],
@@ -354,7 +348,7 @@ def fun(pitch_inv):
         opts.surf_batch_size,
     )
     assert out.ndim == 1
-    data["Gamma_c Velasco"] = (
-        grid.expand(out) / data["V_psi"] / (opts.num_field_periods / grid.NFP * 2**0.5)
-    )
+
+    scalar = opts.num_field_periods / grid.NFP * 2**0.5
+    data["Gamma_c Velasco"] = grid.expand(out / scalar) / data["V_psi"]
     return data
diff --git a/desc/compute/_metric.py b/desc/compute/_metric.py
index a9ff9685a0..123d39b4ea 100644
--- a/desc/compute/_metric.py
+++ b/desc/compute/_metric.py
@@ -1797,6 +1797,24 @@ def _gradrho(params, transforms, profiles, data, **kwargs):
     return data
 
 
+@register_compute_fun(
+    name="|grad(rho)|*kappa_g",
+    label="|\\nabla\\rho| \\kappa_g",
+    units="m^{-2}",
+    units_long="Inverse square meters",
+    description="Radial coordinate gradient magnitude times geodesic curvature.",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["|grad(rho)|", "kappa_g"],
+)
+def _gradrho_mag_times_kappa_g(params, transforms, profiles, data, **kwargs):
+    data["|grad(rho)|*kappa_g"] = data["|grad(rho)|"] * data["kappa_g"]
+    return data
+
+
 @register_compute_fun(
     name="<|grad(rho)|>",  # same as S(r) / V_r(r)
     label="\\langle \\vert \\nabla \\rho \\vert \\rangle",
@@ -1841,6 +1859,24 @@ def _gradpsi_mag(params, transforms, profiles, data, **kwargs):
     return data
 
 
+@register_compute_fun(
+    name="|grad(psi)|*kappa_g",
+    label="|\\nabla\\psi| \\kappa_g",
+    units="T",
+    units_long="Tesla",
+    description="Toroidal flux gradient magnitude times geodesic curvature.",
+    dim=1,
+    params=[],
+    transforms={},
+    profiles=[],
+    coordinates="rtz",
+    data=["|grad(psi)|", "kappa_g"],
+)
+def _gradpsi_mag_times_kappa_g(params, transforms, profiles, data, **kwargs):
+    data["|grad(psi)|*kappa_g"] = data["|grad(psi)|"] * data["kappa_g"]
+    return data
+
+
 @register_compute_fun(
     name="|grad(psi)|^2",
     label="|\\nabla\\psi|^{2}",
diff --git a/desc/compute/_neoclassical.py b/desc/compute/_neoclassical.py
index 662b5dc514..b9aee852b6 100644
--- a/desc/compute/_neoclassical.py
+++ b/desc/compute/_neoclassical.py
@@ -7,7 +7,7 @@
 from ..batching import batch_map
 from ..integrals.bounce_integral import Bounce2D, Options
 from ..integrals.surface_integral import surface_integrals
-from ..utils import parse_argname_change, safediv
+from ..utils import apply, parse_argname_change, safediv
 from .data_index import register_compute_fun
 
 
@@ -64,9 +64,8 @@ def _dI_2(data, B, pitch):
     data=[
         "min_tz |B|",
         "max_tz |B|",
-        "kappa_g",
+        "|grad(rho)|*kappa_g",
         "R0",
-        "|grad(rho)|",
         "<|grad(rho)|>",
         "V_psi",
     ]
@@ -126,15 +125,15 @@ def fun(pitch_inv):
     )
     out = Bounce2D.batch(
         eps_32,
-        {"|grad(rho)|*kappa_g": data["|grad(rho)|"] * data["kappa_g"]},
+        apply(data, subset=("|grad(rho)|*kappa_g",)),
         data,
         angle,
         grid,
         opts.surf_batch_size,
     )
     assert out.ndim == 1
-    data["effective ripple 3/2"] = scalar * (
-        (B0 / data["<|grad(rho)|>"]) ** 2 * grid.expand(out) / data["V_psi"]
+    data["effective ripple 3/2"] = (
+        (B0 / data["<|grad(rho)|>"]) ** 2 * grid.expand(out * scalar) / data["V_psi"]
     )
     return data
 
diff --git a/desc/compute/_old.py b/desc/compute/_old.py
index 1361f3170d..415b60411e 100644
--- a/desc/compute/_old.py
+++ b/desc/compute/_old.py
@@ -9,12 +9,12 @@
 from desc.backend import jit, jnp
 
 from ..integrals.bounce_integral import Bounce1D, Options
-from ..utils import safediv
+from ..utils import apply, safediv
 from ._fast_ion import (
     _gamma_c_data,
     _poloidal_drift_periodic,
-    _radial_drift_1,
-    _radial_drift_2,
+    _radial_drift,
+    _radial_drift_wb_inverse,
     _v_tau,
 )
 from ._neoclassical import _dI_1, _dI_2
@@ -43,9 +43,8 @@
     data=[
         "min_tz |B|",
         "max_tz |B|",
-        "kappa_g",
         "R0",
-        "|grad(rho)|",
+        "|grad(rho)|*kappa_g",
         "<|grad(rho)|>",
         "fieldline length",
     ]
@@ -85,7 +84,7 @@ def eps_32(data):
     scalar = jnp.pi / (8 * 2**0.5) * data["R0"] ** 2
     out = Bounce1D.batch(
         eps_32,
-        {"|grad(rho)|*kappa_g": data["|grad(rho)|"] * data["kappa_g"]},
+        apply(data, subset=("|grad(rho)|*kappa_g",)),
         data,
         grid,
         opts.surf_batch_size,
@@ -93,8 +92,7 @@ def eps_32(data):
     assert out.ndim == 1
     data["old effective ripple 3/2"] = (
         (B0 / data["<|grad(rho)|>"]) ** 2
-        * scalar
-        * grid.expand(out)
+        * grid.expand(out * scalar)
         / data["fieldline length"]
     )
     return data
@@ -151,10 +149,9 @@ def _effective_ripple_1D(params, transforms, profiles, data, **kwargs):
         "b",
         "grad(phi)",
         "grad(psi)",
-        "|grad(psi)|",
+        "|grad(psi)|*kappa_g",
         "|grad(rho)|",
         "|e_alpha|r,p|",
-        "kappa_g",
         "iota_r",
         "fieldline length",
     ]
@@ -196,7 +193,7 @@ def Gamma_c(data):
         bounce = Bounce1D(grid, data, opts.quad)
         points = bounce.points(pitch_inv, num_well)
         v_tau, radial_drift, poloidal_drift = bounce.integrate(
-            [_v_tau, _radial_drift_1, _poloidal_drift_periodic],
+            [_v_tau, _radial_drift, _poloidal_drift_periodic],
             pitch_inv,
             data,
             ["|grad(psi)|*kappa_g", "|B|_r|v,p", "K"],
@@ -217,7 +214,7 @@ def Gamma_c(data):
     out = Bounce1D.batch(Gamma_c, _gamma_c_data(data), data, grid, opts.surf_batch_size)
     assert out.ndim == 1
     data["old Gamma_c"] = (
-        grid.expand(out) / data["fieldline length"] / (2**1.5 * jnp.pi)
+        grid.expand(out / (2**1.5 * jnp.pi)) / data["fieldline length"]
     )
     return data
 
@@ -259,7 +256,7 @@ def _Gamma_c_Velasco_1D(params, transforms, profiles, data, **kwargs):
     opts = Options.guess(eta=-1, grid=grid, Y_B=grid.num_zeta, **kwargs)
     num_well = kwargs.get("num_well", -1)
 
-    def _poloidal_drift_secular(data, B, pitch):
+    def _poloidal_drift_secular_wb_inverse(data, B, pitch):
         return safediv(
             data["gbdrift"] * (1 - 0.5 * pitch * B),
             jnp.sqrt(jnp.abs(1 - pitch * B)),
@@ -270,7 +267,7 @@ def Gamma_c(data):
             data["min_tz |B|"], data["max_tz |B|"], opts.pitch_quad
         )
         v_tau, radial_drift, poloidal_drift = Bounce1D(grid, data, opts.quad).integrate(
-            [_v_tau, _radial_drift_2, _poloidal_drift_secular],
+            [_v_tau, _radial_drift_wb_inverse, _poloidal_drift_secular_wb_inverse],
             pitch_inv,
             data,
             ["cvdrift0", "gbdrift"],
@@ -291,6 +288,6 @@ def Gamma_c(data):
     )
     assert out.ndim == 1
     data["old Gamma_c Velasco"] = (
-        grid.expand(out) / data["fieldline length"] / (2**1.5 * jnp.pi)
+        grid.expand(out / (2**1.5 * jnp.pi)) / data["fieldline length"]
     )
     return data
diff --git a/desc/compute/_turbulence.py b/desc/compute/_turbulence.py
new file mode 100644
index 0000000000..9d03825005
--- /dev/null
+++ b/desc/compute/_turbulence.py
@@ -0,0 +1,181 @@
+"""Compute functions for turbulent transport."""
+
+from functools import partial
+
+import numpy as np
+from jax.lax import stop_gradient
+from orthax import orthgauss
+from orthax.recurrence import GeneralizedLaguerre
+
+from desc.backend import jit, jnp
+
+from ..batching import batch_map
+from ..integrals.bounce_integral import Bounce2D, Options
+from ..utils import apply, safediv
+from ._fast_ion import _radial_drift
+from .data_index import register_compute_fun
+
+
+def _G_hat_half(data, B, pitch):
+    return safediv(1.0, jnp.sqrt(jnp.abs(1 - pitch * B)))
+
+
+def _binormal_drift_wb_inverse(data, B, pitch):
+    # TODO (#465), multiply by (omega + zeta) instead of zeta
+    gbdrift_secular = data["gbdrift (secular)/phi"] * data["zeta"]
+    cvdrift = data["cvdrift (periodic)"] + gbdrift_secular
+    gbdrift = data["gbdrift (periodic)"] + gbdrift_secular
+    g = jnp.sqrt(jnp.abs(1 - pitch * B))
+    return (cvdrift - 0.5 * gbdrift) * g + safediv(0.5 * gbdrift, g)
+
+
+def _ae(G, G_ω_α, G_ω_ψ, data, energy):
+    shape = (-1,) + (1,) * G.ndim
+
+    G = G[..., None, :]  # This is sqrt G hat.
+    # scale by conjugate widths
+    G_ω_α = G_ω_α[..., None, :] * data["ae psi width"].reshape(shape)
+    G_ω_ψ = G_ω_ψ[..., None, :] * data["ae alpha width"].reshape(shape)
+    η_n = data["ae grad(density)"].reshape(shape)
+    η_T = data["ae grad(temperature)"].reshape(shape)
+    C = η_n - 1.5 * η_T
+    energy = energy[..., None]
+
+    drift = jnp.hypot(G_ω_α, G_ω_ψ)
+    drive = jnp.hypot(G * (η_T + C / energy) - G_ω_α, G_ω_ψ)
+
+    return G_ω_α * C + (G_ω_α * η_T + safediv(drift * (drive - drift), G)) * energy
+
+
+def _energy_quad(num_energy):
+    # The energy integral has weight E^(5/2) exp(-E), but
+    # ω_* = η_T + C / E makes AE(E) ~ C/E for E near zero.
+    return stop_gradient(orthgauss(num_energy, GeneralizedLaguerre(np.array([1.5]))))
+
+
+@register_compute_fun(
+    name="available energy",
+    label="\\widehat{A}",
+    units="~",
+    units_long="None",
+    description="Dimensionless available energy of trapped electrons",
+    dim=1,
+    params=[],
+    transforms={"grid": []},
+    profiles=[],
+    coordinates="r",
+    data=[
+        "min_tz |B|",
+        "max_tz |B|",
+        "psi_r",
+        "rho",
+        "ne",
+        "ne_r",
+        "Te",
+        "Te_r",
+        "cvdrift (periodic)",
+        "gbdrift (periodic)",
+        "gbdrift (secular)/phi",
+        "|grad(psi)|*kappa_g",
+        "V_psi",
+    ]
+    + Bounce2D.required_names,
+    resolution_requirement="tz",
+    grid_requirement={"can_fft2": True},
+    radial_scale="float : Multiplier for the radial correlation length.",
+    binormal_scale="float : Multiplier for the binormal correlation length.",
+    num_energy=(
+        "int : Resolution for generalized Gauss-Laguerre quadrature over energy."
+    ),
+    energy_quad="tuple : Nodes and weights for the energy quadrature.",
+    **Options._doc,
+)
+@partial(
+    jit,
+    static_argnames=Options._static_argnames + ("num_energy",),
+)
+def _available_energy(params, transforms, profiles, data, **kwargs):
+    """Dimensionless available energy of trapped electrons.
+
+    References
+    ----------
+    .. [1] R. J. J. Mackenbach et al., J. Plasma Phys. 89, 905890513 (2023).
+    .. [2] K. Unalmis et al., "Spectrally accurate, reverse-mode
+           differentiable bounce-averaging algorithm and its applications,"
+           Journal of Plasma Physics.
+
+    Parameters
+    ----------
+    radial_scale, binormal_scale : float
+        Correlation-length multipliers. Default is 1.
+    num_energy : int
+        Resolution for generalized Gauss-Laguerre quadrature over energy.
+
+    """
+    # noqa: unused dependency
+    radial_scale = kwargs.get("radial_scale", 1.0)
+    binormal_scale = kwargs.get("binormal_scale", 1.0)
+    energy_quad = kwargs.get("energy_quad", None)
+    if energy_quad is None:
+        energy_quad = _energy_quad(kwargs.get("num_energy", 16))
+
+    grid = transforms["grid"]
+    opts = Options.guess(-1, grid, **kwargs)
+
+    def foreach_surface(data):
+        bounce = Bounce2D(grid, data, data["angle"], **opts)
+
+        def foreach_pitch(pitch_inv):
+            return (
+                _ae(
+                    *bounce.integrate(
+                        [_G_hat_half, _binormal_drift_wb_inverse, _radial_drift],
+                        pitch_inv,
+                        data,
+                        names,
+                        num_well=opts.num_well,
+                        loop=opts.loop,
+                    ),
+                    data,
+                    energy_quad[0],
+                )
+                .sum(-1)
+                .mean(-3)
+                .dot(energy_quad[1])
+            )
+
+        pitch_inv, weight = Bounce2D.pitch_quad(
+            data["min_tz |B|"], data["max_tz |B|"], opts.pitch_quad
+        )
+        return jnp.sum(
+            batch_map(foreach_pitch, pitch_inv, opts.pitch_batch_size)
+            * weight
+            / pitch_inv**2,
+            axis=-1,
+        )
+
+    names = (
+        "cvdrift (periodic)",
+        "gbdrift (periodic)",
+        "gbdrift (secular)/phi",
+        "|grad(psi)|*kappa_g",
+    )
+    out = Bounce2D.batch(
+        foreach_surface,
+        apply(data, subset=names),
+        data,
+        kwargs["angle"],
+        grid,
+        opts.surf_batch_size,
+        surface_data={
+            "ae grad(density)": radial_scale * safediv(data["ne_r"], data["ne"]),
+            "ae psi width": radial_scale * data["psi_r"],
+            "ae alpha width": binormal_scale * safediv(1.0, data["rho"]),
+            "ae grad(temperature)": radial_scale * safediv(data["Te_r"], data["Te"]),
+        },
+    )
+    assert out.ndim == 1
+
+    scalar = jnp.sqrt(jnp.pi) * grid.NFP / (3 * opts.num_field_periods)
+    data["available energy"] = grid.expand(scalar * out) / data["V_psi"]
+    return data
diff --git a/desc/integrals/__init__.py b/desc/integrals/__init__.py
index 6bb01b77f7..51f10b513a 100644
--- a/desc/integrals/__init__.py
+++ b/desc/integrals/__init__.py
@@ -1,5 +1,6 @@
 """Classes for function integration."""
 
+from ._ae_plot import plot_available_energy
 from ._interp_utils import nufft1d2r, nufft2d2r
 from .bounce_integral import Bounce1D, Bounce2D
 from .singularities import (
diff --git a/desc/integrals/_ae_plot.py b/desc/integrals/_ae_plot.py
new file mode 100644
index 0000000000..ba39c279a4
--- /dev/null
+++ b/desc/integrals/_ae_plot.py
@@ -0,0 +1,594 @@
+"""Available-energy plotting functions."""
+
+from dataclasses import dataclass
+
+import matplotlib.colors as colors
+import matplotlib.pyplot as plt
+import numpy as np
+from interpax import PPoly
+from interpax_fft import PiecewiseChebyshevSeries
+from matplotlib.collections import LineCollection
+from matplotlib.ticker import MaxNLocator
+
+from desc.backend import jnp
+from desc.grid import LinearGrid
+from desc.utils import apply, errorif, safediv, setdefault
+
+from .bounce_integral import Bounce2D, Options
+
+
+@dataclass(frozen=True)
+class _AvailableEnergyWellData:
+    """Available-energy data on one flux surface and one field line.
+
+    Attributes
+    ----------
+    rho : float
+        Flux-surface label.
+    alpha : float
+        Field-line label.
+    zeta : ndarray
+        Toroidal field-line coordinate used for plotting ``B``.
+    B : ndarray
+        Magnetic-field strength on ``zeta``.
+    pitch_inv : ndarray
+        Pitch quadrature nodes.
+    points : tuple of ndarray
+        Lower and upper bounce points. Each array has shape
+        ``(num_pitch, num_well)``.
+    valid : ndarray
+        Boolean mask for wells with ordered bounce points, with shape
+        ``(num_pitch, num_well)``.
+    omega_alpha, omega_psi, ae_per_pitch_well : ndarray
+        Per-pitch, per-well bounce quantities with shape
+        ``(num_pitch, num_well)``.
+    """
+
+    rho: float
+    alpha: float
+    zeta: np.ndarray
+    B: np.ndarray
+    pitch_inv: np.ndarray
+    points: tuple[np.ndarray, np.ndarray]
+    valid: np.ndarray
+    omega_alpha: np.ndarray
+    omega_psi: np.ndarray
+    ae_per_pitch_well: np.ndarray
+
+
+def _surface_quantity(grid, template, value):
+    if value is None:
+        return None
+    value = jnp.asarray(value)
+    if value.ndim == 0:
+        return value * jnp.ones_like(template)
+    if value.size == grid.num_rho:
+        return grid.expand(value)
+    return value * jnp.ones_like(template)
+
+
+def _prepare_ae_data(
+    grid,
+    data,
+    angle,
+    bounce_names,
+    radial_scale,
+    binormal_scale,
+    density_gradient,
+    temperature_gradient,
+):
+    fun_data = apply(data, subset=bounce_names)
+    for name in Bounce2D.required_names:
+        fun_data[name] = data[name]
+    fun_data.pop("iota", None)
+
+    for name in fun_data:
+        fun_data[name] = Bounce2D.fourier(Bounce2D.reshape(grid, fun_data[name]))
+
+    density_gradient = _surface_quantity(grid, data["rho"], density_gradient)
+    if density_gradient is None:
+        density_gradient = safediv(data["ne_r"], data["ne"])
+    temperature_gradient = _surface_quantity(grid, data["rho"], temperature_gradient)
+    if temperature_gradient is None:
+        temperature_gradient = safediv(data["Te_r"], data["Te"])
+    errorif(
+        not np.isfinite(density_gradient).all(),
+        msg=("Density gradient was not set for this equilibrium."),
+    )
+    errorif(
+        not np.isfinite(temperature_gradient).all(),
+        msg=("Temperature gradient was not set for this equilibrium."),
+    )
+
+    surface_data = {
+        "ae grad(density)": radial_scale * density_gradient,
+        "ae psi width": radial_scale * data["psi_r"],
+        "ae alpha width": binormal_scale * safediv(1.0, data["rho"]),
+        "ae grad(temperature)": radial_scale * temperature_gradient,
+    }
+    fun_data.update(apply(surface_data, grid.compress))
+    fun_data["iota"] = grid.compress(data["iota"])
+    fun_data["min_tz |B|"] = grid.compress(data["min_tz |B|"])
+    fun_data["max_tz |B|"] = grid.compress(data["max_tz |B|"])
+    fun_data["angle"] = angle
+    return fun_data
+
+
+def _fieldline_B(bounce, num_zeta, rho_idx=0):
+    def evaluate_ppoly(B):
+        P = PPoly(B.T, bounce._c["knots"])
+        return P(zeta)
+
+    if isinstance(bounce._B, PPoly):
+        zeta = np.linspace(bounce._B.x[0], bounce._B.x[-1], num_zeta)
+        return zeta, bounce._B(zeta)[None, :]
+
+    if isinstance(bounce._B, PiecewiseChebyshevSeries):
+        domain = bounce._B.domain
+        B = bounce._B.cheb
+        if B.ndim == 4:
+            B = B[rho_idx]
+        zeta = np.linspace(domain[0], domain[-1], num_zeta)
+        if B.ndim == 2:
+            B = B[None, ...]
+        return zeta, np.stack(
+            [
+                np.ravel(PiecewiseChebyshevSeries(B_alpha, domain).eval1d(zeta))
+                for B_alpha in B
+            ]
+        )
+
+    B = bounce._B
+    if B.ndim == 4:
+        B = B[rho_idx]
+    zeta = np.linspace(bounce._c["knots"][0], bounce._c["knots"][-1], num_zeta)
+    if B.ndim == 2:
+        B = B[None, ...]
+    return zeta, np.stack([evaluate_ppoly(B_alpha) for B_alpha in B])
+
+
+def _ae_well_data(
+    eq,
+    rho=None,
+    alpha=0.0,
+    grid=None,
+    data=None,
+    angle=None,
+    radial_scale=1.0,
+    binormal_scale=1.0,
+    density_gradient=None,
+    temperature_gradient=None,
+    num_energy=16,
+    num_zeta=2000,
+    **kwargs,
+):
+    """Compute per-pitch, per-well available-energy data for plotting.
+
+    Parameters
+    ----------
+    eq : Equilibrium
+        Equilibrium to evaluate.
+    rho : float, optional
+        Flux-surface label. Exactly one value is supported. If omitted with a
+        supplied ``grid``, the single rho value from ``grid`` is used.
+    alpha : float, optional
+        Field-line label. Exactly one value is supported.
+    grid : Grid, optional
+        Single-rho grid used for equilibrium quantities. If omitted, a
+        ``LinearGrid`` is constructed.
+    data : dict, optional
+        Precomputed equilibrium data on ``grid``.
+    angle : ndarray, optional
+        Bounce2D angle map. If omitted, it is computed from ``eq``.
+    radial_scale, binormal_scale : float, optional
+        Correlation-length scale factors.
+    density_gradient, temperature_gradient : float or ndarray, optional
+        Values replacing ``ne_r / ne`` and ``Te_r / Te`` before multiplication
+        by ``radial_scale``. If omitted, the equilibrium profiles are used.
+    num_energy : int, optional
+        Number of generalized Gauss-Laguerre nodes for the energy integral.
+    num_zeta : int, optional
+        Number of points used to plot ``|B|`` along the field line.
+    **kwargs
+        Additional options forwarded to ``Options.guess`` and ``Bounce2D``.
+
+    Returns
+    -------
+    _AvailableEnergyWellData
+        Per-pitch, per-well available-energy data on one flux surface and one
+        field line.
+    """
+    from desc.compute._fast_ion import _radial_drift
+    from desc.compute._turbulence import (
+        _ae,
+        _binormal_drift_wb_inverse,
+        _energy_quad,
+        _G_hat_half,
+    )
+
+    bounce_names = (
+        "cvdrift (periodic)",
+        "gbdrift (periodic)",
+        "gbdrift (secular)/phi",
+        "|grad(psi)|*kappa_g",
+    )
+    surface_names = (
+        "min_tz |B|",
+        "max_tz |B|",
+        "psi_r",
+        "rho",
+        "ne",
+        "ne_r",
+        "Te",
+        "Te_r",
+    )
+
+    rho = None if rho is None else np.asarray(rho, dtype=float).ravel()
+    alpha = np.asarray(alpha, dtype=float).ravel()
+    errorif(
+        rho is not None and rho.size != 1,
+        msg="available-energy well plots require exactly one rho value.",
+    )
+    errorif(
+        alpha.size != 1,
+        msg="available-energy well plots require exactly one alpha value.",
+    )
+
+    X = kwargs.pop("X", 32)
+    Y = kwargs.pop("Y", 32)
+    if grid is None:
+        if rho is None:
+            rho = np.asarray([0.5])
+        grid = LinearGrid(rho=rho, M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP, sym=False)
+    else:
+        errorif(
+            grid.num_rho != 1,
+            msg="available-energy well plots require a single-rho grid.",
+        )
+        grid_rho = grid.compress(grid.nodes[:, 0])
+        if rho is None:
+            rho = grid_rho
+        else:
+            errorif(
+                not np.allclose(rho, grid_rho),
+                msg="rho must match the single rho value in grid.",
+            )
+    if angle is None:
+        angle = Bounce2D.angle(eq, X=X, Y=Y, rho=rho)
+
+    opts = Options.guess(-1, grid, alpha=alpha, **kwargs)
+    compute_names = bounce_names + surface_names + tuple(Bounce2D.required_names)
+    if data is None:
+        data = eq.compute(list(dict.fromkeys(compute_names)), grid=grid)
+
+    fun_data = _prepare_ae_data(
+        grid,
+        data,
+        angle,
+        bounce_names,
+        radial_scale,
+        binormal_scale,
+        density_gradient,
+        temperature_gradient,
+    )
+    bounce = Bounce2D(grid, fun_data, angle, **opts)
+
+    pitch_inv, _ = Bounce2D.pitch_quad(
+        fun_data["min_tz |B|"], fun_data["max_tz |B|"], opts.pitch_quad
+    )
+    points = bounce.points(pitch_inv, opts.num_well)
+    G, G_ω_α, G_ω_ψ = bounce.integrate(
+        [_G_hat_half, _binormal_drift_wb_inverse, _radial_drift],
+        pitch_inv,
+        fun_data,
+        bounce_names,
+        points,
+        loop=opts.loop,
+    )
+
+    energy, energy_weight = _energy_quad(num_energy)
+    ae_per_pitch_well = jnp.sum(
+        _ae(G, G_ω_α, G_ω_ψ, fun_data, energy) * energy_weight[..., None],
+        axis=-2,
+    )
+
+    shape = (-1,) + (1,) * (G.ndim - 1)
+    G_ω_α = G_ω_α * fun_data["ae psi width"].reshape(shape)
+    G_ω_ψ = G_ω_ψ * fun_data["ae alpha width"].reshape(shape)
+    ω_α = safediv(G_ω_α, G)
+    ω_ψ = safediv(G_ω_ψ, G)
+
+    zeta, B = _fieldline_B(bounce, num_zeta)
+    valid = points[0] < points[1]
+    if B.ndim > 1:
+        B = B[0]
+
+    return _AvailableEnergyWellData(
+        rho=rho[0],
+        alpha=alpha[0],
+        zeta=zeta,
+        B=B,
+        pitch_inv=pitch_inv[0],
+        points=(points[0][0, 0], points[1][0, 0]),
+        valid=valid[0, 0],
+        omega_alpha=ω_α[0, 0],
+        omega_psi=ω_ψ[0, 0],
+        ae_per_pitch_well=ae_per_pitch_well[0, 0],
+    )
+
+
+def _well_segments(well_data):
+    zeta1, zeta2 = well_data.points
+    segments = []
+    values = []
+
+    for pitch_idx, y in enumerate(well_data.pitch_inv):
+        for well_idx in range(zeta1.shape[1]):
+            if not well_data.valid[pitch_idx, well_idx]:
+                continue
+            segments.append(
+                [
+                    (zeta1[pitch_idx, well_idx], y),
+                    (zeta2[pitch_idx, well_idx], y),
+                ]
+            )
+            values.append(well_data.ae_per_pitch_well[pitch_idx, well_idx])
+
+    return segments, np.asarray(values)
+
+
+def _color_norm(norm, values, vmin, vmax):
+    if norm is None or norm == "linear":
+        return colors.Normalize(vmin=vmin, vmax=vmax)
+    if norm == "log":
+        positive = values[np.isfinite(values) & (values > 0)]
+        errorif(
+            not positive.size,
+            msg="LogNorm requires at least one positive available-energy value.",
+        )
+        vmin = 1e-5 if vmax > 1e-5 else np.nanmin(positive)
+        return colors.LogNorm(vmin=vmin, vmax=vmax)
+    return norm
+
+
+def _pitch_spacing_linewidths(fig, ax, pitch_inv, segment_pitch, fallback=1.0):
+    pitch_inv = np.asarray(pitch_inv, dtype=float)
+    pitch_inv = np.unique(pitch_inv[np.isfinite(pitch_inv)])
+    if pitch_inv.size < 2:
+        return fallback
+
+    display_y = ax.transData.transform(
+        np.column_stack((np.zeros_like(pitch_inv), pitch_inv))
+    )[:, 1]
+    spacing = np.diff(np.sort(display_y))
+    if not np.any(spacing > 0):
+        return fallback
+
+    spacing = np.concatenate(
+        ([spacing[0]], np.minimum(spacing[:-1], spacing[1:]), [spacing[-1]])
+    )
+    # Adjacent strokes that exactly touch can still leave subpixel white seams
+    # after antialiasing. Overscan slightly in display space to make filled
+    # pitch bands visually continuous.
+    LINEWIDTH_OVERSCAN_PIXELS = 1.0
+    max_linewidth = (spacing + LINEWIDTH_OVERSCAN_PIXELS) * 72 / fig.dpi
+    idx = np.searchsorted(pitch_inv, segment_pitch)
+    idx = np.clip(idx, 0, pitch_inv.size - 1)
+    prev_idx = np.clip(idx - 1, 0, pitch_inv.size - 1)
+    idx = np.where(
+        np.abs(pitch_inv[prev_idx] - segment_pitch)
+        < np.abs(pitch_inv[idx] - segment_pitch),
+        prev_idx,
+        idx,
+    )
+    return max_linewidth[idx]
+
+
+def _add_well_segments(
+    fig,
+    ax,
+    segments,
+    values,
+    pitch_inv,
+    cmap,
+    normalize_ae,
+    norm,
+    linewidth,
+):
+    if not values.size:
+        return None
+    auto_linewidth = linewidth is None or not np.isfinite(linewidth)
+    initial_linewidth = 1.0 if auto_linewidth else linewidth
+
+    max_value = np.nanmax(values)
+    if normalize_ae and max_value > 0:
+        values = values / max_value
+        cbar_label = r"$\widehat{A}_\lambda / \widehat{A}_{\lambda,\max}$"
+        vmin, vmax = 0.0, 1.0
+    else:
+        cbar_label = r"$\widehat{A}_\lambda$"
+        vmin, vmax = 0.0, max_value
+
+    collection = LineCollection(
+        segments,
+        array=values,
+        cmap=cmap,
+        norm=_color_norm(norm, values, vmin, vmax),
+        linewidths=initial_linewidth,
+        alpha=0.95,
+        zorder=1,
+    )
+    ax.add_collection(collection)
+    cbar = fig.colorbar(collection, ax=ax, pad=0.08, fraction=0.055)
+    cbar.set_label(cbar_label)
+    if not isinstance(collection.norm, colors.LogNorm):
+        cbar.ax.yaxis.set_major_locator(MaxNLocator(5))
+    if auto_linewidth:
+        fig.canvas.draw()
+        segment_pitch = np.asarray(segments, dtype=float)[:, 0, 1]
+        collection.set_linewidths(
+            _pitch_spacing_linewidths(fig, ax, pitch_inv, segment_pitch)
+        )
+    return collection
+
+
+def _drift_samples(well_data):
+    zeta1, zeta2 = well_data.points
+    valid = well_data.valid
+    if not np.any(valid):
+        empty = np.empty(0)
+        return empty, empty, empty
+
+    roots = np.column_stack((zeta1[valid], zeta2[valid])).ravel()
+    omega_alpha = np.repeat(well_data.omega_alpha[valid], 2)
+    omega_psi = np.repeat(well_data.omega_psi[valid], 2)
+    return roots, omega_alpha, omega_psi
+
+
+def _add_drifts(ax, well_data, omega_psi_color, omega_alpha_color):
+    drift_ax = ax.twinx()
+    roots, omega_alpha, omega_psi = _drift_samples(well_data)
+    if roots.size:
+        drift_ax.scatter(
+            roots,
+            omega_psi,
+            color=omega_psi_color,
+            marker=".",
+            s=9,
+            alpha=0.9,
+            label=r"$\widehat{\omega}_\psi$",
+        )
+        drift_ax.scatter(
+            roots,
+            omega_alpha,
+            color=omega_alpha_color,
+            marker=".",
+            s=9,
+            alpha=0.9,
+            label=r"$\widehat{\omega}_\alpha$",
+        )
+        drift_ax.axhline(0.0, color="0.15", linestyle="dotted", linewidth=1.0)
+        legend = drift_ax.legend(
+            loc="lower right",
+            facecolor="white",
+            edgecolor="0.25",
+            framealpha=0.95,
+            markerscale=2.5,
+            scatterpoints=3,
+        )
+        legend.get_frame().set_linewidth(0.8)
+    drift_ax.set_ylabel(r"$\widehat{\omega}_\alpha,\ \widehat{\omega}_\psi$")
+    drift_ax.yaxis.set_major_locator(MaxNLocator(5))
+    return drift_ax
+
+
+def _style_well_axes(ax, zeta, B):
+    zeta_min = np.nanmin(zeta)
+    zeta_max = np.nanmax(zeta)
+    y_min = np.nanmin(B)
+    y_max = np.nanmax(B)
+    y_range = y_max - y_min
+    y_pad = 0.04 * y_range if y_range > 0 else max(0.04 * abs(y_max), 1.0)
+
+    ax.set_xlim(zeta_min, zeta_max)
+    ax.set_ylim(y_min - y_pad, y_max + y_pad)
+    ax.set_ylabel(r"$|B|$")
+    ax.set_xlabel(r"$\zeta$")
+    ax.xaxis.set_major_locator(MaxNLocator(6))
+    ax.yaxis.set_major_locator(MaxNLocator(5))
+    ax.grid(color="0.88", linewidth=0.6)
+    ax.set_axisbelow(True)
+
+
+def plot_available_energy(
+    eq,
+    ax=None,
+    cmap=None,
+    normalize_ae=True,
+    norm="log",
+    include_drifts=True,
+    linewidth=None,
+    omega_psi_color=None,
+    omega_alpha_color=None,
+    return_data=False,
+    **kwargs,
+):
+    """Plot an AEpy-style available-energy map over bounce wells.
+
+    Parameters
+    ----------
+    eq : Equilibrium
+        Equilibrium to evaluate.
+    ax : matplotlib.axes.Axes, optional
+        Axes to draw into. If omitted, a new figure and axes are created.
+    cmap : str or Colormap, optional
+        Colormap for the available-energy well segments. If omitted, ``"turbo"``
+        is used for log normalization and ``"plasma"`` is used otherwise.
+    normalize_ae : bool, optional
+        If True, normalize segment colors by the maximum plotted
+        available-energy value.
+    norm : {"log", "linear"} or matplotlib.colors.Normalize, optional
+        Color normalization for available-energy well segments. ``"log"``
+        uses ``matplotlib.colors.LogNorm`` with cutoff at ``1e-5``.
+    include_drifts : bool, optional
+        If True, overlay drift-frequency samples on a secondary y-axis.
+    linewidth : float or None, optional
+        Width of the available-energy well segments. If ``None`` or ``np.inf``,
+        fill the pitch-grid spacing.
+    omega_psi_color, omega_alpha_color : color, optional
+        Colors for the overlaid drift-frequency markers.
+    return_data : bool, optional
+        If True, return ``well_data`` along with the figure and axes.
+
+    Returns
+    -------
+    fig : matplotlib.figure.Figure
+        Figure containing the plot.
+    ax : matplotlib.axes.Axes
+        Primary axes containing ``|B|`` and the well segments.
+    well_data : _AvailableEnergyWellData, optional
+        Returned only when ``return_data`` is True.
+
+    """
+    well_data = _ae_well_data(eq, **kwargs)
+
+    fig, ax = (
+        plt.subplots(figsize=kwargs.pop("figsize", (7, 4.5)), constrained_layout=True)
+        if ax is None
+        else (ax.figure, ax)
+    )
+
+    zeta = well_data.zeta
+    ax.plot(zeta, well_data.B, color="black", linewidth=2, zorder=3)
+    _style_well_axes(ax, zeta, well_data.B)
+
+    segments, values = _well_segments(well_data)
+    _add_well_segments(
+        fig,
+        ax,
+        segments,
+        values,
+        well_data.pitch_inv,
+        setdefault(
+            cmap,
+            (
+                "turbo"
+                if (norm == "log" or isinstance(norm, colors.LogNorm))
+                else "plasma"
+            ),
+        ),
+        normalize_ae,
+        norm,
+        linewidth,
+    )
+
+    is_log = norm == "log" or isinstance(norm, colors.LogNorm)
+    omega_psi_color = setdefault(
+        omega_psi_color, "tab:purple" if is_log else "springgreen"
+    )
+    omega_alpha_color = setdefault(omega_alpha_color, "tab:gray")
+
+    if include_drifts:
+        _add_drifts(ax, well_data, omega_psi_color, omega_alpha_color)
+
+    return (fig, ax, well_data) if return_data else (fig, ax)
diff --git a/desc/integrals/bounce_integral.py b/desc/integrals/bounce_integral.py
index 4c8fcbcf52..34fb16b292 100644
--- a/desc/integrals/bounce_integral.py
+++ b/desc/integrals/bounce_integral.py
@@ -2,7 +2,7 @@
 
 import warnings
 from abc import ABC, abstractmethod
-from typing import NamedTuple
+from typing import NamedTuple, Union
 
 import equinox as eqx
 from interpax import CubicHermiteSpline, PPoly
@@ -270,7 +270,7 @@ class Bounce2D(_Bounce):
     _modes_t: jax.Array
     _c: dict[str, jax.Array]
     _theta: PiecewiseChebyshevSeries
-    _B: jax.Array or PiecewiseChebyshevSeries
+    _B: Union[jax.Array, PiecewiseChebyshevSeries]
     _nufft_eps: float = eqx.field(static=True)
 
     def __init__(
@@ -363,8 +363,18 @@ def __init__(
                 self._modes_z,
             )
 
+    # TODO(#2152): clean this up when angle is removed.
     @staticmethod
-    def batch(fun, fun_data, desc_data, angle, grid, surf_batch_size=1, sparse=True):
+    def batch(
+        fun,
+        fun_data,
+        desc_data,
+        angle,
+        grid,
+        surf_batch_size=1,
+        sparse=True,
+        surface_data=None,
+    ):
         """Compute function ``fun`` over phase space in batches.
 
         This is a utility method to compute some function of bounce integrals
@@ -405,6 +415,9 @@ def batch(fun, fun_data, desc_data, angle, grid, surf_batch_size=1, sparse=True)
             the final objective of interest is a lower dimensional quantity
             than the output, it may be preferable to delay the vjp
             by setting to ``False``.
+        surface_data : dict[str, jnp.ndarray]
+            Data constant on each flux surface. These are compressed over ``rho``
+            and added to ``fun_data`` instead of reshaped and Fourier transformed.
 
         Returns
         -------
@@ -414,8 +427,12 @@ def batch(fun, fun_data, desc_data, angle, grid, surf_batch_size=1, sparse=True)
         for name in Bounce2D.required_names:
             fun_data[name] = desc_data[name]
         fun_data.pop("iota", None)
+
         for name in fun_data:
             fun_data[name] = Bounce2D.fourier(Bounce2D.reshape(grid, fun_data[name]))
+        if surface_data is not None:
+            fun_data.update(apply(surface_data, grid.compress))
+
         fun_data["iota"] = grid.compress(desc_data["iota"])
         fun_data["min_tz |B|"] = grid.compress(desc_data["min_tz |B|"])
         fun_data["max_tz |B|"] = grid.compress(desc_data["max_tz |B|"])
diff --git a/desc/objectives/__init__.py b/desc/objectives/__init__.py
index 52ca059e98..3d05fa66a8 100644
--- a/desc/objectives/__init__.py
+++ b/desc/objectives/__init__.py
@@ -57,6 +57,7 @@
 from ._power_balance import FusionPower, HeatingPowerISS04
 from ._profiles import Pressure, RotationalTransform, Shear, ToroidalCurrent
 from ._stability import BallooningStability, MagneticWell, MercierStability
+from ._turbulence import AvailableEnergy
 from .getters import (
     get_equilibrium_objective,
     get_fixed_axis_constraints,
diff --git a/desc/objectives/_turbulence.py b/desc/objectives/_turbulence.py
new file mode 100644
index 0000000000..b64db217bf
--- /dev/null
+++ b/desc/objectives/_turbulence.py
@@ -0,0 +1,195 @@
+"""Objectives for turbulence proxies."""
+
+from desc.backend import jnp
+from desc.compute._turbulence import _energy_quad
+from desc.compute.utils import _compute as compute_fun
+from desc.integrals.bounce_integral import Options
+from desc.utils import warnif
+
+from .objective_funs import _Objective, collect_docs, doc_bounce
+from .utils import errorif
+
+
+class AvailableEnergy(_Objective):
+    """Available energy of trapped electrons.
+
+    The available-energy metric estimates the dimensionless free energy available
+    to trapped electrons from density and temperature profile gradients.
+
+    References
+    ----------
+    [1] Mackenbach et al., J. Plasma Phys. 89, 905890513 (2023).
+
+    [2] Spectrally accurate, reverse-mode differentiable bounce-averaging algorithm
+        and its applications. Kaya Unalmis et al. Journal of Plasma Physics.
+
+    """
+
+    __doc__ = (
+        __doc__.rstrip()
+        + doc_bounce
+        + collect_docs(
+            target_default="``target=0``.",
+            bounds_default="``target=0``.",
+            normalize_detail=" Note: Has no effect for this objective.",
+            normalize_target_detail=" Note: Has no effect for this objective.",
+            jac_chunk_size=False,
+        )
+    )
+
+    _static_attrs = _Objective._static_attrs + ["_hyperparam"]
+
+    _coordinates = "r"
+    _units = "~"
+    _print_value_fmt = "Available energy: "
+
+    def __init__(
+        self,
+        eq,
+        *,
+        target=None,
+        bounds=None,
+        weight=1,
+        normalize=True,
+        normalize_target=True,
+        loss_function=None,
+        deriv_mode="auto",
+        name="Available energy",
+        grid=None,
+        X=32,
+        Y=32,
+        Y_B=None,
+        alpha=None,
+        num_field_periods=20,
+        num_well=None,
+        num_quad=32,
+        num_pitch=65,
+        num_energy=16,
+        radial_scale=1.0,
+        binormal_scale=1.0,
+        pitch_batch_size=None,
+        surf_batch_size=1,
+        nufft_eps=1e-7,
+        spline=True,
+    ):
+        errorif(
+            deriv_mode == "fwd",
+            ValueError,
+            "Reverse mode should be used for the objective: AvailableEnergy.",
+        )
+        try:
+            import jax_finufft  # noqa: F401
+        except Exception:
+            warnif(
+                nufft_eps >= 1e-14,
+                msg="\njax-finufft is not installed properly.\n"
+                "Setting parameter nufft_eps to zero.\n"
+                "Performance may be somewhat slower.\n",
+            )
+            nufft_eps = 0.0
+        nufft_eps = float(nufft_eps)
+
+        if target is None and bounds is None:
+            target = 0.0
+
+        self._grid = grid
+        if alpha is None:
+            alpha = jnp.zeros(1)
+        self._constants = {"quad_weights": 1.0, "alpha": alpha}
+        self._hyperparam = {
+            "X": X,
+            "Y": Y,
+            "Y_B": Y_B,
+            "num_field_periods": num_field_periods,
+            "num_well": num_well,
+            "num_quad": num_quad,
+            "num_pitch": num_pitch,
+            "num_energy": num_energy,
+            "radial_scale": radial_scale,
+            "binormal_scale": binormal_scale,
+            "pitch_batch_size": pitch_batch_size,
+            "surf_batch_size": surf_batch_size,
+            "nufft_eps": nufft_eps,
+            "spline": spline,
+        }
+
+        super().__init__(
+            things=eq,
+            target=target,
+            bounds=bounds,
+            weight=weight,
+            normalize=normalize,
+            normalize_target=normalize_target,
+            loss_function=loss_function,
+            deriv_mode=deriv_mode,
+            name=name,
+            jac_chunk_size=None,
+        )
+
+    def build(self, use_jit=True, verbose=1):
+        """Build constant arrays.
+
+        Parameters
+        ----------
+        use_jit : bool, optional
+            Whether to just-in-time compile the objective and derivatives.
+        verbose : int, optional
+            Level of output.
+
+        """
+        Options._build_objective(self, "available energy", eta=-1)
+        self._constants["energy_quad"] = _energy_quad(
+            self._hyperparam.pop("num_energy")
+        )
+        super().build(use_jit=use_jit, verbose=verbose)
+
+    def compute(self, params, constants=None):
+        """Compute the available energy.
+
+        Parameters
+        ----------
+        params : dict
+            Dictionary of equilibrium degrees of freedom, e.g.
+            ``Equilibrium.params_dict``.
+        constants : dict
+            Dictionary of constant data, e.g. transforms, profiles etc.
+            Defaults to ``self.constants``.
+
+        Returns
+        -------
+        available_energy : ndarray
+            Available energy as a function of the flux surface label.
+
+        """
+        constants = self._get_deprecated_constants(constants)
+        eq = self.things[0]
+
+        data = compute_fun(
+            eq, "iota", params, constants["transforms"], constants["profiles"]
+        )
+        delta = eq._map_poloidal_coordinates(
+            constants["transforms"]["grid"].compress(data["iota"]),
+            constants["x"],
+            constants["y"],
+            params["L_lmn"],
+            constants["lambda"],
+            outbasis="delta",
+            # TODO (#1034): Use old theta values as initial guess.
+            tol=1e-8,
+        )[..., ::-1]
+
+        data = compute_fun(
+            eq,
+            "available energy",
+            params,
+            constants["transforms"],
+            constants["profiles"],
+            data,
+            angle=delta,
+            alpha=constants["alpha"],
+            quad=constants["quad"],
+            energy_quad=constants.get("energy_quad", None),
+            _vander=constants["_vander"],
+            **self._hyperparam,
+        )
+        return constants["transforms"]["grid"].compress(data["available energy"])
diff --git a/docs/api_objectives.rst b/docs/api_objectives.rst
index 72b6468e78..a16206a4d5 100644
--- a/docs/api_objectives.rst
+++ b/docs/api_objectives.rst
@@ -73,6 +73,16 @@ Neoclassical
     desc.objectives.EffectiveRipple
 
 
+Turbulence
+----------
+.. autosummary::
+    :toctree: _api/objectives
+    :recursive:
+    :template: class.rst
+
+    desc.objectives.AvailableEnergy
+
+
 Omnigenity
 ----------
 .. autosummary::
diff --git a/publications/unalmis2025/spectral_reverse_bounce_supplement.pdf b/publications/unalmis2025/spectral_reverse_bounce_supplement.pdf
index c37eed3250..c2a19dca20 100644
Binary files a/publications/unalmis2025/spectral_reverse_bounce_supplement.pdf and b/publications/unalmis2025/spectral_reverse_bounce_supplement.pdf differ
diff --git a/requirements.txt b/requirements.txt
index c0621a9f43..deacd1c7fe 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,7 +12,7 @@ netcdf4 >= 1.5.4, !=1.7.4, <= 1.7.5
 numpy >= 1.20.0, <= 2.5
 nvidia-ml-py >=12.535.77
 optax < 0.3.0
-orthax < 0.3
+orthax >= 0.2.7, < 0.3
 plotly >= 5.16, <= 6.7.0
 psutil <= 7.2.2
 pylatexenc >= 2.0, <= 2.10
diff --git a/tests/baseline/test_plot_available_energy.png b/tests/baseline/test_plot_available_energy.png
new file mode 100644
index 0000000000..6557dc379f
Binary files /dev/null and b/tests/baseline/test_plot_available_energy.png differ
diff --git a/tests/inputs/master_compute_data_rpz.pkl b/tests/inputs/master_compute_data_rpz.pkl
index 608a3407d5..f69a0e96f9 100644
Binary files a/tests/inputs/master_compute_data_rpz.pkl and b/tests/inputs/master_compute_data_rpz.pkl differ
diff --git a/tests/test_compute_everything.py b/tests/test_compute_everything.py
index 7a7e021e16..7573b69787 100644
--- a/tests/test_compute_everything.py
+++ b/tests/test_compute_everything.py
@@ -311,7 +311,7 @@ def fft_grid_data(p):
         return {}
 
     # TODO: can automate this later to add omnigeneity, boozer transform, etc.
-    fft_names = ["effective ripple", "Gamma_c", "Gamma_c Velasco"]
+    fft_names = ["effective ripple", "Gamma_c", "Gamma_c Velasco", "available energy"]
 
     eq = get("W7-X")
     rho = np.linspace(0, 1, 10)
diff --git a/tests/test_integrals.py b/tests/test_integrals.py
index 52828abf16..859f1bd462 100644
--- a/tests/test_integrals.py
+++ b/tests/test_integrals.py
@@ -15,6 +15,7 @@
 
 from desc.backend import jnp, vmap
 from desc.basis import FourierZernikeBasis
+from desc.compute._turbulence import _binormal_drift_wb_inverse
 from desc.equilibrium import Equilibrium
 from desc.equilibrium.coords import get_rtz_grid
 from desc.examples import get
@@ -1614,18 +1615,6 @@ def _not_part_of_tutorial_test(
         l, m = 1, 0
         _, _ = bounce.plot(l, m, pitch_inv[l], show=False)
 
-    @staticmethod
-    def drift_num_integrand(data, B, pitch):
-        """Integrand of numerator of bounce averaged binormal drift."""
-        g = jnp.sqrt(jnp.abs(1 - pitch * B))
-        cvdrift = (
-            data["cvdrift (periodic)"] + data["gbdrift (secular)/phi"] * data["zeta"]
-        )
-        gbdrift = (
-            data["gbdrift (periodic)"] + data["gbdrift (secular)/phi"] * data["zeta"]
-        )
-        return (cvdrift * g) - (0.5 * g * gbdrift) + (0.5 * gbdrift / g)
-
     @pytest.mark.unit
     @pytest.mark.mpl_image_compare(remove_text=True, tolerance=tol_1d)
     @pytest.mark.parametrize(
@@ -1664,7 +1653,7 @@ def test_binormal_drift_bounce2d(self, nufft_eps, spline, Y_B):
 
         data = {name: Bounce2D.reshape(grid, grid_data[name]) for name in names}
         drift_numerical_num, drift_numerical_den = bounce.integrate(
-            [TestBounce2D.drift_num_integrand, TestBounce.drift_den_integrand],
+            [_binormal_drift_wb_inverse, TestBounce.drift_den_integrand],
             pitch_inv,
             data,
             points=points,
@@ -1681,7 +1670,7 @@ def test_binormal_drift_bounce2d(self, nufft_eps, spline, Y_B):
         )
 
         TestBounce._test_bounce_autodiff(
-            bounce, TestBounce2D.drift_num_integrand, data, nufft_eps
+            bounce, _binormal_drift_wb_inverse, data, nufft_eps
         )
 
         fig, ax = plt.subplots()
diff --git a/tests/test_objective_funs.py b/tests/test_objective_funs.py
index 18f02612a5..341bbcbbad 100644
--- a/tests/test_objective_funs.py
+++ b/tests/test_objective_funs.py
@@ -43,6 +43,7 @@
 )
 from desc.objectives import (
     AspectRatio,
+    AvailableEnergy,
     BallooningStability,
     BootstrapRedlConsistency,
     BoundaryError,
@@ -2083,6 +2084,11 @@ def test(field, grid, regularization):
     def test_objective_against_compute_bounce(self):
         """Test objectives are built properly."""
         eq = get("W7-X")
+        eq.pressure = None
+        eq.electron_density = PowerSeriesProfile([1e19, 0, -5e18])
+        eq.electron_temperature = PowerSeriesProfile([1e3, 0, -5e2])
+        eq.ion_temperature = PowerSeriesProfile([1e3, 0, -5e2])
+        eq.atomic_number = 1.0
         rho = np.linspace(0.1, 1, 3)
         obj_grid = LinearGrid(rho=rho, M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP, sym=False)
         X = 16
@@ -2112,6 +2118,15 @@ def test_objective_against_compute_bounce(self):
         np.testing.assert_allclose(
             obj.compute(eq.params_dict), grid.compress(data[names[1]])
         )
+        data = eq.compute("available energy", grid, angle=angle, num_energy=8, **opts)
+        obj = AvailableEnergy(
+            eq, grid=obj_grid, nufft_eps=1e-7, X=X, Y=Y, num_energy=8, **opts
+        )
+        obj.build()
+        assert obj._hyperparam["num_well"] == opts["num_well"]
+        np.testing.assert_allclose(
+            obj.compute(eq.params_dict), grid.compress(data["available energy"])
+        )
 
     @pytest.mark.unit
     def test_objective_against_compute_ballooning(self):
@@ -3249,13 +3264,16 @@ def test_loss_function_asserts():
 
 def _reduced_resolution_objective(eq, objective, **kwargs):
     """Speed up testing suite by defining rules to reduce objective resolution."""
-    if objective in {EffectiveRipple, GammaC}:
+    if objective in {AvailableEnergy, EffectiveRipple, GammaC}:
         kwargs["X"] = 16
         kwargs["Y"] = 24
         kwargs["num_field_periods"] = 10
         kwargs["num_well"] = 15 * kwargs["num_field_periods"] // eq.NFP
         kwargs["num_pitch"] = 24
         kwargs["num_quad"] = 16
+    if objective is AvailableEnergy:
+        kwargs["num_energy"] = 8
+        kwargs.setdefault("Y_B", 16)
     return objective(eq=eq, **kwargs)
 
 
@@ -3274,6 +3292,7 @@ class TestComputeScalarResolution:
     ]
     specials = [
         # these require special logic
+        AvailableEnergy,
         BootstrapRedlConsistency,
         BoundaryError,
         CoilArclengthVariance,
@@ -3311,6 +3330,32 @@ class TestComputeScalarResolution:
     eq = get("HELIOTRON")
     res_array = np.array([2, 2.5, 3])
 
+    @pytest.mark.regression
+    def test_compute_scalar_resolution_available_energy(self):
+        """AvailableEnergy."""
+        eq = self.eq.copy()
+        eq.pressure = None
+        eq.electron_density = PowerSeriesProfile([1e19, 0, -5e18])
+        eq.electron_temperature = PowerSeriesProfile([1e3, 0, -5e2])
+        eq.ion_temperature = PowerSeriesProfile([1e3, 0, -5e2])
+        eq.atomic_number = 1.0
+        f = np.zeros_like(self.res_array, dtype=float)
+        for i, res in enumerate(self.res_array):
+            grid = LinearGrid(
+                M=int(self.eq.M * res),
+                N=int(self.eq.N * res),
+                NFP=self.eq.NFP,
+                rho=np.array([0.3, 0.7]),
+                sym=False,
+            )
+            obj = ObjectiveFunction(
+                _reduced_resolution_objective(eq, AvailableEnergy, grid=grid),
+                use_jit=False,
+            )
+            obj.build(verbose=0)
+            f[i] = obj.compute_scalar(obj.x())
+        np.testing.assert_allclose(f, f[-1], rtol=1e-1)
+
     @pytest.mark.regression
     def test_compute_scalar_resolution_plasma_vessel(self):
         """PlasmaVesselDistance."""
@@ -3792,6 +3837,7 @@ class TestObjectiveNaNGrad:
     ]
     specials = [
         # these require special logic
+        AvailableEnergy,
         BallooningStability,
         BootstrapRedlConsistency,
         BoundaryError,
@@ -3827,6 +3873,33 @@ class TestObjectiveNaNGrad:
     ]
     other_objectives = list(set(objectives) - set(specials))
 
+    @pytest.mark.unit
+    def test_objective_no_nangrad_available_energy(self):
+        """AvailableEnergy."""
+        eq = get("ESTELL")
+        with pytest.warns(UserWarning, match="Reducing radial"):
+            eq.change_resolution(2, 2, 2, 4, 4, 4)
+        eq.pressure = None
+        eq.electron_density = PowerSeriesProfile([1e19, 0, -5e18])
+        eq.electron_temperature = PowerSeriesProfile([1e3, 0, -5e2])
+        eq.ion_temperature = PowerSeriesProfile([1e3, 0, -5e2])
+        eq.atomic_number = 1.0
+
+        obj_0 = ObjectiveFunction(
+            _reduced_resolution_objective(eq, AvailableEnergy, nufft_eps=0)
+        )
+        obj_0.build(verbose=0)
+        g_0 = obj_0.grad(obj_0.x())
+        assert not np.any(np.isnan(g_0))
+
+        obj = ObjectiveFunction(
+            _reduced_resolution_objective(eq, AvailableEnergy, nufft_eps=1e-8)
+        )
+        obj.build(verbose=0)
+        g = obj.grad(obj.x())
+        assert not np.any(np.isnan(g))
+        np.testing.assert_allclose(g, g_0, atol=5e-5)
+
     @pytest.mark.unit
     def test_objective_no_nangrad_plasma_vessel(self):
         """PlasmaVesselDistance."""
diff --git a/tests/test_plotting.py b/tests/test_plotting.py
index 525111ba32..10f8453a3b 100644
--- a/tests/test_plotting.py
+++ b/tests/test_plotting.py
@@ -21,6 +21,7 @@
 from desc.geometry import FourierRZToroidalSurface, FourierXYZCurve
 from desc.grid import ConcentricGrid, Grid, LinearGrid, QuadratureGrid
 from desc.integrals import surface_averages
+from desc.integrals._ae_plot import plot_available_energy
 from desc.io import load
 from desc.magnetic_fields import (
     OmnigenousField,
@@ -1063,6 +1064,23 @@ def test_plot_b_mag():
     return fig
 
 
+@pytest.mark.unit
+@pytest.mark.mpl_image_compare(remove_text=True, tolerance=tol_2d)
+def test_plot_available_energy():
+    """Test plotting available-energy wells for NCSX."""
+    eq = get("NCSX")
+    grid = LinearGrid(rho=[0.5], M=eq.M_grid, N=eq.N_grid, NFP=eq.NFP, sym=False)
+    fig, ax = plot_available_energy(
+        eq=eq,
+        grid=grid,
+        density_gradient=-3.0,
+        temperature_gradient=0.0,
+        num_quad=16,
+        num_field_periods=10,
+    )
+    return fig
+
+
 @pytest.mark.unit
 @pytest.mark.mpl_image_compare(remove_text=True, tolerance=tol_2d)
 def test_plot_coefficients():