InverseJacobian & BroydenFirst as preconditioner: what is it exactly going on under the hoods? [scipy]

Question

I am struggling to precisely understand how scipy.optimize.InverseJacobian & scipy.optimize.BroydenFirst operate under the hood. In particular, when used as inner_M in scipy.optimize.newton_krylov, what happens?

For instance, if GMRES is selected, then the class should solve the equivalent Newton problem of the kind: J(x^k)delta_x^k=-F(x^k), but by approximating the matrix vector product J(x^k)delta_x^k, instead of directly computing it. I have read that having a preconditioner that is (an approximation of) the inverse Jacobian, helps in converging, so I am assuming that scipy.optimize.InverseJacobian computes some sort of approximation of J(x^k), but I was wondering how, and if this approximation is updated for every outer iteration of GMRES, or if it is computed once and for all.

Thanks in advance for your help.

I am using the following:

        jac = BroydenFirst()
        x0 = scipy.optimize.newton_krylov(fun, 
                                  x0, 
                                  ..
                                  inner_M=InverseJacobian(jac))

successfully, but I would like to better understand what is going on mathematically.