  High Harmonic Generation (HHG), which has been widely used in atomic gas, has
recently expanded to solids as a means to study highly nonlinear electronic
response in condensed matter and produce coherent high frequency radiation with
new properties. Most recently, attention has turned to Topological Materials
(TMs) and the use of HHG to characterize topological bands and invariants.
Theoretical interpretation of nonlinear electronic response in TMs, however,
presents many challenges. In particular, the Bloch wavefunction phase of TMs
has undefined points in the Brillouin Zone. This leads to singularities in
calculating the inter-band and intra-band transition dipole matrix elements of
Semiconductor Bloch Equations (SBEs). Here, we use the laser-electromagnetic
velocity gauge ${\boldsymbol p}\cdot {\bf A}(t)$ to numerically integrate the
SBEs and treat the singularity in the production of the electrical currents and
HHG spectra. We use a prototype of Chern Insulators (CIs), the Haldane model,
to demonstrate our approach. We find good qualitative agreement of the velocity
gauge compared to the length gauge and the Time-Dependent Density Functional
theory in the case of topologically trivial materials such as MoS$_2$. For
velocity gauge and length gauge, our two-band Haldane model reproduces key HHG
spectra features: ($\textit i$) The selection rules for linear and circular
light drivers, ($\textit ii$) The linear cut-off law scaling and ($\textit
iii$) The anomalous circular dichroism. We conclude that the velocity-gauge
approach captures experimental observations and provides theoretical tools to
investigate topological materials.
