How is the 2-point and 3-point jacobian implemented for BFGS minimization?

Question

I'd like to ask about the implementation or equation of the default 2-point jacobian and the 3-point jacobian for the BFGS minimization in scipy.optimize.

For context, I am using this to optimize crystal structure lattices.

I tried using it and had success with the 3-point jacobian. Here is a short input/output of the run.

result = optimize.minimize(fun = energy_function,
                           x0 = np.array([bulk_lc]),
                           method = 'BFGS',
                           jac = '3-point',
                           options = dict(disp = True,
                                         maxiter = 100,
                                         return_all = True,
                                         gtol = 1.e-8,
                                        )
                          )

BFGS
Optimization terminated successfully.
         Current function value: -56.362077
         Iterations: 2
         Function evaluations: 18
         Gradient evaluations: 6
Method:            BFGS
Message:           Optimization terminated successfully.
Iterations:        2
DFT calls:         18
Optimal input:     [4.0863827]
Function value:    -56.362077189626994
Tested value:      [array([4.046]), array([4.08419912]), array([4.0863827])]

The problem is that when I use 2-point for jacobian, it doesn't iterate at all and is stuck at the initial value for some reason.

May I ask for some insight on how the 2-point and 3-point jacobian works?