How to get gradients from scipy.optimize.least_squares for each iteration?

Im working with a package that uses scipy.optimize.least_squares and offers the possibility to use a callback function. However, only objective function parameters are passed to the callback function (for each function evaluation). I need to get the gradient of the function with respect to the parameters, though. I know that scipy.optimize.least_squares computes the gradients for each evaluation via J.T.dot(f), where J is the Jacobian and f is the function. But the gradient is only returned, after the optimization is finished, via the results and shows the gradient for the last iteration / evaluation (solution). Is there a way to get the gradient of the cost function with respect to the parameters for each evaluation?