  We have derived an expression for the magnetic susceptibility of
topologically trivial insulators, however an important consideration for any
response tensor is whether it is gauge-invariant. By this we refer to the
gauge-freedom in choosing the Bloch functions one uses in a computation since
multiplication by an arbitrary $\textbf{k}$-dependent complex phase produces
equally valid Bloch functions. Additionally Wannier functions, which can be
constructed by a unitary mixing of the Bloch functions, are a useful basis for
computation due to their localized properties, but are strongly non-unique.
Therefore, we show that the theoretical expression we derived is independent of
the choice of Bloch functions or Wannier functions used, and thus is
gauge-invariant in this sense.
