Function rkmk_dae_step_first_order_root 

pub fn rkmk_dae_step_first_order_root<Field, Tab, F, J, Z, T>(
    rate: &mut impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<Derivative<Field::Increment, T>, String>,
    function: impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<F, String>,
    jacobian: impl FnMut(Quantity<T>, &Field::Point, &Z) -> Result<J, String>,
    solver: &impl FirstOrderRootFinding<F, J, Z>,
    point: &Field::Point,
    z: &Z,
    t: Quantity<T>,
    dt: Quantity<T>,
    scratch: &mut Vec<Field::Increment>,
    first_rate: Option<&Derivative<Field::Increment, T>>,
    equality_constraint: impl FnMut(Quantity<T>) -> EqualityConstraint,
) -> Result<(Field::Point, Z, Option<Derivative<Field::Increment, T>>), IntegrationError>
where
    Field: Integrable,
    Tab: ButcherTableau,
    Field::Point: Clone,
    Field::Increment: Clone + Differentiable<T>,
    Z: Clone,
    T: Copy,
    Quantity<T>: Mul<Scalar, Output = Quantity<T>>,
    for<'a> &'a Derivative<Field::Increment, T>: Mul<Quantity<T>, Output = Field::Increment>,

rkmk_dae_step with the algebraic unknown resolved by first-order root-finding at every stage abscissa, built from function/jacobian/ solver exactly as ExplicitDaeVariableStepExplicitFirstOrderRoot builds its solution closure for the legacy flat DAE solver — the split between root-finding and minimization is orthogonal to which field the state lives on, so this is the one place that wrapping happens for the RKMK-DAE path. Any super::StateEvolution model that also supplies a residual and its Jacobian in terms of the whole field state gets the manifold-aware stage-equilibrium step for free, without hand-rolling this closure itself.