The projected Newton method is currently using the constraint JF * d = 0 instead of F + JF * d = 0. We can modify the current Lagrangian system to have this constraint instead, which may have some significant impact on later convergence.
It is worth revisiting this problem and trying to see if convergence improves with this constraint.