\section{Magnetic response and gauge invariance in real-space grids}
%Miguel, Xavier

In the presence of a magnetic field \(\vec{B}(\vec{r}, t)\), generated
by a vector potential \(\vec{A}(\vec{r}, t)\), additional terms describing the coupling of 
the electrons to the magnetic field must be included in the Hamiltonian
%
\begin{equation}
  \label{eq:magnetic}
  \hat{H} = \frac12\left({\hat{\vec{p}}} -\frac1c{\vec{A}}\right)^2 + \hat{v} + \vec{B}\cdot{\hat{\vec{S}}}\ .
\end{equation}
%
The first part describes the orbital interaction with the field, and the
second one is the Zeeman term that represents the coupling of the
electronic spin with the magnetic field.

As our main interest is the evaluation of the magnetic susceptibility, 
in the following, we consider a perturbative uniform static magnetic
field \(\vec{B}\) applied to a finite system with zero total spin. In the Coulomb gauge the
corresponding vector potential, \(\vec{A}\), is given as
%
\begin{equation}
  \label{eq:unifmag}
  \vec{{A}}(\vec{r}) = \frac12\vec{B}\times{\vec{r}}\ .
\end{equation}
%
In orders of \(\vec{B}\) the perturbing potentials are
%
\begin{equation}
  \label{eq:mag1}
  \delta\hat{v}_i^{\mathrm{mag}} = \frac1{2c}({\vec{r}}\times{\hat{\vec{p}}})_i = \frac1{2c}\hat{L}_i\ ,
\end{equation}
%
with \(\vec{\hat{L}}\) the angular momentum operator, and
%
\begin{equation}
  \label{eq:mag2}
  \delta^2\hat{v}_{ij}^{\mathrm{mag}} = \frac1{8c^2}(\delta_{ij}r^2 - r_i r_j)\ .
\end{equation}

The induced magnetic moment can be expanded in terms of the external magnetic field which, to first order, reads
%
\begin{equation}
  \label{eq:bchi}
  m_i=m^0_i+\sum_j\chi_{ij}B^{ext}_j\ ,
\end{equation}
%
where \(\vec{\chi}\) is the magnetic susceptibility tensor. For finite systems the
permanent magnetic moment can be calculated directly from the
ground-state wave-functions as
%
\begin{equation}
  \label{eq:magneticmoment}
  \vec{m}^0=
  \sum_n\langle\varphi_n|\delta{\hat{\vec{v}}}^{\mathrm{mag}}|\varphi_{n}\rangle\ .
\end{equation}
%
For the susceptibility, we need to calculate the first-order response
functions in the presence of a magnetic field. This can be done in
practice by using the magnetic perturbation, eq.~(\ref{eq:mag1}), in
the Sternheimer formalism described in 
section \ref{sec:sternheimer}. If the system is time-reversal symmetric, since the perturbation is anti-symmetric under time-reversal, it does not induce a change
in the density and the Sternheimer equation does not need to be solved
self-consistently. From there we find
%
\begin{multline}
  \label{eq:susc}
  \chi_{ij} = \sum_n\Big[
    \langle{\varphi_n}|\delta\hat{v}^{\mathrm{mag}}_j|{\delta\varphi_{n,i}}\rangle + {\rm cc.}
    %+\langle{\delta\varphi_{n\,,\,i}}|\delta\hat{v}^{\mathrm{mag}}_j|{\varphi_{n}}\rangle\\
    +\langle{\varphi_n}|\delta^2\hat{v}^{\mathrm{mag}}_{ij}|{\varphi_{n}}\rangle
  \Big]\ .
\end{multline}
%
Before applying this formalism in a calculation, however, we must make sure that our calculation is gauge invariant.

In numerical implementations, the gauge freedom in choosing the vector
potential might lead to poor convergence with the quality of the
discretization, and to a dependence of the magnetic response on the
origin of the simulation cell. In other words, an arbitrary
translation of the molecule could introduce an nonphysical change in
the calculated observables. This broken gauge-invariance is well known
in molecular calculations with all-electron methods that make use of
localized basis sets. In this case,
the error can be traced to the finite-basis-set representation of the
wave-functions~\cite{Wolinski_1990,Schindler_1982}. A simple measure of the error
is to check for the fulfillment of the hyper-virial
relation~\cite{Bouman_1977}.
%
\begin{equation}
  \label{eq:hypervirial}
  i\langle\varphi_j\vert{\hat{\vec{p}}}\vert\varphi_n\rangle = 
  (\epsilon_n-\epsilon_j)\langle\varphi_j\vert{\hat{\vec{r}}}\vert\varphi_n\rangle\ ,
\end{equation}
where \(\epsilon_n\) is the eigenvalue of the state \(\varphi_n\).

