NonlinearSolve uses AD Jacobian even when I specify it analytically

Specific Domains Numerics
3

So here is an MWE that demonstrates unexpected behaviour (for me) at multple places:

using NonlinearSolve

function r(gamma, params_nonlinear)
    println("r") # Does show up
    return gamma^params_nonlinear.exponent
end

function dr(gamma, params_nonlinear)
    println("dr") # Does not show up
    return params_nonlinear.exponent * gamma^(params_nonlinear.exponent - 1.0)
end

h = NonlinearFunction(r; jac = dr)
gamma_0 = 1.0
params_nonlinear = (exponent = 2.0, )
sol = NonlinearSolve.solve(NonlinearProblem(h, gamma_0, params_nonlinear), 
                           NewtonRaphson(;concrete_jac = true, autodiff = false), 
                           abstol = 2 * eps(Float64))
println(sol.stats) # SciMLBase.NLStats(27, 26, 0, 26, 26) => 26 Jacobians created