The RHF magnetic inducibility
For the RHF wave function model, we can calculate the total (probability) current density as (at the converged, and current orbitals denoted by \(\boldsymbol{\kappa}^\star\)): \[\begin{multline} \require{physics} \vb{j} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}}, \vb{G}), \vb{B}_{\text{ext}}, \vb{G} ) \\ = - i \sum_i \Bigg( \phi_i^*(\vb{r}) \thinspace [ \grad{ \phi_i(\vb{r}) } ] - [ \grad{ \phi_i^*(\vb{r}) } ] \thinspace \phi_i(\vb{r}) + i \thinspace \phi_i^*(\vb{r}) \phi_i(\vb{r}) \thinspace \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}) \Bigg) \end{multline}\] if we do not employ any field-dependent spatial orbitals. We can also find this formula in (Stevens, Pitzee, and Lipscomb 1963) and (Keith and Bader 1993) in equation (1) (although they list the charge current density).
Without London orbitals
Derivation
The magnetic inducibility has implicit and explicit contributions. In this derivation, we will assume that the spin-orbitals are not London orbitals and are thus field-free.
Since the explicit perturbational derivative of the current density matrix elements are given by \[\begin{equation} \eval{ \pdv{ j_{u, \thinspace pq} ( \vb{r}; \vb{B}_{\text{ext}}, \vb{G} ) }{ B_{\text{ext}, m} } }_{\vb{B}_{\text{ext}, \thinspace 0}} = \frac{1}{2} \rho_{pq}(\vb{r}) \thinspace \epsilon_{umv} [\vb{r} - \vb{G}]_v \thinspace , \end{equation}\] the explicit contribution to the (first-order) magnetic inducibility can be calculated as: \[\begin{equation} \mathcal{J}_{um, \thinspace \text{explicit}} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{G} ) = \frac{1}{2} \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace \epsilon_{umv} [\vb{r} - \vb{G}]_v \thinspace , \end{equation}\]
The implicit part of the first-order induced response current density is given through a contraction of the first-order (response) derivative of the 1-DM and the field-free current density matrix elements. Specifically, for RHF theory, we apply linear response theory for a real perturbation-free optimization and a purely imaginary response and use the appropriate partial derivative of the RHF \(1\)-DM to find: \[\begin{equation} \mathcal{J}_{um, \thinspace \text{implicit}} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{G} ) = - 2i \sum_{ai} \Bigg( \phi_i(\vb{r}) [ \grad_u{ \phi_a(\vb{r}) } ] - [ \grad_u{ \phi_i(\vb{r}) } ] \thinspace \phi_a(\vb{r}) \Bigg) \thinspace (-i x_{ai, m} (\vb{G}) ) \thinspace . \end{equation}\] Summing up both contributions, we find that the RHF (first-order) magnetic inducibility can be written as: \[\begin{multline} \mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{G} ) \\ = - 2i \sum_{ai} \Bigg( \phi_i(\vb{r}) [ \grad_u{ \phi_a(\vb{r}) } ] - [ \grad_u{ \phi_i(\vb{r}) } ] \thinspace \phi_a(\vb{r}) \Bigg) \thinspace (-i {}^{\text{I}} x_{ai, m} (\vb{G}) ) + \frac{1}{2} \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace \epsilon_{umv} [\vb{r} - \vb{G}]_v \thinspace . \end{multline}\]
Continuous set of gauge transformations
Since the magnetic inducibility is still gauge origin dependent (both explicitly due to the last term and implicitly due to the linear response), we will now follow Keith and Bader (Keith and Bader 1993), using a continuous set of gauge transformations to resolve the gauge origin issue. Viewing the first-order magnetizability perturbationally dependent only on the external magnetic field \(\vb{B}_{\text{ext}}\) is incomplete. Rather, we must acknowledge that, if a magnetic field is present, a perturbation in the gauge origin may also cause a change in the first-order magnetizability: \[\begin{equation} \mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d} ) = \mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d}_0 ) + \mathcal{J}_{um, \thinspace \text{gauge}} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d} ) \thinspace . \end{equation}\] with \(\mathcal{J}_{um, \thinspace \text{gauge}}(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d})\) the contribution to the (first-order) magnetic inducibility due to the gauge translation \(\vb{d}\). If we introduce the corresponding response coefficient as \[\begin{equation} \Gamma_{umn} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}) ) = \eval{ \dv{ \mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d} )}{d_n} }_{\vb{d}_0} \thinspace , \end{equation}\] we may write the contribution to the magnetic inducibility due to the gauge translation as \[\begin{equation} \mathcal{J}_{um, \thinspace \text{gauge}} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d} ) = \sum_n \Gamma_{umn} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}) ) \thinspace d_n \thinspace . \end{equation}\] Since there is no gauge origin dependence in the field-free case, we know that \[\begin{equation} \eval{ \pdv{ {}^{\text{R}/\text{I}} \kappa_{ai}^\star ( \vb{B}_{\text{ext}, \thinspace 0} ) }{d_n} }_{\vb{d}_0} = 0 \end{equation}\] and thus the application of the chain rule in the calculation of \(\Gamma_{umn}(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}))\) only yields a contribution from the explicit partial derivative: \[\begin{align} &\Gamma_{umn} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}) ) \\ &\hspace{12pt} = \eval{ \pdv{ \mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d} )}{d_n} }_{\vb{d}_0} \\ &\hspace{12pt} = - 2i \sum_{ai} \Bigg( \phi_i(\vb{r}) [ \grad_u{ \phi_a(\vb{r}) } ] - [ \grad_u{ \phi_i(\vb{r}) } ] \thinspace \phi_a(\vb{r}) \Bigg) \thinspace (-i {}^{\text{I}} \varepsilon_{ai, mn}) - \frac{1}{2} \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace \epsilon_{umn} \thinspace , \end{align}\] where \({}^{\text{I}} \varepsilon_{ai, mn}\) is a mixed second-order response. Summing up the contributions, the total first-order magnetic inducibility can be calculated as \[\begin{align} &\mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{d} ) \\ &\hspace{12pt} = - 2i \sum_{ai} \Bigg( \phi_i(\vb{r}) [ \grad_u{ \phi_a(\vb{r}) } ] - [ \grad_u{ \phi_i(\vb{r}) } ] \thinspace \phi_a(\vb{r}) \Bigg) \thinspace \qty[ -i {}^{\text{I}} x_{ai, m} (\vb{G}_0) + \sum_n (-i {}^{\text{I}} \varepsilon_{ai, mn}) \thinspace d_n ] \\ &\hspace{24pt} + \frac{1}{2} \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace \epsilon_{umn} [\vb{r} - \vb{G}_0 - \vb{d}]_n \thinspace . \end{align}\] The CSGT approach (Keith and Bader 1993) then chooses the gauge translation \(\vb{d}\) to be different for each point \(\vb{r}\) where the first-order magnetizability is to be evaluated. More specifically, the choice \[\begin{equation} \vb{d}(\vb{r}) = \vb{r} - \vb{G}_0 \end{equation}\] is made, in order to make the density-related contribution vanish.
Connection with (Keith and Bader 1993)
We have now determined the (first-order) magnetic inducibility, which is the perturbational response coefficient of the (first-order) induced current density. In order to obtain the (first-order) induced current density, we only have to contract the magnetic inducibility with the corresponding perturbations.
If we only consider the external magnetic field (with a fixed gauge origin \(\vb{G}_0\)) as a perturbation, we find: \[\begin{align} &\vb{j}^{(1)} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}}, \vb{G}_0 ) \\ &\hspace{12pt} = \sum_m \mathcal{J}_{um} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{G}_0 ) \thinspace B_{\text{ext}, \thinspace m} \\ &\hspace{12pt} = - 2i \sum_{ai} \Bigg( \phi_i(\vb{r}) [ \grad{ \phi_a(\vb{r}) } ] - [ \grad{ \phi_i(\vb{r}) } ] \thinspace \phi_a(\vb{r}) \Bigg) \thinspace C^{(1)}_{ai} ( \vb{B}_{\text{ext}}, \vb{G}_0 ) + \frac{1}{2} \rho ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}) ) \thinspace \vb{B} \cross (\vb{r} - \vb{G}_0) \thinspace , \end{align}\] where \(C^{(1)}_{ai} (\vb{B}_{\text{ext}}, \vb{G}_0 )\) are first-order coefficients. If we then introduce the purely imaginary first-order perturbed orbitals as \[\begin{equation} \phi^{(1)}_i ( \vb{r}; \vb{B}_{\text{ext}}, \vb{G}_0 ) = \sum_a C^{(1)}_{ai} ( \vb{B}_{\text{ext}}, \vb{G}_0 ) \thinspace \phi_a(\vb{r}) \thinspace , \end{equation}\] we may finally write the first-order induced response current density as: \[\begin{align} &\vb{j}^{(1)} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}}, \vb{G}_0 ) \\ &\hspace{12pt} = - 2i \sum_i \Bigg( \phi_i(\vb{r}) \thinspace [ \grad{ \phi^{(1)}_i ( \vb{r}; \vb{B}_{\text{ext}}, \vb{G}_0 )} ] - [ \grad{ \phi_i(\vb{r}) } ] \thinspace \phi^{(1)}_i ( \vb{r}; \vb{B}_{\text{ext}}, \vb{G}_0 ) \Bigg) \\ &\hspace{24pt} + \frac{1}{2} \thinspace \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace [ \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}_0) ] \thinspace , \end{align}\] which is exactly the result found in (Stevens, Pitzee, and Lipscomb 1963) on p. 554 and the result found in (Keith and Bader 1993) in equation (3) in the case that \(\vb{G}_0 = \vb{0}\) (although they list the charge current density).12
If the gauge origin translation \(\vb{d}\) is also regarded as a perturbation, we must also calculate the associated contribution to the (first-order) induced response current density. Using a similar contraction procedure and realizing that the first-order perturbed orbitals become \[\begin{equation} \phi^{(1)}_i ( \vb{r}; \vb{B}_{\text{ext}}, \vb{d} ) = \sum_a \qty[ C^{(1)}_{ai} ( \vb{B}_{\text{ext}}, \vb{G}_0 ) + \delta C^{(1)}_{ai} ( \vb{B}_{\text{ext}}, \vb{d} ) ] \thinspace \phi_a(\vb{r}) \thinspace , \end{equation}\] where \(\delta C^{(1)}_{ai} ( \vb{B}_{\text{ext}}, \vb{d} )\) are the changes in the first-order coefficients, we may write the (first-order) induced response current density as \[\begin{align} &\vb{j}^{(1)} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0}), \vb{B}_{\text{ext}}, \vb{d} ) \\ &\hspace{12pt} = - 2i \sum_i \Bigg( \phi_i(\vb{r}) \thinspace [ \grad{ \phi^{(1)}_i ( \vb{r}; \vb{B}_{\text{ext}}, \vb{d} ) }] - [ \grad{ \phi_i(\vb{r}) } ] \thinspace \phi^{(1)}_i ( \vb{r}; \vb{B}_{\text{ext}}, \vb{d} ) \Bigg) \\ &\hspace{24pt} + \frac{1}{2} \thinspace \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace [ \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}_0 - \vb{d}) ] \thinspace . \end{align}\] Specifically, when applying the magnetic field along the \(x\)-axis, we obtain \[\begin{align} &\vb{j} ( \vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_0), \vb{B}, \vb{d} ) \\ &\hspace{12pt} = -2i B \sum_i \Bigg( \phi_i(\vb{r}) [ \grad{ \phi_i^{(1), \thinspace L_x} (\vb{r}; \vb{G}_0) } ] - [ \grad{ \phi_i(\vb{r})} ] \phi_i^{(1), \thinspace L_x} (\vb{r}; \vb{G}_0) \Bigg) \\ &\hspace{24pt} +2i B d_y \sum_i \Bigg( \phi_i(\vb{r}) [ \grad{ \phi_i^{(1), \thinspace p_z} (\vb{r})} ] - [ \grad{ \phi_i(\vb{r})} ] \phi_i^{(1), \thinspace p_z} (\vb{r}) \Bigg) \\ &\hspace{24pt} -2i B d_z \sum_i \Bigg( \phi_i(\vb{r}) [ \grad{ \phi_i^{(1), \thinspace p_y} (\vb{r}) } ] - [ \grad{ \phi_i(\vb{r})} ] \phi_i^{(1), \thinspace p_y} (\vb{r}) \Bigg) \\ &\hspace{24pt} + \frac{1}{2} \thinspace \rho(\vb{r}; \boldsymbol{\kappa}^\star(\vb{B}_{\text{ext}, \thinspace 0})) \thinspace [ \vb{B}_{\text{ext}} \cross (\vb{r} - \vb{G}_0 - \vb{d}) ] \thinspace , \end{align}\] which is equation (13) in (Keith and Bader 1993) (although they list the charge current density). The realizations of the (purely imaginary) perturbed orbitals are as follows: \[\begin{equation} \phi_i^{(1), \thinspace L_m} ( \vb{r}; \vb{G}_0 ) = \sum_a (-i {}^{\text{I}} x_{ai, m}(\vb{G}_0)) \thinspace \phi_a(\vb{r}) \end{equation}\] and \[\begin{equation} \phi_i^{(1), \thinspace p_o} ( \vb{r} ) = -\epsilon_{mno} \sum_a (-i {}^{\text{I}} \varepsilon_{ai, mn}) \thinspace \phi_a(\vb{r}) \thinspace . \end{equation}\]
References
Keith and Bader (Keith and Bader 1993) have called the term involving the density the diamagnetic contribution but, as far as I am aware, this has no connection with the diamagnetic part of the scalar kinetic operator. In fact, following the derivation of the Pauli current density, it is actually the paramagnetic contribution to the scalar kinetic operator which gives rise to this term.↩︎