You can get closer to second order if your error is a sum of squares of the form 0.5 * f^T * f, in which case you case compute the jacobian J of the fitness f, and approximate your hessian as J^T * J. That's the Gauss-Newton method and it's very effective if a linear approximation of f holds locally.

That's an approximation but that's arguably a second order method, since it acknowledges the nonlinearity coming from the squared norm of f.