When working with a real-space mesh, this problem also appears, though
it is milder, because the standard operator representation in the grid is
not gauge-invariant. In this case the error can be controlled by
reducing the spacing of the mesh. On the other hand, real-space grids usually require the use of the pseudo-potential approximation, where the electron-ion interaction is
described by a non-local potential
\(\hat v_{\text{nl}}\). This, or any other non-local
potential, introduces a fundamental problem when describing the
interaction with magnetic fields or vector potentials in general. To
preserve gauge invariance, this term must be adequately coupled to the
external electromagnetic field, otherwise the results will strongly
depend on the origin of the gauge. For example, an extra term has to be
included in the hyper-virial expression, eq.~(\ref{eq:hypervirial}), resulting in
%
\begin{equation}
 \label{eq:hypervirial_nl}
 i\langle\varphi_j\vert\hat{\vec{p}}\vert\varphi_n\rangle = 
 (\epsilon_n-\epsilon_j)\langle\varphi_j\vert\hat{\vec{r}}\vert\varphi_n\rangle +
 \langle\varphi_j\vert[\hat{\vec{r}},\hat{v}_{\text{nl}}]\vert\varphi_n\rangle\ .
\end{equation}
%

In general, the gauge-invariant non-local potential is given by
%
\begin{equation}
 \label{eq:gauge-methods}
 \langle \vec{r} | \hat{v}^{\vec{{A}}}_{\text{nl}} | \vec{r}' \rangle =
 \langle \vec{r} | \hat{v}_{\text{nl}}| \vec{r}' \rangle \mathrm{exp}\left(\frac{i}{c}\int_{\vec{r}}^{\vec{r}'} \vec{A}(\vec{x},t) \cdot \mathrm{d}\vec{x}\right)
 \ .
\end{equation}
The integration path can be any one that connects the two
points \(\vec{r}\) and \(\vec{r}'\), so an infinite number of choices
is possible. 

In order to calculate the corrections required to the magnetic perturbation operators, we use two different integration paths that have been
suggested in the literature. The first was proposed by Ismail-Beigi,
Chang, and Louie (ICL)~\cite{Ismail_Beigi_2001} who give the
following correction to the first-order magnetic perturbation term
%
\begin{equation}
  \label{eq:iclmag1}
  \delta \hat{\vec{v}}^{\mathrm{ICL}} = \delta \vec{\hat{v}}^{\mathrm{mag}} - \frac{i}{2c}\vec{\hat{r}}\times\left[\hat{\vec{r}}, v_{\text{nl}}\right]\ ,
\end{equation}
%
and a similar term for the second-order perturbation. Using a different integration path, Pickard and Mauri~\cite{Pickard_2003} proposed the GIPAW method, that has the form
%GIPAW
%
\begin{equation}
  \label{eq:gipawmag1}
  \delta \hat{\vec{v}}^{\mathrm{GIPAW}} = \delta \hat{\vec{v}}^{\mathrm{mag}} - \frac{i}{2c}\sum_\alpha\vec{R}_\alpha\times\left[\hat{\vec{r}}, \hat{v}^\alpha_{\text{nl}}\right]\ ,
\end{equation}
%
where \(\vec{R}_\alpha\) and \(\hat{v}^\alpha_{\text{nl}}\) are, respectively, the position and non-local potential of atom \(\alpha\).
%
With the inclusion of either one of these methods, both implemented in Octopus, we recover gauge invariance in our formalism when pseudo-potentials are used. This allows us to predict the magnetic susceptibility and other properties that depend on magnetic observables, like optical activity~\cite{Varsano_2009}.

A class of systems with interesting magnetic susceptibilities are fullerenes. For example, it is known that the C\(_{60}\) fullerene has a very small magnetic susceptibility due to the cancellation of the paramagnetic and diamagnetic responses~\cite{Haddon_1991,Haddon_1995}. Botti~\emph{et al.}~\cite{Botti_2009} used the real-space implementation of Octopus to study the magnetic response of the boron fullerenes depicted in Fig.~\ref{fig:boronfullerenes}. As shown in table~\ref{tab:boronmagnetic}, they found that, while most clusters are diamagnetic, B\(_{80}\) is paramagnetic, with a strong cancellation of the paramagnetic and diamagnetic terms.

\begin{table}
\centering
\begin{tabular}{lc}
\mbox{} & $\bar\chi$ \\ \hline
B$_{20}$     & -250.2  \\
B$_{38}$     & -468.3  \\
B$_{44}$     & -614.4  \\
B$_{80}$     &  219.3  \\
B$_{92}$     & -831.3 
\end{tabular}
\caption{
\label{tab:boronmagnetic}
Calculated magnetic susceptibilities ($\chi$
in cgs\,ppm/mol) per number of boron atoms
for the selected boron clusters shown in Fig.~\ref{fig:boronfullerenes}. Results from Ref.~\cite{Botti_2009}}
\end{table}