  We present a calculation of Third Harmonic Generation (THG) for two-band
systems using the length gauge that avoids unphysical divergences otherwise
present in the evaluation of the third order current density response. The
calculation is applied to bulk and monolayer black Phosphorus (bP) using a
non-orthogonal tight-binding model. Results show that the low energy response
is dominated by mixed inter-intraband processes and estimates of the magnitude
of THG susceptibility are comparable to recent experimental reports for bulk bP
samples.
