I had initially performed tests on an 8x8x8 k-grid and obtained a positive Im[Chi^3] in the mentioned energy range. However, I’ve now performed tests in denser k-grids and obtained something similar to what you have. Both 24x24x24 and 32x32x32 results show that Im[Chi^3] becomes negative at the energy of the gap and positive afterwards. So, it does not surprise me that your case (28x28x28) shows the same. See example for Silicon 32x32x32.
I do not think there is anything unphysical about Im[Chi^3] being negative. You can find examples of it in the literature, be it related to simulated (Stokes) Raman scattering, or self-defocusing nonlinearities (negative n_2). Also the paper you cite for the expressions for beta and n_2 mentions a negative Im[Chi^3] in connection to absorption saturation. In relation to the cited examples (liquid water and h-BN) where Im[chi^3] is positive, I would say they are wide-gap systems, so I would expect significant differences with Si or GaAs.
Finally, let me give a minor recommendation for your calculations going forwards. I’ve noticed you use bands “1 | 40 |”, which seems to be default you get from Abinit. In my experience, you can get away with using much fewer bands for Silicon (e.g., 8), which will significantly reduce the computational burden in CPU time and memory (I saw your runs took over 2 days, but Silicon is supposedly a ‘small’ system). Ultimately, the number of bands included should be determined though convergence tests. Also, be sure to perform convergence tests for NLtime and NLstep. My tests indicate that you could safely use half the NLtime you are using, which would halve the CPU time required.
Hope that helps,