# Magnetic response of metallic nanoparticles: geometrical and weakly-relativistic effects

Mauricio Gómez Viloria, Guillaume Weick, Dietmar Weinmann, Rodolfo A. Jalabert

► To cite this version:

Mauricio Gómez Viloria, Guillaume Weick, Dietmar Weinmann, Rodolfo A. Jalabert. Magnetic response of metallic nanoparticles: geometrical and weakly-relativistic effects. Physical Review B, 2021, 104, pp.245428. 10.1103/PhysRevB.104.245428. hal-03281633

HAL Id: hal-03281633

https://hal.science/hal-03281633v1

Submitted on 26 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Magnetic response of metallic nanoparticles: Geometric and weakly relativistic effects

Mauricio Gómez Viloria, Guillaume Weick, Dietmar Weinmann, and Rodolfo A. Jalabert

Université de Strasbourg, CNRS, Institut de Physique et Chimie des Matériaux de Strasbourg, UMR 7504, F-67000 Strasbourg, France

While the large paramagnetic response measured in certain ensembles of metallic nanoparticles has been assigned to orbital effects of conduction electrons, the spin-orbit coupling has been pointed out as a possible origin of the anomalously large diamagnetic response observed in other cases. Such a relativistic effect, arising from the inhomogeneous electrostatic potential seen by the conduction electrons, might originate from the host ionic lattice, impurities, or the self-consistent confining potential. Here we theoretically investigate the effect of the spin-orbit coupling arising from the confining potential, quantifying its contribution to the zero-field magnetic susceptibility and gauging it against the ones generated by other weakly-relativistic corrections. Two ideal geometries are considered in detail, the sphere and the half-sphere, focusing on the expected increased role of the spin-orbit coupling upon a symmetry reduction, and the application of these results to actual metallic nanoparticles is discussed. The matrix elements of the different weakly-relativistic corrections are obtained and incorporated in a perturbative treatment of the magnetic field, leading to tractable semi-analytical and semiclassical expressions for the case of the sphere, while a numerical treatment becomes necessary for the half-sphere. The correction to the zero-field susceptibility arising from the spin-orbit coupling in a single sphere is quite small, and it is dominated by the weakly-relativistic kinetic energy correction, which in turn remains considerably smaller than the typical values of the nonrelativistic zero-field susceptibility. Moreover, the spin-orbit contribution to the average response for ensembles of nanoparticles with a large size dispersion is shown to vanish. The symmetry reduction in going from the single sphere to the half-sphere does not translate into a significant increase of the spin-orbit contribution to the zero-field susceptibility.

# I. INTRODUCTION

The spin-orbit coupling (SOC) is a relativistic effect having a decisive role in certain properties of unconfined and confined condensed matter systems. A celebrated example is the spin-orbit driven change from weak localization to anti-weak localization observed in the electronic transport through metallic films [1] or ballistic quantum dots at a semiconductor heterojunction [2, 3]. In ferromagnetic materials the SOC underlines the phenomena of magnetic anisotropy and the anomalous Hall effect [4]. The domain of spintronics addresses numerous cases where the SOC influences spin dynamics and spin relaxation [5-7]. At the fundamental level, the SOC is important because it alters the symmetry properties of the electronic Hamiltonian. As a consequence, the statistics of energy levels in time-reversal symmetric (chaotic or disordered) confined systems changes from orthogonal to symplectic (Gaussian or circular) distributions when going from vanishing to strong SOC [8, 9].

The tunneling resonances of disordered metal nanoparticles and the magnetic response of an ensemble of metallic nanoparticles are two examples of physical properties depending on the level statistics of a confined system, and thus on the SOC.

In the first case, the extracted *g*-factor of the discrete energy levels was found [10, 11] to be below the free-electron value of $g_0 = 2$. Such a reduction could be explained by the fact that the energy eigenstates in the presence of SOC, not being purely spin up or spin down, respond more weakly to an applied magnetic field than pure spin states [12, 13]. The statistical distribution of the *g*-factors has been obtained from random matrix theory, using the spin-orbit scattering rate as a phenomenological parameter in order to describe the transition between statistical ensembles [14, 15].

In the second case, the zero-field susceptibility (ZFS) aver-

aged over a nanoparticle ensemble is determined [16, 17] by the magnetic field dependence of the variance in the number of energy levels below the chemical potential [see Eq. (10c) below]. Therefore, in the disordered or chaotic regimes, it depends on the transition between statistical ensembles driven by the influence on the SOC. In particular, the SOC has been invoked to be responsible for the large diamagnetism measured in an ensemble of gold nanorods [18, 19]. This is a very interesting proposal, since among the anomalous magnetic responses observed in ensembles of metallic nanoparticles, only the paramagnetic behavior [20-28] has been accounted for [29], while no satisfying theory existed to describe the observed [18, 22, 26, 30, 31] large diamagnetic response in metallic nanoparticles.<fnref n="1" />,<fnref n="2" />

In the previously presented cases of the tunneling resonances of disordered metal nanoparticles and the magnetic response of an ensemble of metallic nanoparticles, the strength of the SOC is a key parameter that needs to be determined by microscopic theories. Towards this goal, the main source of SOC must be identified. The genesis of the SOC for the conduction electrons of a metallic nanoparticle lies in the existence of an inhomogeneous electrostatic potential, which may have an intrinsic origin (the host ionic lattice) or an extrinsic origin (impurities or the confining potential). Since gold is a

<footnote n="1">¹ Superconducting fluctuations that persist beyond the critical temperature were shown to result in a relatively large diamagnetic response [32], which nevertheless remains one to two orders of magnitude smaller than the one reported in the experiments of Ref. [18].</footnote>

<footnote n="2">² The most spectacular among the anomalous responses, which is the ferromagnetic behavior [25, 26, 28, 30, 31, 33-39], requires to have ferromagnetism at the level of individual nanoparticles or through the interparticle dipolar interaction [29, 40], and will not be addressed in this work. (See, e.g., Ref. [41] for a discussion of the different magnetic behaviors exhibited by ensembles of gold nanoparticles.)</footnote>
heavy atom, SOC plays an important role in its band structure [42], but the effect for the conduction electrons is mainly seen by the Bloch part of the wave function, while the smooth part remains unaffected. This observation is consistent with the g-factor $g_{Au} = 2.1$ measured by electron spin resonance in macroscopic gold samples [43] and the fact that the spin diffusion length in gold is typically limited by extrinsic effects [44]. Impurities have been invoked to be responsible for the SOC of Ag nanoparticles intentionally doped with Au [10], but they are expected to play a lesser role in ballistic nanoparticles where the magnetic susceptibility has been measured. In this last setup, the electronic confinement remains as the chief source of SOC. The effect of this latter mechanism in the ZFS is the goal of this paper, where we analyze model systems of noninteracting electrons with different kinds of confinement, making the link with experimentally relevant cases of metallic nanoparticles.

The SOC yields a contribution to the fine-structure of atomic spectra which is of the same order as those arising from other weakly-relativistic corrections, namely the kinetic energy and Darwin terms [45]. Hence, we contrast the SOC with the previously cited weakly-relativistic corrections, as well as with the angular magneto-electric coupling relevant at finite magnetic field, for a confined electron gas in the limit of Fermi velocities $v_F$ much smaller than the speed of light in vacuum c.

The ZFS of a nonrelativistic, three dimensional, unconfined (bulk) degenerate electron gas is $\chi_b^{(nr)} = \chi_L + \chi_P$, resulting from the combined effect of the diamagnetic Landau susceptibility [46]

$$ \chi_L = -\frac{1}{12\pi^2} \frac{e^2 k_F}{m^* c^2}, \qquad (1) $$

arising from the orbital motion, and the paramagnetic Pauli susceptibility

$$ \chi_P = 3 |\chi_L|, \qquad (2) $$

originating from the Zeeman effect over the electron spin. We note $-e$ the electron charge, $k_F$ the Fermi wave vector, and $m^*$ the effective mass. In metals, the difference between $m^*$ and the free electron mass $m$ is very small (i.e., $m_{Au}^* = 1.1 m$) and therefore we will identify the two electron masses.<fnref n="3" /> Throughout the paper, we use cgs units, which for the case of gold result in a magnetic susceptibility $\chi_L = -2.9 \times 10^{-7}$. The ZFS of macroscopic gold samples is $\chi_{Au} \simeq -9.3 |\chi_L|$, as the contribution $\chi_b^{(nr)}$ from the conduction electrons is dominated by the Larmor diamagnetic response arising from the closed-shell ion-core electrons [47].

The geometrical confinement of a metal to the nano- or micro-scales does not affect the response of the core electrons, while the effect on conduction electrons might be important,

leading to a diamagnetic or paramagnetic ZFS with typical values considerably larger than $|\chi_L|$, and the Larmor susceptibility as well, depending on the temperature, the Fermi wave vector, and the size of the containing box. The enhancement of the magnetic response is a consequence of the orbital motion of the conduction electrons [29], and can be estimated with the help of semiclassical expansions (see Appendix A).

In the relativistic case the diamagnetic and paramagnetic contributions cannot be separated. The bulk ZFS of a weakly-relativistic electron gas is given, in the limit $v_F \ll c$, by [48]<fnref n="4" />

$$ \chi_b^{(wr)} = \left[ 2 - \frac{1}{3} \left( \frac{v_F}{c} \right)^2 \right] |\chi_L|. \qquad (3) $$

Thus, in normal metals the weakly-relativistic correction with respect to $\chi_b^{(nr)}$ is very small. Given the dramatic increase of the ZFS with respect to $\chi_L$ induced by the electronic confinement in the nonrelativistic case, we might ask if a similar effect occurs for the weakly-relativistic susceptibility, and in particular for the contribution arising from the SOC. This is the goal of this work, where we analyze the case of a model system and apply it to metallic nanoparticles.

Our paper is organized as follows: In Sec. II we recall the thermodynamical quantum-mechanical and semiclassical formalisms which enable us to obtain the ZFS of a confined system. In Sec. III we present our model which we employ in order to assess the relevance of weakly-relativistic effects for the ZFS of spherical (Sec. IV) and half-spherical nanoparticles (Sec. V). We conclude in Sec. VI. Six appendixes presenting some of the details of our calculations, as well as alternative semiclassical derivations, complete the manuscript.

## II. MAGNETIC RESPONSE OF A FINITE FERMION SYSTEM

The magnetic moment of a system of noninteracting fermions in a volume $V$, at a temperature $T$, and with a chemical potential $\mu$ is given by [50, 51]

$$ \begin{aligned} \mathcal{M} &= -\frac{\partial \Omega}{\partial H} \\ &= -\sum_{\{\lambda\}} f_{\mu}(E_{\lambda}) \frac{\partial E_{\lambda}}{\partial B}, \end{aligned} \qquad (4) $$

with

$$ \Omega(\mu, T, B) = -k_B T \int_0^\infty dE \, \varrho(E, B) \ln \left( 1 + e^{(\mu - E)/k_B T} \right) \qquad (5) $$

the thermodynamic potential associated with the grand canonical ensemble (GCE), $f_\mu(E) = \{\exp([E - \mu]/k_B T) + 1\}^{-1}$

<footnote n="3">³ The situation considerably changes when going from the case of metals to that of semiconductors, where such an identification is not valid and Eq. (2) does not hold.</footnote>

<footnote n="4">⁴ Other values of the ZFS have been proposed in the literature [49]. However, the result (3) is easily obtainable in the grand canonical ensemble by using the relativistic Landau levels, and it is in agreement with the numerical calculations that we present in Sec. IV C.</footnote>
the Fermi-Dirac distribution, and $k_B$ the Boltzmann constant. The eigenenergies $E_\lambda$ of the system are characterized by the set of quantum numbers $\{\lambda\}$, while $\varrho(E, B) = \sum_{\{\lambda\}} \delta(E - E_\lambda)$ is the field-dependent single-particle density of states (DOS). Here, $\delta(\zeta)$ denotes the Dirac delta function. In the regime of a small magnetic response of the system that we are interested in, the magnitude of the magnetic field $H$ or the magnetic induction $B$ can be indistinctly used when taking the derivatives in Eq. (4).

The ZFS is obtained from Eq. (4) as

$$ \chi = \frac{1}{\bar{V}} \left. \frac{\partial \mathcal{M}}{\partial H} \right|_{H=0} \\ = -\frac{1}{\bar{V}} \sum_{\{\lambda\}} \left[ f'_{\mu}(E_{\lambda}) \left( \frac{\partial E_{\lambda}}{\partial B} \right)^2 + f_{\mu}(E_{\lambda}) \frac{\partial^2 E_{\lambda}}{\partial B^2} \right] \Big|_{B=0} . \quad (6) $$

Note that the unconstrained ZFS of the bulk, discussed in the introduction, follows from Eq. (6) while taking the limit $\mathcal{V} \to \infty$. This last procedure is not exempted of some subtleties [52–54].

If the fermion system is not able to exchange particles with a reservoir, and then has a fixed number of particles $N$ (as is the case of a metallic nanoparticle), the free energy $F$, associated with the canonical ensemble (CE), should be used instead of the thermodynamic potential $\Omega$ [cf. Eq. (5)] in the definitions (4) and (6) of the $N$-fixed magnetic moment and ZFS, $M_N$ and $\chi_N$, respectively. For large $N$ the difference between $\chi$, evaluated in the GCE at the Fermi energy $E_F$, and $\chi_N$ is generally very small. However, such a correction becomes essential whenever the GCE magnetic response vanishes, as is the case of the average ZFS over an ensemble of nanoparticles with a large size dispersion [29].

The above corrections can be readily incorporated in a semiclassical formalism in which the DOS is decomposed into a smooth (Weyl) and an oscillating (in energy or over the nanoparticle ensemble) part as

$$ \varrho(E, B) = \bar{\varrho}(E, B) + \varrho^{\text{osc}}(E, B). \quad (7) $$

The trace formula expresses $\varrho^{\text{osc}}(E, B)$ as a sum running over the periodic trajectories $\xi$ of the system,

$$ \varrho^{\text{osc}}(E, B) = \sum_{\xi} \varrho_{\xi}(E, B), \quad (8) $$

where the specific form of the contributions $\varrho_\xi(E, B)$ depends on whether $\xi$ represents isolated periodic orbits or degenerate families of periodic orbits [55]. For temperatures such that $k_B T$ is larger than the typical level spacing, the free energy can be approximated as [17, 53]

$$ F(N, T, B) \approx \bar{F} + \Delta F^{(1)} + \Delta F^{(2)}, \quad (9) $$

where

$$ \bar{F} = \bar{\Omega}(\bar{\mu}, T, B) + \bar{\mu}N, \quad (10a) $$

$$ \Delta F^{(1)} = \Omega^{\text{osc}}(\bar{\mu}, T, B), \quad (10b) $$

$$ \Delta F^{(2)} = \frac{1}{2\bar{\varrho}(\bar{\mu}, 0)} \left[ \int_{0}^{\infty} dE \, \varrho^{\text{osc}}(E, B) f_{\bar{\mu}}(E) \right]^2 . \quad (10c) $$

In the equations above, $\bar{\Omega}$ and $\Omega^{\text{osc}}$ are obtained, respectively, when $\bar{\varrho}$ and $\varrho^{\text{osc}}$ are used in Eq. (5) instead of $\varrho$. The mean chemical potential $\bar{\mu}$ is determined in such a way as to ensure accommodating the $N$ electrons while using the mean density of states $\bar{\varrho}(E, 0)$, that is,

$$ N = \int_{0}^{\infty} dE \, \bar{\varrho}(E, 0) f_{\bar{\mu}}(E). \quad (11) $$

When only the orbital motion is considered, the term $\bar{F}$ of Eq. (10a) leads to the Landau susceptibility presented in Eq. (1). The calculation of $\Omega^{\text{osc}}$ in Eq. (10b) involves the energy integration of rapidly oscillating functions of $E$ (appearing through the classical action of the periodic orbits) that results in [53]

$$ \chi^{\text{osc}} = -\frac{1}{\bar{V}} \sum_{\xi} R(\tau_{\xi}/\tau_{T}) \left(\frac{\hbar}{\tau_{\xi}}\right)^2 \left. \frac{\partial^2 \varrho^{\text{osc}}(E_F, B)}{\partial B^2} \right|_{B=0}, \quad (12) $$

where $\hbar$ is Planck's (reduced constant), $\tau_\xi$ is the period of the classical orbit, and $\tau_T = \hbar/\pi k_B T$. The thermal factor

$$ R(\zeta) = \frac{\zeta}{\sinh(\zeta)} \quad (13) $$

exponentially cuts off the long trajectories of the semiclassical expansion. The two previously discussed contributions to $\chi_N$, arising from $\bar{F}$ and $\Delta F^{(1)}$, give together a good account of the ZFS for an individual nanoparticle [29], when compared to a quantum-mechanical perturbative (in the magnetic field $B$) calculation [50, 51]. The magnetic field dependence of $\bar{F}$ is much weaker than that of $\Delta F^{(1)}$ (which is given by the magnetic flux enclosed by the periodic orbits). Thus, $\chi^{\text{osc}}$ (related with the confinement) dominates over the bulk contribution (1). The finite-$N$ correction arising from $\Delta F^{(2)}$ in Eq. (10c) becomes crucial when considering an ensemble of nanoparticles, and it is important in order to assess the relevance of the relativistic corrections that we study in this work.

The above formulation has been applied to the problems of the persistent currents in metallic rings [16, 17, 56–58] and the orbital magnetism in low-dimensional ballistic systems [53, 54] for weakly-interacting Landau quasiparticles, as well as including the particularly relevant Cooper-like electronic correlations [59, 60]. In the absence of spin-orbit interaction, the DOS oscillations characterized by $\varrho^{\text{osc}}$ [i.e., Eq. (8)] exhibit a diminishing amplitude as the magnetic field is turned on, in line with the enhanced spectral rigidity associated with the transition from orthogonal to unitary Gaussian ensembles describing chaotic or disordered systems. The dependence of measurable physical effects on the transition between statistical ensembles, discussed in the introduction, appears in the case of the average response of a large number of nanostructures, where the contribution arising from $\Delta F^{(1)}$ vanishes, and thus $\Delta F^{(2)}$ provides the leading response. The average ZFS resulting from Eqs. (4), (6), and (10c) is then expected to be paramagnetic. On the one hand, the previous reasoning is confirmed by the paramagnetic ZFS measured in ensembles of $10^5$ phase-coherent, ballistic squares lithographically defined on a high mobility GaAs heterojunction [61].
On the other hand, the orbital response of $10^5$ microscopic silver rings was found to be diamagnetic [62]. The resulting discrepancy with the previous theory has not been settled, and points to the possible role of attractive electron-electron interaction above the superconducting transition temperature, as proposed in Refs. [32, 63].<fnref n="5" />

For a system with a strong SOC, the magnetic field drives the transition from symplectic to unitary Gaussian ensembles. From the larger spectral rigidity of the former, one might conclude that there is an SOC-induced inversion of the average magnetic response. However, as pointed out in Ref. [19], such a reasoning does not always hold, since in the presence of a strong SOC the breaking of the Kramers degeneracy occurs at magnetic fields much smaller than those lifting the Zeeman splitting in the absence of SOC. While Ref. [19] predicted, within a noninteracting model, a switch from paramagnetic to diamagnetic ZFS in the case of metallic nanoparticles, the average persistent current in disordered wires was shown to retain its paramagnetic character even if a strong SOC is present for noninteracting electrons [58, 65], as well as in the case of a repulsive electron-electron interaction [56, 59]. The nontrivial effect of SOC and the angular magneto-electric coupling for the persistent current in conducting rings has also been analyzed in terms of a Berry phase [66]. The persistent currents in ballistic mesoscopic rings wish Rashba SOC have been studied [67]. The coupling strength of the previous mechanism, often taken as a phenomenological parameter, has been linked with measurable features of the flux dependence of the persistent currents.

### III. WEAKLY-RELATIVISTIC DESCRIPTION

The weakly-relativistic Hamiltonian for an electron subject to the action of static electric and magnetic induction fields, respectively $\mathbf{E}(\mathbf{r}) = -\nabla\phi(\mathbf{r})$ and $\mathbf{B}(\mathbf{r}) = \nabla \times \mathbf{A}(\mathbf{r})$ (with $\phi$ and $\mathbf{A}$ the scalar and vector potentials, respectively), is given as a $2 \times 2$ matrix by [45, 68]

$$ \mathcal{H}^{(\text{wr})} = \mathcal{H}^{(\text{nr})} + \Delta\mathcal{H}. \quad (14) $$

The nonrelativistic Hamiltonian $\mathcal{H}^{(\text{nr})}$ is

$$ \mathcal{H}^{(\text{nr})} = \mathcal{H}^{(\text{orb})} + \mathcal{H}^{(\text{Z})}, \quad (15) $$

with

$$ \mathcal{H}^{(\text{orb})} = \frac{1}{2m} \left[ \mathbf{p} + \frac{e}{c} \mathbf{A}(\mathbf{r}) \right]^2 - e\phi(\mathbf{r}), \quad (16\text{a}) $$

$$ \mathcal{H}^{(\text{Z})} = \frac{e\hbar}{2mc} \boldsymbol{\sigma} \cdot \mathbf{B}(\mathbf{r}), \quad (16\text{b}) $$

representing, respectively, the spin-independent Hamiltonian associated with the orbital motion (orb), and the Zeeman energy of the spin in the magnetic field (Z). We note $\mathbf{p}$ the

electron momentum and $\boldsymbol{\sigma}$ the vector of Pauli matrices, related with the electron spin angular momentum operator by $\mathbf{S} = (\hbar/2)\boldsymbol{\sigma}$. The weakly-relativistic correction $\Delta\mathcal{H}$ in Eq. (14) writes

$$ \Delta\mathcal{H} = \mathcal{H}^{(\text{k})} + \mathcal{H}^{(\text{so-ame})} + \mathcal{H}^{(\text{D})} + \mathcal{H}^{(\text{r})}, \quad (17) $$

with

$$ \mathcal{H}^{(\text{k})} = -\frac{1}{8m^3c^2} \left\{ \boldsymbol{\sigma} \cdot \left[ \mathbf{p} + \frac{e}{c} \mathbf{A}(\mathbf{r}) \right] \right\}^4, \quad (18\text{a}) $$

$$ \mathcal{H}^{(\text{so-ame})} = \frac{e\hbar}{4m^2c^2} \boldsymbol{\sigma} \cdot \left\{ \mathbf{E}(\mathbf{r}) \times \left[ \mathbf{p} + \frac{e}{c} \mathbf{A}(\mathbf{r}) \right] \right\}, \quad (18\text{b}) $$

$$ \mathcal{H}^{(\text{D})} = \frac{e\hbar^2}{8m^2c^2} \nabla \cdot \mathbf{E}(\mathbf{r}), \quad (18\text{c}) $$

$$ \mathcal{H}^{(\text{r})} = \frac{ie\hbar^2}{8m^2c^2} \boldsymbol{\sigma} \cdot [\nabla \times \mathbf{E}(\mathbf{r})]. \quad (18\text{d}) $$

$\mathcal{H}^{(\text{k})}$ contains the first relativistic correction to the kinetic energy, $\mathcal{H}^{(\text{so-ame})}$ combines the spin-orbit (so) and the angular magneto-electric (ame) [69] couplings, $\mathcal{H}^{(\text{D})}$ stands for the Darwin term responsible for the zitterbewegung effect in atomic physics, while the last term $\mathcal{H}^{(\text{r})}$ with the curl of $\mathbf{E}$ vanishes for conservative potentials.

Our idealized nanoparticle is defined by a confining potential $V(\mathbf{r}) = -e\phi(\mathbf{r})$ acting on noninteracting electrons. We assume an applied magnetic field along the $z$-direction, and thus a magnetic induction $\mathbf{B} = B\hat{e}_z$ that in the symmetric gauge can be represented by $\mathbf{A}(\mathbf{r}) = \frac{1}{2}\mathbf{B} \times \mathbf{r}$. Under the previous hypothesis the different components of the Hamiltonian (14) admit simpler expressions that we describe below. The nonrelativistic spinless term (16a) describing the electron orbital motion can be written as

$$ \mathcal{H}^{(\text{orb})} = \mathcal{H}^{(\text{0})} + \mathcal{H}^{(\text{para})} + \mathcal{H}^{(\text{dia})}, \quad (19) $$

with

$$ \mathcal{H}^{(0)} = \frac{\mathbf{p}^2}{2m} + V(\mathbf{r}), \quad (20\text{a}) $$

$$ \mathcal{H}^{(\text{para})} = \frac{\omega_c}{2} L_z, \quad (20\text{b}) $$

$$ \mathcal{H}^{(\text{dia})} = \frac{m\omega_c^2}{8} (x^2 + y^2). \quad (20\text{c}) $$

Here, $\mathcal{H}^{(0)}$ is the zero-field Hamiltonian describing a spinless particle in the confining box, while the terms $\mathcal{H}^{(\text{para})}$ and $\mathcal{H}^{(\text{dia})}$ represent paramagnetic and diamagnetic contributions, respectively, depending of the B-field through the cyclotron frequency $\omega_c = eB/m^*c$ (as we discussed in Sec. I we will identify $m^*$ and $m$). We denote $\mathbf{L}$ the orbital angular momentum of the electron and $L_z$ its $z$ component.

Denoting $\mu_B = e\hbar/2mc$ the Bohr magneton, it is convenient to group the paramagnetic (20b) and Zeeman (16b) components in a term representing the coupling of the total magnetic moment

$$ \boldsymbol{\mu} = -\frac{\mu_B}{\hbar} (\mathbf{L} + g_0 \mathbf{S}) \quad (21) $$

<footnote n="5">⁵ For a complete account of theory and experiments on the persistent current of normal metal rings, see Ref. [64].</footnote>
to the magnetic induction, leading to

$$ \begin{aligned} \mathcal{H}^{(\mu)} &= \mathcal{H}^{(\text{para})} + \mathcal{H}^{(\text{Z})} \\ &= -\boldsymbol{\mu} \cdot \mathbf{B}. \end{aligned} \qquad (22) $$

In the ideal system that we are describing, we use the g-factor $g_0 = 2$.

The correction (18b) takes the form

$$ \mathcal{H}^{(\text{so-ame})} = \mathcal{H}^{(\text{so})} + \mathcal{H}^{(\text{ame})}, \qquad (23) $$

with

$$ \mathcal{H}^{(\text{so})} = \frac{1}{2m^2c^2} \mathbf{S} \cdot [\nabla V(\mathbf{r}) \times \mathbf{p}], \qquad (24a) $$

$$ \mathcal{H}^{(\text{ame})} = \frac{e}{4m^2c^3} \mathbf{S} \cdot [\nabla V(\mathbf{r}) \times (\mathbf{B} \times \mathbf{r})], \qquad (24b) $$

while the Darwin term (18c) can be written as

$$ \mathcal{H}^{(\text{D})} = \frac{\hbar^2}{8m^2c^2} \nabla^2 V(\mathbf{r}). \qquad (25) $$

Our main interest is the effect of the SOC correction (24a) on the magnetic response of different kinds of nanoparticles. In order to assess its relevance, we need to also examine the role of the other weakly-relativistic corrections.

## IV. SPHERICAL NANOPARTICLES

### A. Nanoparticle modeling

We here consider a spherically-symmetric confinement defined by the potential

$$ V(r) = V_0 \Theta(r - a), \qquad (26) $$

where we adopt the spherical coordinates $(r, \theta, \varphi)$ with $r = |\mathbf{r}|$. We denote $\Theta(\zeta)$ the Heaviside step function, $V_0$ the height of the confining wall, and $a$ the nanoparticle radius. The spherical symmetry of the problem allows one to further simplify the expressions of Eq. (24) for the spin-orbit and angular magneto-electric couplings as

$$ \mathcal{H}^{(\text{so})} = \frac{1}{2m^2c^2} \frac{1}{r} \left( \frac{dV}{dr} \right) \mathbf{S} \cdot \mathbf{L}, \qquad (27a) $$

$$ \mathcal{H}^{(\text{ame})} = -\frac{\mu_B B}{2mc^2\hbar} r \left(\frac{dV}{dr}\right) \sin(\theta) \mathbf{S} \cdot \hat{\mathbf{e}}_{\theta}, \qquad (27b) $$

together with that of Eq. (25) for the Darwin term, which can be written as

$$ \mathcal{H}^{(\text{D})} = \frac{\hbar^2}{8m^2c^2} \frac{1}{r^2} \frac{d}{dr} \left( r^2 \frac{dV}{dr} \right). \qquad (28) $$

The previously exposed model of free electrons contained in a hard-wall box constitutes a first approximation for describing noble-metal nanoparticles. Within this scheme, the jellium approximation is adopted for the description of the

conduction electrons and $V(r)$ in Eq. (26) should be understood as the self-consistent electrostatic potential resulting from the confinement due to the positive background from the ionic lattice, together with the mean-field treatment of the electron-electron interaction. Density functional theory calculations [70, 71] indicate that, for not too small nanoparticles and in the absence of a magnetic field, the spherical jellium model (26) is a good approximation for the confining potential, where $V_0 = E_F + W$, with $W$ the work function of the nanoparticle.

The advantages and disadvantages of the jellium model have been thoroughly discussed in the literature concerning the electronic properties of metallic nanoparticles and the physics of surface plasmon resonances [72, 73]. For very small metallic nanoparticles (or clusters) the corrections to the ionic structure have been approached by the use of pseudopotential perturbation theory as a multipole expansion [74] and the effect of surface irregularities and the underlying crystalline lattice has been addressed with numerical methods [75, 76]. The influence of smooth disorder in low-dimensional systems has been addressed with the help of semiclassical methods [77, 78]. We will leave aside the case of very small nanoparticles while neglecting the effect of surface and/or bulk disorder, and we will then not be concerned with such corrections.

The finite height $V_0$ of the confining wall is responsible for the spill-out effect describing the nonzero probability to find an electron density outside the nanoparticle [71]. The proper description of the spill-out, as well as the lifetime of the surface plasmon resonance necessitates to go beyond the discontinuous form (26) of the confining potential, including the abruptness of the potential jump at the nanoparticle surface. In the case of the ZFS, the abruptness of the potential jump is not a crucial parameter, and moreover, the careful analysis of the wave-function behavior close to a potential discontinuity, that we perform in the sequel, shows that the precise value of $V_0$ is not determinant for the nonrelativistic ZFS, nor for the weakly-relativistic corrections (with the exception of the Darwin contribution discussed in Appendix C), justifying the limit of $V_0 \to \infty$ commonly adopted in the nonrelativistic case [29, 50, 51, 79].

When going from the description of the ideal model system to that of a gold nanoparticle, we should in principle use $g_{Au}$ instead of $g_0$ in the expression (21) for the total magnetic moment. However, the small difference between $g_{Au}$ and $g_0$, like the one between $m^*$ and $m$, only induce very small corrections. Another effect to take into account when treating metallic nanoparticles is the Larmor diamagnetic response of the core electrons, which are not considered in the idealized model. While the effect of the core electrons on the ZFS is not altered by the confinement, its contribution must be confronted with the corrections under study in order to gauge the relevance of the latter.
## B. Nonrelativistic susceptibility

The magnetization and the ZFS in the nonrelativistic case have been calculated in the case of a spherical geometry by quantum perturbation theory [50, 51], as well as through the use of semiclassical expansions [29] for the DOS. We show in Appendix A the connection between the two approaches. The quantum procedure starts from the eigenstates of the unperturbed ($B=0$) Hamiltonian $\mathcal{H}^{(0)}$ of Eq. (20a), describing nonrelativistic electrons in a spherical potential (of radius $a$), which are characterized in the product basis by the set $\{\lambda\} = \{n,l,m_z,m_s\}$ of quantum numbers, with $n > 0$ the principal quantum number, $l \ge 0$ the azimuthal quantum number, $m_z \in [-l, l]$ the magnetic quantum number, and $m_s = \pm 1/2$ associated with the spin component along the $z$ direction. The corresponding wave functions are given by the two-component spinors

$$ \Psi_{n,l,m_z, \pm \frac{1}{2}}^{(0)}(\mathbf{r}) = \psi_{n,l,m_z}^{(0)}(\mathbf{r}) \begin{bmatrix} 1 \\ 0 \end{bmatrix}, \quad (29a) $$

$$ \Psi_{n,l,m_z, -\frac{1}{2}}^{(0)}(\mathbf{r}) = \psi_{n,l,m_z}^{(0)}(\mathbf{r}) \begin{bmatrix} 0 \\ 1 \end{bmatrix}, \quad (29b) $$

with the orbital wave function

$$ \psi_{n,l,m_z}^{(0)}(\mathbf{r}) = R_{n,l}(r) Y_l^{m_z}(\vartheta). \quad (30) $$

We note $Y_l^{m_z}(\vartheta)$ the spherical harmonic of degree $l$ and order $m_z$ as a function of the solid angle $\vartheta = (\theta, \varphi)$, while $R_{n,l}(r)$ stands for the associated radial wave function and $E_{n,l} = E_0(k_{n,l} a)^2$ is the corresponding eigenenergy, with $E_0 = \hbar^2/2ma^2$.

For the confining potential (26), $R_{n,l}$ can be expressed in terms of Bessel functions (see Appendix C). In the limiting case of a hard-wall potential ($V_0 \to \infty$), according to Eq. (C6a), the eigenenergy of the $m_z$ and spin-degenerate states characterized by the quantum numbers $l$ and $n$ is

$$ E_{n,l}^{(0)} = E_0 \zeta_{n,l}^2, \quad (31) $$

where $\zeta_{n,l}$ is the $n$th root of the spherical Bessel function of the first kind $j_l(\zeta)$, while

$$ R_{n,l}(r) = \sqrt{\frac{2}{a^3} \frac{j_l(\zeta_{n,l} r/a)}{|j_{l+1}(\zeta_{n,l})|}}. \quad (32) $$

The action of a magnetic field through the term $\mathcal{H}^{(\mu)}$ of Eq. (22) is, one the one hand, to break the above-mentioned degeneracy without giving rise to a new basis of eigenstates. On the other hand, the states $\Psi_{n,l,m_z,m_s}^{(0)}$ [cf. Eq. (29)] are no longer eigenstates once the diamagnetic term $\mathcal{H}^{(\text{dia})}$ of Eq. (20c) is considered under a finite magnetic field. The matrix elements of $\mathcal{H}^{(\text{dia})}$ in the product basis are given in Appendix F 3. Since they are quadratic in $B$, we can treat $\mathcal{H}^{(\text{dia})}$ by first-order perturbation theory. Thus, up to terms of order $B^2$, the perturbed energies are

$$ E_{n,l,m_z,m_s}^{(\text{nr})} = E_{n,l}^{(0)} + \delta E_{n,l,m_z,m_s}^{(\text{nr})}, \quad (33) $$

where $\delta E_{n,l,m_z,m_s}^{(\text{nr})}$ is the magnetic-field correction, which reads

$$ \delta E_{n,l,m_z,m_s}^{(\text{nr})} = E_{m_z, m_s}^{(\mu)} + E_{n,l,m_z}^{(\text{dia})}. \quad (34) $$

Consistently with Eq. (22), we denote

$$ E_{m_z, m_s}^{(\mu)} = E_{m_z}^{(\text{para})} + E_{m_s}^{(Z)}, \quad (35) $$

with

$$ E_{m_z}^{(\text{para})} = \frac{\hbar\omega_c}{2} m_z, \quad (36a) $$

$$ E_{m_s}^{(Z)} = g_0 \mu_B m_s B. \quad (36b) $$

For the case of a hard-wall confinement [50, 51],

$$ E_{n,l,m_z}^{(\text{dia})} = \frac{m\omega_c^2 a^2}{8} \mathcal{R}_{n,l} Y_l^{m_z}, \quad (37) $$

with

$$ \mathcal{R}_{n,l} = \frac{1}{3} \left[ 1 + \frac{2}{\zeta_{nl}^2} \left( l - \frac{1}{2} \right) \left( l + \frac{3}{2} \right) \right], \quad (38a) $$

$$ Y_l^{m_z} = \frac{1}{2} \left[ \frac{m_z^2 + l(l+1) - 1}{(l-1/2)(l+3/2)} \right]. \quad (38b) $$

According to Eq. (6), the nonrelativistic ZFS is determined from the following parameters directly obtained from Eqs. (36) and (37),

$$ E_{n,l,m_z,m_s}^{(\text{nr})} \Big|_{B=0} = E_{n,l}^{(0)}, \quad (39a) $$

$$ \frac{\partial E_{n,l,m_z,m_s}^{(\text{nr})}}{\partial B} \Big|_{B=0} = \mu_B (m_z + 2m_s), \quad (39b) $$

$$ \frac{\partial^2 E_{n,l,m_z,m_s}^{(\text{nr})}}{\partial B^2} \Big|_{B=0} = \frac{\mu_B^2}{2E_0} \mathcal{R}_{n,l} Y_l^{m_z}. \quad (39c) $$

Treating separately each of the field-dependent energy contributions (36), the nonrelativistic ZFS following from Eqs. (6) and (39) can be written as

$$ \chi^{(\text{nr})} = \chi^{(Z)} + \chi^{(\text{orb})}, \quad (40) $$

in terms of a spin-dependent susceptibility $\chi^{(Z)}$ and the orbital component $\chi^{(\text{orb})}$. The latter admits the decomposition

$$ \chi^{(\text{orb})} = \chi^{(\text{para})} + \chi^{(\text{dia})}. \quad (41) $$

Performing the $m_z$ and $m_s$ sums in Eq. (6), we have

$$ \chi^{(Z)} = - \frac{2\mu_B^2}{V} \sum_{l=0}^{\infty} (2l+1) \sum_{n=1}^{\infty} f'_{\bar{\mu}_0} (E_{nl}^{(0)}), \quad (42) $$

as well as [50, 51]

$$ \frac{\chi^{(\text{para})}}{|\chi_{\text{L}}|} = - \frac{3\pi E_0}{k_F a} \sum_{l=0}^{\infty} l(l+1)(2l+1) \sum_{n=1}^{\infty} f'_{\bar{\mu}_0} (E_{nl}^{(0)}), \quad (43a) $$

$$ \frac{\chi^{(\text{dia})}}{|\chi_{\text{L}}|} = - \frac{3\pi}{k_F a} \sum_{l=0}^{\infty} (2l+1) \sum_{n=1}^{\infty} \mathcal{R}_{n,l} f_{\bar{\mu}_0} (E_{nl}^{(0)}). \quad (43b) $$
The mean chemical potential $\bar{\mu}_0$ is associated with the spectrum of the Hamiltonian $\mathcal{H}^{(0)}$ presented in Eq. (20a). Trading the sum over the principal quantum number $n$ by an integral over energy allows to recast Eq. (42) as

$$ \frac{\chi^{(Z)}}{|\chi_L|} = -\frac{9\pi E_0}{k_F a} \sum_{l=0}^{\infty} (2l+1) \int_0^{\infty} dE \, \varrho_l(E) f'_{\mu_0}(E), \quad (44) $$

where we have introduced the $l$-fixed density of states $\varrho_l(E)$ corresponding to the radial problem, itself related to the zero-field DOS by

$$ \varrho(E, 0) = 2 \sum_{l=0}^{l_{\max}} (2l+1) \varrho_l(E), \quad (45) $$

with $l_{\max}$ given by Eq. (A6). The prefactor of 2 in the equation above takes into account the spin degeneracy, as we follow the standard convention of using a spinless one-dimensional $\varrho_l(E)$ and a spinful three-dimensional $\varrho(E, B)$.

For a degenerate electron gas, where $E_F \gg k_B T$, we use $f'_\mu(E) = -\delta(E - \mu)$, and thus Eq. (44) leads to the standard result of the spin-dependent susceptibility [80]

$$ \chi^{(Z)} = \frac{\mu_B^2}{V} \varrho(E_F, 0). \quad (46) $$

In the unconstrained case of $a \to \infty$, the use of the DOS per unit volume for the zero-field, three-dimensional, free electron gas $g^{(3D)}(E) = (m/\pi^2\hbar^2)\sqrt{2mE/\hbar^2}$ in Eq. (46), results in the form (2) of the Pauli susceptibility. In the constrained case of a finite $a$, the separation (7) and the fact that $\bar{\varrho}(E, 0) = g^{(3D)}(E)V$ result in a mean spin-dependent susceptibility $\bar{\chi}^{(Z)} = \chi_P$. Thus, the confinement only adds a small contribution to the bulk susceptibility, which is associated with $\varrho^{\text{osc}}(E, 0)$ [81].

In the constrained case, the numerical implementation of Eqs. (42) and (43) leads to the nonrelativistic ZFS $\chi^{(nr)}$ presented in Fig. 1 (in black) as a function of $k_F a$. In the shown

interval, these numerical results are almost indistinguishable from those in which $\chi_b^{(nr)} = 2|\chi_L|$ is added to the semiclassical ZFS $\chi^{(orb)-osc}$ of Eq. (A19). The oscillations of $\chi^{(nr)}$ as a function of $k_F a$ are typically much larger than $|\chi_L|$ and can be understood as a shell structure [82]. The suppression of these oscillations for large sizes can be understood, at the semiclassical level, by the thermal damping (13) acting on the contribution of each family of classical periodic orbits. For very large $a$ (inset), the oscillations of $\chi^{(nr)}$ are quite reduced, and we can see that they are around the bulk ZFS $\chi_b^{(nr)}$, given by Eqs. (1) and (2), as the confinement becomes irrelevant in such a limit. The typical values of $\chi^{(nr)}$ are important in order to assess the relevance of the relativistic corrections to be calculated in the sequel. Also for comparison purposes, we show in Fig. 1 (in red) the ZFS $\chi^{(2)} + \chi_b^{(nr)}$ arising from finite-$N$ corrections to the free energy (10c) using the semiclassical expression for the oscillating part of the density of states $\varrho^{\text{osc}}(E, B)$ given by Eqs. (A7), (A8), and (A11), as detailed in Ref. [29].<fnref n="6" /> In addition, we show in Fig. 1 (blue line) the ZFS $\chi_{\text{ens}}^d + \chi_b^{(nr)}$ of an ensemble of metallic nanoparticles with an important size dispersion. Here, $\chi_{\text{ens}}^d$ is the ensemble average of $\chi^{(2)}$ over a Gaussian probability distribution of the size parameter $a$ (cf. Eq. (30) in Ref. [29]).

### C. Kinetic correction

The zero-field Dirac equation in a spherical potential box admits an exact solution [83, 84]. While the inclusion of an infinitesimal magnetic field allowing to address the ZFS can in principle be implemented as a perturbation [85], it is simpler to proceed from the weakly-relativistic Hamiltonian (14). The kinetic term (18a) can be written as $\mathcal{H}^{(k)} = -\{\mathcal{H}^{(nr)} + e\phi\}^2/2mc^2$. Therefore, in the case of a hard-wall confinement, and up to quadratic terms in $B$, it leads to the energy correction

$$ \begin{aligned} E_{n,l,m_z,m_s}^{(k)} &= -\frac{1}{2mc^2} \langle \Psi_{n,l,m_z,m_s}^{(0)} \left| \left(\frac{\mathbf{p}^2}{2m}\right)^2 + \frac{\mathbf{p}^2}{m} \mathcal{H}^{(\mu)} + (\mathcal{H}^{(\mu)})^2 + \frac{\mathbf{p}^2}{2m} \mathcal{H}^{(\text{dia})} + \mathcal{H}^{(\text{dia})} \frac{\mathbf{p}^2}{2m} |\Psi_{n,l,m_z,m_s}^{(0)} \rangle \\ &= -\frac{(E_{n,l}^{(0)})^2}{2mc^2} + \delta E_{n,l,m_z,m_s}^{(k)}, \end{aligned} \quad (47) $$

where the $B$-dependent component is given by

$$ \delta E_{n,l,m_z,m_s}^{(k)} = -\frac{1}{2mc^2} \left[ 2E_{n,l}^{(0)} E_{m_z,m_s}^{(\mu)} + (E_{m_z,m_s}^{(\mu)})^2 + 2E_{n,l}^{(0)} E_{n,l,m_z}^{(\text{dia})} \right], \quad (48) $$

with $E_{n,l}^{(0)}$, $E_{m_z,m_s}^{(\mu)}$, and $E_{n,l,m_z}^{(\text{dia})}$ given, respectively, by Eqs. (31), (35), and (37). Since the product basis is constituted of eigenvectors of $\mathbf{p}^4$ and $\mathcal{H}^{(\mu)}$, the off-diagonal matrix elements of $\mathcal{H}^{(k)}$ are of quadratic order in $B$, and thus do not need to be considered.

<footnote n="6">⁶ Up to the constant term $\chi_b^{(nr)}$, Fig. 1 presents equivalent results to those of</footnote>

<footnote>Fig. 9(b) of Ref. [29], correcting a factor of $\pi/4$ that erroneously overrated the contribution $\chi^{(2)}$.</footnote>
<figure><img src="image_2.png" /><figcaption>FIG. 1. In black: nonrelativistic ZFS $\chi^{(nr)}$ for $k_B T/\bar{\mu}_0 = 5 \times 10^{-3}$ (corresponding to room temperature in the case of gold) obtained from Eqs. (42)–(43) (in units of the absolute value of the Landau susceptibility $\chi_L$), as a function of the nanoparticle radius $a$ (scaled with the Fermi wave vector $k_F$). In red: $\chi^{(2)} + \chi_b^{(nr)}$, with $\chi^{(2)}$ obtained from the finite-$N$ correction (10c) to the free energy and the oscillating part of the $B$-dependent density of states $g^{\text{osc}}(E, B)$. In blue: ZFS of an ensemble of nanoparticles with an important size dispersion, $\chi_{\text{ens}}^d + \chi_b^{(nr)}$, where $\chi_{\text{ens}}^d$ is obtained by taking the average of $\chi^{(2)}$ over a Gaussian probability distribution (in size). Inset: Corresponding ZFSs for larger nanoparticle radii, showing the approach to the bulk value given by Eqs. (1) and (2) indicated by a dashed line.</figcaption></figure>

Even if the modification of the nonrelativistic ZFS $\chi^{(nr)}$ due to the kinetic correction (47) does not have a physical meaning by itself, it is nevertheless interesting to calculate it in view of weighting its importance against the other weakly-relativistic modifications. Moreover, in the case of a very large radius $a$, where the role of the confining potential $V(r)$ should become irrelevant, the ZFS $\chi^{(nr-k)}$ taking only into account the modification of $\chi^{(nr)}$ due to the correction (47), can be compared with the bulk weakly-relativistic ZFS $\chi_b^{(wr)}$ of Eq. (3).

An important aspect of the correction (47) is that it induces at $B = 0$ an energy shift of all levels. In particular, within the grand canonical ensemble, the shift of the Fermi level translates into a renormalization

$$ \Delta\mu^{(k)} \approx -\frac{\bar{\mu}_0^2}{2mc^2} \qquad (49) $$

of the zero-field nonrelativistic mean chemical potential $\bar{\mu}_0$.

The kinetic correction (47) results in the eigenenergies $E_{n,l,m_z,m_s}^{(nr-k)} = E_{n,l,m_s,m_s}^{(nr)} + E_{n,l,m_z,m_s}^{(k)}$ from which the ZFS $\chi^{(nr-k)}$ can be obtained by using the parameters

$$ E_{n,l,m_z,m_s}^{(nr-k)} \Big|_{B=0} = E_{n,l}^{(0)} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right], \qquad (50a) $$

$$ \frac{\partial E_{n,l,m_z,m_s}^{(nr-k)}}{\partial B} \Big|_{B=0} = \mu_B (m_z + 2m_s) \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right], \qquad (50b) $$

$$ \frac{\partial^2 E_{n,l,m_z,m_s}^{(nr-k)}}{\partial B^2} \Big|_{B=0} = \mu_B^2 \left( \frac{1}{2E_0} \mathcal{R}_{n,l} \mathcal{Y}_l^{m_z} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] - \frac{1}{mc^2} (m_z + 2m_s)^2 \right). \qquad (50c) $$
Performing the $m_z$ and $m_s$ sums in Eq. (6), while working up to linear order in $E_{n,l}^{(0)}/mc^2$, we have

$$ \frac{\chi^{(\text{nr-k})}}{|\chi_L|} = - \frac{6\pi E_0}{k_F a} \sum_{l=0}^{\infty} (l+1/2) \sum_{n=0}^{\infty} \left\{ \begin{aligned} & \left[ l(l+1) + 3 \right] f'_{\bar{\mu}} \left( E_{n,l}^{(0)} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right] \right) \left[ 1 - 2 \frac{E_{n,l}^{(0)}}{mc^2} \right] \\ & + \frac{\mathcal{R}_{n,l}}{E_0} f_{\bar{\mu}} \left( E_{n,l}^{(0)} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right] \right) \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] \\ & - \frac{l(l+1)+3}{mc^2} f_{\bar{\mu}} \left( E_{n,l}^{(0)} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right] \right) \end{aligned} \right\} . \quad (51) $$

As discussed above, we are interested in the correction

$$ \Delta\chi^{(k)} = \chi^{(\text{nr-k})} - \chi^{(\text{nr})}, \quad (52) $$

where $\chi^{(\text{nr})}$ is defined in Eq. (40), and thus associated with the mean chemical potential $\bar{\mu}_0$ of the nonrelativistic problem, while $\chi^{(\text{nr-k})}$ is associated with the renormalized mean chemical potential $\bar{\mu} = \bar{\mu}_0 + \Delta\mu^{(k)}$ [cf. Eq. (49)]. We then write

$$ \frac{\Delta\chi^{(k)}}{|\chi_L|} = - \frac{6\pi E_0}{k_F a mc^2} \sum_{l=0}^{\infty} (l+1/2) \sum_{n=1}^{\infty} \left\{ \begin{aligned} & \left[ l(l+1) + 3 \right] f''_{\bar{\mu}_0} (E_{n,l}^{(0)}) + \frac{\mathcal{R}_{n,l}}{E_0} f'_{\bar{\mu}_0} (E_{n,l}^{(0)}) \right\} \left[ \frac{\left(E_{n,l}^{(0)}\right)^2}{2mc^2} - \Delta\mu^{(k)} \right] \\ & + \frac{6\pi E_0}{k_F a mc^2} \sum_{l=0}^{\infty} (l+1/2) \sum_{n=1}^{\infty} \left\{ \begin{aligned} & 2[l(l+1)+3] E_{n,l}^{(0)} f'_{\bar{\mu}_0} (E_{n,l}^{(0)}) + \left( \frac{\mathcal{R}_{n,l}}{E_0} E_{n,l}^{(0)} + l(l+1)+3 \right) f_{\bar{\mu}_0} (E_{n,l}^{(0)}) \end{aligned} \right\} . \quad (53) $$

The first sum in the above equation results from the zero-field component of the correction $E_{n,l,m_z,m_s}^{(k)}$, and thus it could have alternatively been derived by simply implementing in the ZFS expressions (42)–(43) the shift of $E_{n,l}^{(0)}$ to $E_{n,l}^{(0)}[1 - E_{n,l}^{(0)}/2mc^2]$, with the corresponding renormalization $\Delta\mu^{(k)}$ of the chemical potential. The second sum in (53) represents the nontrivial effect of the magnetic-field dependent contribution $\delta E_{n,l,m_z,m_s}^{(k)}$ to the kinetic correction, given by Eq. (48).

In Fig. 2 we present the numerical evaluation of the kinetic correction $\Delta\chi^{(k)}$ to the nonrelativistic ZFS, according to Eq. (53), as a function of $k_F a$. For large $a$ (inset), the values of $\Delta\chi^{(k)}$ oscillate around $-1/3|\chi_L| (v_F/c)^2$, consistently with the result (3) for the bulk ZFS $\chi_b^{(\text{wr})}$, since in such a limit the confinement becomes irrelevant and the kinetic correction is the only weakly-relativistic effect that needs to be taken into account. The oscillations as a function of $k_F a$ are much larger than the bulk value, but remain considerably smaller than the typical values of $\chi^{(\text{nr})}$ exhibited in Fig. 1 (and even of the finite-$N$ correction characterizing the response of an ensemble of nanoparticles, see the blue line in Fig. 1).

The correction $\Delta\chi^{(k)}$ of Eq. (53) admits a one-dimensional semiclassical treatment, analogous to that of the nonrelativistic ZFS treated in Appendix A. In particular, the smooth part of $\Delta\chi^{(k)}$ can be evaluated along the lines of Eqs. (A14)–(A16), and we find in the leading order in $k_F a$, that the smooth part of the first sum behaves as $1/5 (k_F a)^2 (v_F/c)^2$, while the second sum cancels the previous contribution. The next-leading term of the smooth part of $\Delta\chi^{(k)}$ is of order $(k_F a)^0$, but for the same reasons discussed in Appendix A, the correct asymptotic value $\chi_b^{(\text{wr})}$ of the bulk yielded by the numerical calculation, is not recovered by our semiclassical approach.

### D. Spin-orbit coupling in a spherical nanoparticle

While the product eigenbasis of $\mathcal{H}^{(0)}$ used in Sec. IV B, and characterized by the set $\{\lambda\} = \{n, l, m_z, m_s\}$, remained the appropriate one once the term $\mathcal{H}^{(\mu)}$ of Eq. (22) was taken into account, the inclusion of the SOC term $\mathcal{H}^{(\text{so})}$ of Eq. (27a) makes it more convenient to change to the eigenbasis of the total angular momentum $\mathbf{J} = \mathbf{L} + \mathbf{S}$, characterized by the set $\{\tilde{\lambda}\} = \{n, j, m_j, l\}$. The rule of addition of angular momentum results in $[l] \otimes [1/2] = [l - 1/2, l] \oplus [l + 1/2, l]$, where the left-hand side is the tensor product of the orbital angular momentum $l \neq 0$ subspace with that of the spin $1/2$ and the right-hand side represents the direct sum of the subspaces with total angular momenta $j = l \mp 1/2$ and orbital angular momentum $l$. For radially symmetric electrostatic potentials the separability between radial and angular coordinates allows to write the eigenstates of the coupled basis in terms of the spinors [83]

$$ \Phi_{n,j,m_j}^{(\pm)}(\mathbf{r}) = R_{n,j \pm 1/2}(r) \Upsilon_{j,m_j}^{(\pm)}(\vartheta), \quad (54) $$

where we have defined the spinor spherical harmonics

$$ \Upsilon_{j,m_j}^{(\pm)}(\vartheta) = \frac{1}{\sqrt{2(j \pm 1/2) + 1}} \times \left[ \frac{\mp \sqrt{j \pm 1/2 \mp m_j + 1/2} Y_{j \pm 1/2}^{m_j - 1/2}(\vartheta)}{\sqrt{j \pm 1/2 \pm m_j + 1/2} Y_{j \pm 1/2}^{m_j + 1/2}(\vartheta)} \right] \quad (55) $$

in terms of the usual spherical harmonics $Y_l^{m_z}$ introduced in Eq. (30). The label $(\pm)$ corresponds to $l = j \pm 1/2$, and sets the parity $(-1)^{j \pm 1/2}$ of the state. The associated eigenenergies of $\mathcal{H}^{(0)}$ in the case of a spherical confining potential
<figure><img src="image_3.png" /><figcaption>FIG. 2. Kinetic correction $\Delta\chi^{(k)}$ to the ZFS for $k_B T/\bar{\mu}_0 = 5 \times 10^{-3}$ obtained from Eq. (53) [in units of the absolute value of the Landau susceptibility $\chi_L$ and multiplied by the scaling factor $(c/v_F)^2$], as a function of the nanoparticle radius $a$ (scaled with the Fermi wave vector $k_F$). Inset: $\Delta\chi^{(k)}$ (in the same above-mentioned units as in the main panel) corresponding to larger sizes, showing the approach to the weakly-relativistic bulk value $-1/3$ [cf. Eq. (3)] and indicated by a dashed line.</figcaption></figure>

are $\mathcal{E}_{n,j,(\pm)}^{(0)} = E_{n,j\pm 1/2}^{(0)}$, and therefore there is degeneracy between the $2j+1$ dimensional subspace $\{n, j, (+)\}$ and the $2j+3$ dimensional subspace $\{n, j+1, (-)\}$, as they are characterized by the same quantum number $l$.

The coupled basis remains an eigenbasis of the subspace $\{n, j, (+)\} \oplus \{n, j+1, (-)\}$ once $\mathcal{H}^{(\text{so})}$ is taken into account, while the eigenenergies change according to

$$ 
\begin{aligned}
\mathcal{E}_{n,j,(\pm)}^{(\text{so})} &= \langle \Phi_{n,j,m_j}^{(\pm)} | \mathcal{H}^{(\text{so})} | \Phi_{n,j,m_j}^{(\pm)} \rangle \\
&= \frac{\hbar^2}{4m^2c^2a^2} \left[ \mp \left(j + \frac{1}{2}\right) - 1 \right] I_{n,j\pm 1/2}^{(\text{so})}, \quad (56)
\end{aligned}
$$

with the radial matrix element

$$ I_{n,l}^{(\text{so})} = a^2 \int_0^\infty dr \, r [R_{n,l}(r)]^2 V'(r). \quad (57) $$

For the potential (26), we have

$$ I_{n,l}^{(\text{so})} = V_0 a^3 [R_{n,l}(a)]^2. \quad (58) $$

In the limit where the confining potential approaches a hard wall, the product $V_0 [R_{n,l}(a)]^2$ remains finite [86]. Using the limiting expressions (C6) and (C7c) we obtain

$$ I_{n,l}^{(\text{so})} = \frac{\hbar^2}{ma^2} \zeta_{n,l}^2, \quad (59) $$

and therefore

$$ \mathcal{E}_{n,j,(\pm)}^{(\text{so})} = \left[ \mp \left(j + \frac{1}{2}\right) - 1 \right] \zeta_{n,j\pm 1/2}^2 \frac{E_0^2}{mc^2}, \quad (60) $$

independently of $V_0$.

Obviously, $\mathcal{E}_{n,1/2,(-)}^{(\text{so})} = 0$, since the spherical symmetry of the s-states ($l=0$) renders the SOC ineffective. The degeneracy between the subspaces $\{n, j, (+)\}$ and $\{n, j+1, (-)\}$ for $B=0$ is broken by $\mathcal{H}^{(\text{so})}$. For $B \neq 0$, the remaining degeneracy within each subspace is lifted according to the different values of $m_j$.

### E. Perturbative treatment of the magnetic field

Once the term $\mathcal{H}^{(\mu)}$ of Eq. (22) is taken into account, $m_j$ is still a good quantum number, but the coupled basis is no longer an eigenbasis of the subspace $\{n, j, (+)\} \oplus \{n, j+1, (-)\}$. Therefore, in order to treat the terms $\mathcal{H}^{(\mu)}$, $\mathcal{H}^{(\text{dia})}$, $\mathcal{H}^{(k)}$, $\mathcal{H}^{(\text{ame})}$ in the two lowest orders in $B$, the perturbative approach in magnetic field of Sec. IV B yielding the nonrelativistic ZFS has to be extended using the decomposition in
subspaces of fixed $m_j$. These are represented by

$$ \{n, j, (+)\} \oplus \{n, j + 1, (-)\} = S_{n,j+1/2,(d)}^{\mathrm{P}} \\ \oplus \bigcup_{m_j = -j}^{m_j = +j} S_{n,j+1/2,m_j}^{\mathrm{e}} \oplus S_{n,j+1/2,(u)}^{\mathrm{P}} \quad (61) $$

in terms of the down (d) and up (u) one-dimensional subspaces

$$ S_{n,l,(d/u)}^{\mathrm{P}} = \{n, l + 1/2, (-)\} \big|_{m_j = \mp(l+1/2)}, \quad (62) $$

and the two-dimensional subspaces $S_{n,l,m_j}^e$ subtended by the vectors $|\Phi_{n,l \mp 1/2, m_j}^{(\pm)}\rangle$ of the coupled basis, associated with the quantum numbers $\{n, l \mp 1/2, m_j, l\}$ with $l \neq 0$ and $|m_j| \le l - 1/2$. The labels p (e) stand for "product" ("entangled"), characterizing the one (two)-dimensional subspaces where the coupled basis does (does not) coincide with the product basis. The choice of using the index $l$ (instead of $j$) in order to label the subspaces is motivated in view of the book-keeping for the sum over states. We notice that the definition (62) is also valid for the subspaces with $l = 0$, which are not considered in the decomposition (61), but should be included when taking the sum over states yielding a thermodynamic quantity like the ZFS.

The Hamiltonian (14), restricted to the subspaces $S_{n,l,m_j}^e$, can be expressed by the $2 \times 2$ matrix

$$ \mathcal{H}_{n,l,m_j} = \begin{pmatrix} E_{n,l,m_j}^{(+)} & -\mu_B B \mathfrak{N}_{l,m_j} \\ -\mu_B B \mathfrak{N}_{l,m_j} & E_{m,l,m_j}^{(-)} \end{pmatrix} . \quad (63) $$

In the diagonal matrix elements $E_{n,l,m_j}^{(\pm)}$ we separate the field-independent and the field-dependent contributions as

$$ E_{n,l,m_j}^{(\pm)} = E_{n,l,(\pm)}^{(0)} + \delta E_{n,l,m_j}^{(\pm)}, \quad (64) $$

with

$$ E_{n,l,(\pm)}^{(0)} = E_{n,l}^{(0)} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right] + E_{n,l}^{(\mathrm{D})} + E_{n,l,(\pm)}^{(\mathrm{so})} . \quad (65) $$

The first term in the right-hand side of the above equation represents the $B=0$ nonrelativistic eigenvalue together with its kinetic energy correction, as expressed in Eq. (50a). The second term is the Darwin correction, given for the case of $V_0 \gg E_0$ by Eq. (C13). Recalling Eq. (60), it is convenient to express the SOC contribution as

$$ \begin{aligned} E_{n,l,(\pm)}^{(\mathrm{so})} &= \mathcal{E}_{n,l \mp 1/2,(\pm)}^{(\mathrm{so})} \\ &= \left[ \mp \left( l + \frac{1}{2} \right) - \frac{1}{2} \right] \zeta_{n,l}^2 \frac{E_0^2}{mc^2} . \end{aligned} \quad (66) $$

The spin splitting in the subspace $S_{n,l,m_j}^e$ is

$$ \begin{aligned} \Delta_{n,l} &= E_{n,l,(-)}^{(0)} - E_{n,l,(+)}^{(0)} \\ &= E_{n,l,(-)}^{(\mathrm{so})} - E_{n,l,(+)}^{(\mathrm{so})} \\ &= 2 \left( l + \frac{1}{2} \right) \zeta_{n,l}^2 \frac{E_0^2}{mc^2} . \end{aligned} \quad (67) $$

In reducing the Hamiltonian (14) to its $2 \times 2$ form (63), we are assuming an infinitesimal field $B$, and neglecting the coupling between the subspaces $S_{n,l,m_j}^e$ and $S_{n',l,m_j}^e$. This last approximation requires $\Delta_{n,l} \ll E_{n+1,l}^{(0)} - E_{n,l}^{(0)} \simeq 2\pi\zeta_{n,l}E_0$, where we have used the asymptotic form of the zeros of $j_l(\zeta)$. The previous condition translates into $(l+1/2)(v_F/c)^2 \ll 2\pi k_F a$, which is always verified in the weakly-relativistic limit since $l_{max} + 1/2 \simeq k_F a$ (see Appendix A). Nevertheless, our perturbative treatment of $\mathcal{H}^{(\mathrm{so})}$ can become problematic in cases of quasi-degeneracies between eigenstates $E_{n,l}^{(0)}$ and $E_{n',l'}^{(0)}$ corresponding to different quantum numbers. We remark that $\mathcal{H}^{(\mathrm{so})}$ does not induce a renormalization of the chemical potential, since the shift $-E_0 E_{n,l}^{(0)}/mc^2$ associated with each of the subspaces $S_{n,l,m_j}^e$ is compensated by the correction $lE_0 E_{n,l}^{(0)}/mc^2$ characterizing each of the two subspaces $S_{n,l,(d/u)}^P$.

The $B$-dependent components of the diagonal matrix elements (64) are given by

$$ \begin{aligned} \delta E_{n,l,m_j}^{(\pm)} ={}& E_{l,m_j,(\pm)}^{(\mu)} + E_{n,l,m_j,(\pm)}^{(\mathrm{dia})} \\ &+ \delta E_{n,l,m_j,(\pm)}^{(\mathrm{k})} + E_{n,l,m_j,(\pm)}^{(\mathrm{ame})} . \end{aligned} \quad (68) $$

The first term can be obtained with the aid of the Wigner-Eckart theorem, which allows us to write

$$ \begin{aligned} \mathcal{E}_{j,m_j,(\pm)}^{(\mu)} &= \langle \Phi_{n,j,m_j}^{(\pm)} | \mathcal{H}^{(\mu)} | \Phi_{n,j,m_j}^{(\pm)} \rangle \\ &= \frac{j + 1/2}{j + 1/2 \pm 1/2} \mu_B B m_j, \end{aligned} \quad (69) $$

independently of $n$, and therefore

$$ \begin{aligned} E_{l,m_j,(\pm)}^{(\mu)} &= \mathcal{E}_{l \mp 1/2, m_j, (\pm)}^{(\mu)} \\ &= \left( 1 \mp \frac{1/2}{l + 1/2} \right) \mu_B B m_j . \end{aligned} \quad (70) $$

The diamagnetic contributions in Eq. (68) follow from the matrix elements

$$ \begin{aligned} \mathcal{E}_{n,j,m_j,(\pm)}^{(\mathrm{dia})} &= \langle \Phi_{n,j,m_j}^{(\pm)} | \mathcal{H}^{(\mathrm{dia})} | \Phi_{n,j,m_j}^{(\pm)} \rangle \\ &= \left[ 1 + \frac{m_j^2}{j(j+1)} \right] \frac{\mu_B^2 B^2}{8E_0} \mathcal{R}_{n,j \pm 1/2}, \end{aligned} \quad (71) $$

where $\mathcal{R}_{n,l}$ has been defined in Eq. (38a). Thus,

$$ \begin{aligned} E_{n,l,m_j,(\pm)}^{(\mathrm{dia})} &= \mathcal{E}_{n,l \mp 1/2, m_j, (\pm)}^{(\mathrm{dia})} \\ &= \left[ 1 + \frac{m_j^2}{(l \mp 1/2)(l \mp 1/2 + 1)} \right] \frac{\mu_B^2 B^2}{8E_0} \mathcal{R}_{n,l}. \end{aligned} \quad (72) $$

The diagonal matrix elements of $\mathcal{H}^{(\mathrm{k})}$ in the subspace $S_{n,l,m_j}^e$ can be obtained, up to quadratic order in $B$, by a similar expression to that of Eq. (47), where the vectors
$|\Psi_{n,l,m_z,m_s}^{(0)}\rangle$ of the product basis have to be replaced by the ones of the coupled basis, $|\Phi_{n,j,m_j}^{(\pm)}\rangle$, leading to

$$E_{n,l,m_j,(\pm)}^{(k)} = -\frac{\left(E_{n,l}^{(0)}\right)^2}{2mc^2} + \delta E_{n,l,m_j,(\pm)}^{(k)}, \quad (73)$$

with

$$\delta E_{n,l,m_j,(\pm)}^{(k)} = -\frac{E_{n,l}^{(0)}}{mc^2} \left(E_{l,m_j,(\pm)}^{(\mu)} + E_{n,l,m_j,(\pm)}^{(\text{dia})}\right) - \frac{(\mu_B B)^2}{2mc^2} \left[m_j^2 \left(1 \mp \frac{1}{l+1/2}\right) + \frac{1}{4}\right]. \quad (74)$$

The last contribution to $\delta E_{n,l,m_j}^{(\pm)}$ in Eq. (68) results from the diagonal matrix element of $\mathcal{H}^{(\text{ame})}$ in the coupled basis, which with the help of Eqs. (D1), (59), and (D7), can be written as

$$\begin{aligned} E_{n,l,m_j,(\pm)}^{(\text{ame})} &= \mathcal{E}_{n,l\mp 1/2,m_j,(\pm)}^{(\text{ame})} \\ &= \mp \mu_B B \frac{E_{n,l}^{(0)}}{2mc^2} \frac{m_j(l \mp 1/2 + 1/2)}{(l \mp 1/2)(l \mp 1/2 + 1)}. \end{aligned} \quad (75)$$

We remark that the expression (70) of $E_{l,m_j,(\pm)}^{(\mu)}$ does not simply follow from the result (35) for $E_{m_z,m_s}^{(\mu)}$, since the latter represents the exact energy shift associated with $\mathcal{H}^{(\mu)}$ in the product basis, while the former is just the first-order perturbative correction in the coupled basis. Similarly, the perturbative correction (72) for $E_{n,l,m_j,(\pm)}^{(\text{dia})}$ does not simply follow from the analogous correction (37) for $E_{n,l,m_z}^{(\text{dia})}$, nor does $\delta E_{n,l,m_j,(\pm)}^{(k)}$ in Eq. (74) from $E_{n,l,m_z,m_s}^{(k)}$ in Eq. (47). In the same vein, the angular magneto-electric correction $E_{n,l,m_j,(\pm)}^{(\text{ame})}$ of Eq. (75) follows from the matrix element $\mathcal{E}_{n,l\mp 1/2,m_j,(\pm)}^{(\text{ame})}$ in the coupled basis.

The off-diagonal matrix elements of the restricted Hamiltonian (63) need to be only obtained up to linear order in B. Thus, the nonrelativistic component does not have a term associated with $\mathcal{H}^{(\text{dia})}$, but consists only of

$$\langle \Phi_{n,l+1/2,m_j}^{(-)} | \mathcal{H}^{(\mu)} | \Phi_{n,l-1/2,m_j}^{(+)} \rangle = - \mu_B B \frac{\sqrt{(l+1/2)^2 - m_j^2}}{2(l+1/2)}, \quad (76)$$

independently of n. The only off-diagonal matrix element of $\mathcal{H}^{(k)}$ that we need to consider is

$$\langle \Phi_{n,l+1/2,m_j}^{(-)} | \left(-\frac{1}{mc^2}\right) \frac{\mathbf{p}^2}{2m} \mathcal{H}^{(\mu)} | \Phi_{n,l-1/2,m_j}^{(+)} \rangle = \mu_B B \frac{E_{n,l}^{(0)}}{2mc^2} \frac{\sqrt{(l+1/2)^2 - m_j^2}}{l+1/2}. \quad (77)$$

The remaining contribution to the off-diagonal matrix element arising from the angular magneto-electric coupling is given by Eqs. (D8) and (D12) as

$$\langle \Phi_{n,l+1/2,m_j}^{(-)} | \mathcal{H}^{(\text{ame})} | \Phi_{n,l-1/2,m_j}^{(+)} \rangle = - \mu_B B \frac{E_{n,l}^{(0)}}{4mc^2} \frac{\sqrt{(l+1/2)^2 - m_j^2}}{l+1/2}. \quad (78)$$

The off-diagonal matrix element of $\mathcal{H}_{n,l,m_j}$ with the form of Eq. (63) is determined by the B-independent dimensionless parameter

$$\mathcal{N}_{l,m_j} = \frac{\sqrt{(l+1/2)^2 - m_j^2}}{2(l+1/2)} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right]. \quad (79)$$

The diagonalization of $\mathcal{H}_{n,l,m_j}$ yields the B-dependent low (l) and high (h) eigenenergies within the subspaces $S_{n,l,m_j}^e$, respectively given by

$$E_{n,l,m_j}^{(\text{l/h})} = \frac{E_{n,l,m_j}^{(-)} + E_{n,l,m_j}^{(+)}}{2} \mp \sqrt{\left(\frac{E_{n,l,m_j}^{(-)} - E_{n,l,m_j}^{(+)}}{2}\right)^2 + \mu_B^2 B^2 \mathcal{N}_{l,m_j}^2}. \quad (80)$$

In analogy with the nonrelativistic parameters of Eq. (39), the contributions to the ZFS stemming from the (l/h) eigenstates within the subspaces $S_{n,l,m_j}^e$ are determined by

$$E_{n,l,m_j}^{(\text{l/h})} \Big|_{B=0} = E_{n,l,m_j}^{(\pm)} \Big|_{B=0} = E_{n,l,(\pm)}^{(0)}, \quad (81a)$$

$$\left. \frac{\partial E_{n,l,m_j}^{(\text{l/h})}}{\partial B} \right|_{B=0} = \left. \frac{\partial E_{n,l,m_j}^{(\pm)}}{\partial B} \right|_{B=0} = \mu_B m_j \left\{ \left( 1 \mp \frac{1}{l+1/2} \right) \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] \mp \frac{E_{n,l}^{(0)}}{2mc^2} \frac{l \mp 1/2 + 1/2}{(l \mp 1/2)(l \mp 1/2 + 1)} \right\}, \quad (81b)$$
$$ \left. \frac{\partial^2 E_{n,l,m_j}^{(\text{l/h})}}{\partial B^2} \right|_{B=0} = \left. \frac{\partial^2 E_{n,l,m_j}^{(\pm)}}{\partial B^2} \right|_{B=0} \mp \frac{2\mu_B^2 \mathcal{R}_{l,m_j}^2}{\Delta_{n,l}} = \frac{\mu_B^2}{4} \left\{ \left( 1 + \frac{m_j^2}{(l \mp 1/2)(l + 1 \mp 1/2)} \right) \frac{\mathcal{R}_{n,l}}{E_0} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] \right. \\ \left. - \frac{1}{mc^2} \left( 4m_j^2 \left( 1 \mp \frac{1}{l+1/2} \right) + 1 \right) \mp \frac{(l+1/2)^2 - m_j^2}{(E_0/mc^2) E_{n,l}^{(0)} (l+1/2)^3} \left[ 1 - \frac{E_{n,l}^{(0)}}{2mc^2} \right]^2 \right\}. \quad (81c) $$

The eigenenergies of the down (up) state of the subspace $\mathcal{S}_{n,l,(d/u)}^p$ can be written as

$$ E_{n,l}^{(d/u)} = E_{n,l,(-)}^{(0)} + \delta E_{n,l}^{(d/u)}, \quad (82) $$

with

$$ \delta E_{n,l}^{(d/u)} = E_{l, \mp(l+1/2), (-)}^{(\mu)} + E_{n,l, \mp(l+1/2), (-)}^{(\text{dia})} + \delta E_{n,l, \mp(l+1/2), (-)}^{(\text{k})} + E_{n,l, \mp(l+1/2), (-)}^{(\text{ame})}. \quad (83) $$

Therefore, in addition to the parameters (81), we have to consider the contributions to the ZFS associated with the (d/u) states, which follow from

$$ E_{n,l}^{(d/u)} \Big|_{B=0} = E_{n,l,(-)}^{(0)}, \quad (84a) $$

$$ \frac{\partial E_{n,l}^{(d/u)}}{\partial B} \Big|_{B=0} = \mp \mu_B (l+1) \left\{ 1 - \frac{E_{n,l}^{(0)}}{mc^2} + \frac{E_{n,l}^{(0)}}{2mc^2} \frac{1}{l+3/2} \right\}, \quad (84b) $$

$$ \frac{\partial^2 E_{n,l}^{(d/u)}}{\partial B^2} \Big|_{B=0} = \mu_B^2 \left\{ \frac{l+1}{l+3/2} \frac{\mathcal{R}_{n,l}}{2E_0} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] - \frac{1}{mc^2} (l+1)^2 \right\}. \quad (84c) $$

The weakly-relativistic ZFS follows from the evaluation of Eq. (6) using the parameters of Eqs. (81) and (84).

## F. Relativistic correction to the zero-field susceptibility

The weakly-relativistic ZFS will have contributions from p and e eigenstates, and it can therefore be written as $\chi^{(\text{wr})} = \chi_P + \chi_e$. Performing the $m_j$ summation in Eq. (6) for each subset of eigenstates while keeping up to linear terms in $E_{n,l}^{(0)}/mc^2$, we have

$$ \frac{\chi_P}{|\chi_L|} = - \frac{9\pi E_0}{k_F a} \sum_{l=0}^{\infty} \sum_{n=1}^{\infty} \left\{ (l+1)^2 \left[ 1 - \frac{2E_{n,l}^{(0)}}{mc^2} \frac{l+1}{l+3/2} \right] f_{\bar{\mu}}'(E_{n,l,(-)}^{(0)}) \right. \\ \left. + \left( \frac{l+1}{l+3/2} \frac{\mathcal{R}_{n,l}}{2E_0} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] - \frac{1}{mc^2} (l+1)^2 \right) f_{\bar{\mu}}(E_{n,l,(-)}^{(0)}) \right\}, \quad (85a) $$

$$ \frac{\chi_e}{|\chi_L|} = - \frac{3\pi E_0}{k_F a} \sum_{l=1}^{\infty} \sum_{n=1}^{\infty} l \left\{ \frac{l-1/2}{l+1/2} \left( l^2 \left[ 1 - \frac{2E_{n,l}^{(0)}}{mc^2} \frac{l}{l-1/2} \right] f'_{\bar{\mu}}(E_{n,l,(+)}) + (l+1)^2 \left[ 1 - \frac{2E_{n,l}^{(0)}}{mc^2} \frac{l+1}{l+3/2} \right] f'_{\bar{\mu}}(E_{n,l,(-)}) \right) \right. \\ \left. + \left[ \left( \frac{\mathcal{R}_{n,l}}{E_0} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] - \frac{l^2-l+1}{mc^2} \right) f_{\bar{\mu}}(E_{n,l,(+)}) + \left( \frac{l+1}{l+3/2} \frac{\mathcal{R}_{n,l}}{E_0} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] - \frac{l(l+1)}{mc^2} \right) f_{\bar{\mu}}(E_{n,l,(-)}) \right] \right. \\ \left. - \frac{mc^2}{2E_0 E_{n,l}^{(0)}} \left[ 1 - \frac{E_{n,l}^{(0)}}{mc^2} \right] \frac{l+1}{(l+1/2)^2} \left[ f_{\bar{\mu}}(E_{n,l,(+)}) - f_{\bar{\mu}}(E_{n,l,(-)}) \right] \right\}. \quad (85b) $$

In the right-hand side of Eq. (85) we identify the first term of the sums, containing $f'_{\bar{\mu}}(E_{n,l,(\pm)}^{(0)})$, which generalize the nonrelativistic paramagnetic contributions of Eqs. (42) and (43a) including the corresponding weakly-relativistic corrections.

Similarly, the second term of the sums, containing $f'_{\bar{\mu}}(E_{n,l,(\pm)}^{(0)})$, generalizes the nonrelativistic diamagnetic contribution of Eq. (43b). The third term in the sum of Eq. (85b),
containing $f_{\bar{\mu}}(E_{n,l,(+)}^{(0)}) - f_{\bar{\mu}}(E_{n,l,(-)}^{(0)})$, yields a paramagnetic contribution of the van Vleck kind, since it follows from the second-order perturbation correction in $\mathcal{H}^{(\mu)}$ appearing in the coupled basis, together with a weakly-relativistic correction.

Equation (85) can be directly evaluated from the knowledge of the zero-field, weakly-relativistic eigenenergies (65). However, since we are interested in the relativistic corrections to the ZFS, we will consider

$$ \Delta\chi = \chi^{(\text{wr})} - \chi^{(\text{nr})}, \quad (86) $$

where $\chi^{(\text{nr})}$ is defined in Eq. (40), and thus associated with the mean chemical potential $\bar{\mu}_0$ of the nonrelativistic problem, while the weakly relativistic $\chi^{(\text{wr})}$ is associated with the renormalized mean chemical potential

$$ \bar{\mu} = \bar{\mu}_0 + \Delta\mu^{(\text{k})} + \Delta\mu^{(\text{D})}, \quad (87) $$

where $\Delta\mu^{(\text{k})}$ and $\Delta\mu^{(\text{D})}$ are, respectively, given by Eqs. (49) and (C14).

Furthermore, for calculational purposes, it is convenient to treat separately the corrections from the p and e contributions. Thus, we write

$$ \Delta\chi = \Delta\chi_p + \Delta\chi_e, \quad (88) $$

where $\chi_p^{(\text{nr})}$ and $\chi_e^{(\text{nr})}$ have been worked out in Appendix B. In the weakly-relativistic limit, Eqs. (85), (B3), and (B4), together with the form (65) of the zero-field, weakly-relativistic energy correction and that of the renormalized mean chemical potential (87), yield to first-order in $E_0/mc^2$

$$ \frac{\Delta\chi_p}{|\chi_L|} = - \frac{9\pi E_0}{k_F a} \sum_{l=0}^{\infty} (l+1) \sum_{n=1}^{\infty} \left\{ (l+1) f_{\bar{\mu}_0}''(E_{n,l}^{(0)}) + \frac{1}{l+3/2} \frac{\mathcal{R}_{n,l}}{2E_0} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) \right\} \left[ -\frac{(E_{n,l}^{(0)})^2}{2mc^2} + E_{n,l}^{(\text{D})} - \Delta\mu + \frac{E_0 E_{n,l}^{(0)}}{mc^2} l \right] \\ + \frac{9\pi E_0}{k_F a m c^2} \sum_{l=0}^{\infty} (l+1) \sum_{n=1}^{\infty} \left\{ \frac{2(l+1)^2}{l+3/2} E_{n,l}^{(0)} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) + \left( \frac{1}{l+3/2} \frac{\mathcal{R}_{n,l}}{2E_0} E_{n,l}^{(0)} + l+1 \right) f_{\bar{\mu}_0}(E_{n,l}^{(0)}) \right\}, \quad (89\text{a}) $$

$$ \frac{\Delta\chi_e}{|\chi_L|} = - \frac{3\pi E_0}{k_F a} \sum_{l=1}^{\infty} l \sum_{n=1}^{\infty} \left\{ \left[ (2l^2 + 1) f_{\bar{\mu}_0}''(E_{n,l}^{(0)}) + \frac{\mathcal{R}_{n,l}}{E_0} \frac{2l + 5/2}{l + 3/2} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) \right] \left[ -\frac{(E_{n,l}^{(0)})^2}{2mc^2} + E_{n,l}^{(\text{D})} - \Delta\mu \right] \right. \\ \left. + \left[ (l^2 - 1) f_{\bar{\mu}_0}''(E_{n,l}^{(0)}) - \frac{3}{2} \frac{\mathcal{R}_{n,l}}{E_0} \frac{l+1}{l+3/2} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) \right] \frac{E_0 E_{n,l}^{(0)}}{mc^2} \right\} \\ + \frac{3\pi E_0}{k_F a m c^2} \sum_{l=1}^{\infty} l \sum_{n=1}^{\infty} \left\{ E_{n,l}^{(0)} \frac{4l^3 + 6l^2 + l + 1}{l + 3/2} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) + \left[ \frac{2l + 5/2}{l + 3/2} \frac{\mathcal{R}_{n,l}}{E_0} E_{n,l}^{(0)} + 2l^2 + 1 \right] f_{\bar{\mu}_0}(E_{n,l}^{(0)}) \right\}. \quad (89\text{b}) $$

Grouping the two components of Eq. (88), we have

$$ \frac{\Delta\chi}{|\chi_L|} = - \frac{6\pi E_0}{k_F a} \sum_{l=0}^{\infty} \left(l + \frac{1}{2}\right) \sum_{n=1}^{\infty} \left\{ (l^2 + l + 3) f_{\bar{\mu}_0}''(E_{n,l}^{(0)}) + \frac{\mathcal{R}_{n,l}}{E_0} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) \right\} \left[ -\frac{(E_{n,l}^{(0)})^2}{2mc^2} + E_{n,l}^{(\text{D})} - \Delta\mu \right] \\ - \frac{3\pi E_0}{k_F a m c^2} \sum_{l=0}^{\infty} \sum_{n=1}^{\infty} \left\{ 4l \left(l + \frac{1}{2}\right) (l+1) E_0 E_{n,l}^{(0)} f_{\bar{\mu}_0}''(E_{n,l}^{(0)}) - 2(2l^3 + 3l^2 + 5l + 2) E_{n,l}^{(0)} f_{\bar{\mu}_0}'(E_{n,l}^{(0)}) \right. \\ \left. - \left[ 2\left(l + \frac{1}{2}\right) \frac{\mathcal{R}_{n,l}}{E_0} E_{n,l}^{(0)} + 2l^3 + 3l^2 + 7l + 3 \right] f_{\bar{\mu}_0}(E_{n,l}^{(0)}) \right\}. \quad (90) $$

Similarly to the discussion presented after Eq. (53), we remark that the first sum in Eq. (90) corresponds to the zero-field weakly-relativistic correction of the eigenenergies arising from the kinetic and Darwin terms. It could then be directly obtained by implementing, in the nonrelativistic ZFS expressions (42)–(43), the shift from $E_{n,l}^{(0)}$ to $E_{n,l}^{(0)}[1-$

$E_{n,l}^{(0)}/2mc^2] + E_{n,l}^{(\text{D})}$, with the corresponding renormalization $\Delta\mu$ of the chemical potential. The first contribution in the second sum of Eq. (90) corresponds to the effect of $\mathcal{H}^{(\text{so})}$, while the remaining contributions arise from the magnetic-field dependence of $\mathcal{H}^{(\text{k})}$ and $\mathcal{H}^{(\text{ame})}$.
G. Numerical evaluation of the relativistic corrections

When we evaluate numerically the weakly-relativistic correction $\Delta\chi$, according to Eq. (90), as a function of $k_F a$, we find no noticeable difference with the results for the kinetic correction $\Delta\chi^{(k)}$ on the scale of Fig. 2. We therefore conclude that $\Delta\chi$ is dominated by $\Delta\chi^{(k)}$.

As discussed at the end of the last section, it is straightforward to disentangle in $\Delta\chi$ the different contributions arising

from the various Hamiltonian components of Eq. (17), and thus we write

$$ \Delta\chi = \Delta\chi^{(k)} + \Delta\chi^{(D)} + \Delta\chi^{(so)} + \Delta\chi^{(ame)}, \quad (91) $$

where the kinetic correction $\Delta\chi^{(k)}$ is given by Eq. (53). The corrections arising from the energy shift associated with the Darwin term, the spin-orbit coupling, and the angular magneto-electric effect are given, respectively, by

$$ \frac{\Delta\chi^{(D)}}{|\chi_L|} = -\frac{6\pi E_0}{k_F a} \sum_{l=0}^{\infty} \left(l + \frac{1}{2}\right) \sum_{n=1}^{\infty} \left\{ (l^2 + l + 3) f'_{\mu_0}(E_{n,l}^{(0)}) + \frac{R_{n,l}}{E_0} f'_{\mu_0}(E_{n,l}^{(0)}) \right\} [E_{n,l}^{(D)} - \Delta\mu^{(D)}], \quad (92a) $$

$$ \frac{\Delta\chi^{(so)}}{|\chi_L|} = -\frac{12\pi E_0^2}{k_F a mc^2} \sum_{l=1}^{\infty} l \left(l + \frac{1}{2}\right) (l + 1) \sum_{n=1}^{\infty} E_{n,l}^{(0)} f'_{\mu_0}(E_{n,l}^{(0)}), \quad (92b) $$

$$ \frac{\Delta\chi^{(ame)}}{|\chi_L|} = -\frac{12\pi E_0}{k_F a mc^2} \sum_{l=1}^{\infty} \left(l + \frac{1}{2}\right) \sum_{n=1}^{\infty} E_{n,l}^{(0)} f'_{\mu_0}(E_{n,l}^{(0)}). \quad (92c) $$

The behavior of $\Delta\chi^{(k)}$, including its asymptotic dependence in $(k_F a)^0$, has been discussed in Sec. IV C. The numerical evaluation of the other weakly-relativistic corrections to the ZFS, given by Eq. (92), is presented in the three panels of Fig. 3. In all three cases, the typical values are much smaller than those of $\Delta\chi^{(k)}$, and we observe oscillations as a function of $k_F a$ around mean values. Using a one-dimensional semiclassical treatment, analogous to that of the nonrelativistic ZFS (see Appendix A) and the kinetic correction (see Sec. IV C), we can evaluate the leading-order corrections in $k_F a$ in the three cases. For the Darwin, spin-orbit, and angular magneto-electric contributions, in leading-order in $k_F a$, we have, respectively,

$$ \frac{\overline{\Delta\chi}^{(D)}}{|\chi_L|} = -\frac{1}{5} (u^{-1} - 2) \left(\frac{v_F}{c}\right)^2, \quad (93a) $$

$$ \frac{\overline{\Delta\chi}^{(so)}}{|\chi_L|} = -\left(\frac{v_F}{c}\right)^2, \quad (93b) $$

$$ \frac{\overline{\Delta\chi}^{(ame)}}{|\chi_L|} = \left(\frac{v_F}{c}\right)^2, \quad (93c) $$

with $u = \sqrt{E_0/V_0}$.

Using the physical parameters of gold, we have $\overline{\Delta\chi}^{(D)}/|\chi_L| = -0.27 (v_F/c)^2 (k_F a)$, which is in good agreement with the slope associated to the secular behavior of the very small oscillations present in Fig. 3(a). The unbounded behavior obtained for large $a$ is unphysical, since it prevents of achieving the bulk value $\chi^{(wr)}$. As discussed in Appendix C, our perturbative approach is problematic for the Darwin term when treating the discontinuous electric field resulting

from an abrupt electrostatic potential that confines fermions in a reduced (bag) region of space [83].

According to Eqs. (93b) and (93c), we have $\overline{\Delta\chi}^{(so-ame)} = 0$ to leading order in $k_F a$, in agreement with the numerical results [see inset in Fig. 3(c)] and with the expectation that in the infinite-volume limit the spin-orbit and magneto-electric couplings become irrelevant, and Eq. (3) accounts for the weakly-relativistic ZFS. The cancelation between the mean values of $\Delta\chi^{(so)}$ and $\Delta\chi^{(ame)}$ occurs despite the fact that the oscillations of the former are one order of magnitude larger than those of the latter. The comparison between Figs. 2 and 3(b) indicates that $\Delta\chi^{(so)}$ is typically more than an order of magnitude smaller than $\Delta\chi^{(k)}$.

The smallness of the SOC contribution is understandable in view of the high symmetry of the spherical geometry. In Appendix E we present an alternative semiclassical evaluation of this correction, which confirms the results obtained in the previous sections and allows us to highlight the dependence of the SOC effects on the different system parameters. In particular, the spin-orbit correction to the ZFS (E21) exhibits the expected $(v_F/c)^2$ dependence, as well as oscillations between diamagnetic and paramagnetic values with typical amplitudes scaling as $(k_F a)^{1/2}$, which are out-of-phase with respect to the oscillations of the nonrelativistic ZFS. The previous typical amplitudes should be contrasted with the $(k_F a)^{3/2}$ dependence of the nonrelativistic ZFS of an individual nanoparticle and the $k_F a$ dependence resulting from an ensemble average [29]. The suppression of the oscillations of $\Delta\chi^{(so)}$ for large nanoparticle sizes is faster than that of $\Delta\chi^{(nr)}$ since in the first case, in addition to the thermal damping (13) acting on the contribution of each family of classical periodic orbits, we have the fact that the SOC becomes comparatively weaker as $a$ increases. Indeed, from Eqs. (31) and (66) we have $E_{n,l,(\pm)}^{(so)}/E_{n,l}^{(0)} = [\mp(l+1/2)-1/2]E_0/mc^2$, indicating
<figure><img src="image_4.png" /><figcaption>FIG. 3. Weakly-relativistic corrections to the zero-field susceptibility arising from the (a) Darwin, (b) spin-orbit, and (c) angular magneto- electric components, for $k_B T/\bar{\mu}_0 = 5 \times 10^{-3}$, obtained from Eq. (92), as a function of the nanoparticle radius $a$, using the same scaling and physical parameters as in Fig. 2. The dashed lines represent the mean (smooth), leading-order in $k_F a$, values of the different corrections according to Eq. (93). The insets present the corrections $\Delta\chi^{(\text{so})}$ [panel (b)] and $\Delta\chi^{(\text{so-ame})} = \Delta\chi^{(\text{so})} + \Delta\chi^{(\text{ame})}$ [panel (c)] corresponding to large sizes (in the same units as in the main panels), showing the approach to the values predicted by the one-dimensional semiclassical approach given in Eq. (93).</figcaption></figure>

that the relative importance of the SOC decreases with $a$.

## H. Average zero-field susceptibility

In the previous section we established that the SOC con- tribution to the ZFS of an individual nanoparticle is much smaller than its nonrelativistic counterpart, as well as the kinetic-energy correction. Since the measurements are typi- cally not done on individual objects, but on ensembles with

a large number of nanoparticles, we need to consider the av- erage magnetic response, which, as discussed in Sec. II, fol- lows from the finite-$N$ correction to the free energy $\Delta F^{(2)}$ [cf. Eq. (10c)]. The semiclassical description of the spin-orbit coupling, developed in Appendix E as an alternative approach to quantum perturbation theory, becomes extremely useful for calculating $\Delta F^{(2)}$, expressed in terms of the oscillating part of the DOS. Since we are interested in the effect of the SOC, we ignore in the remaining of this section the other weakly- relativistic corrections. Performing the energy integration in Eq. (10c) we write
$$ \Delta F^{(2)} = \frac{1}{2\bar{\rho}(\bar{\mu}, 0)} \left(\frac{\hbar v_F}{E_0}\right)^2 \frac{k_F a}{\pi} \left\{ \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \frac{1}{L_{\nu\eta}} \cos \varphi_{\nu\eta} \sin^{3/2} \varphi_{\nu\eta} R(L_{\nu\eta}/L_T) \right. \\ \left. \times \left[-I_{\nu\eta,c}(\bar{\mu}, B) \sin(\theta_{\nu\eta}(k)) + I_{\nu\eta,s}(\bar{\mu}, B) \cos(\theta_{\nu\eta}(k))\right] \right\}^2, \quad (94) $$

where the symbols $I_{\nu\eta,c}$ and $I_{\nu\eta,s}$ are defined in Eq. (E18). The associated ZFS follows from the magnetic-field dependence of these symbols. To leading order in $v_F/c$ we have $I_{\nu\eta,c}(\bar{\mu}, 0) = 4$, $I_{\nu\eta,s}(\bar{\mu}, 0) = 0$, and $(\partial I_{\nu\eta,c}/\partial B)|_{B=0} = (\partial I_{\nu\eta,s}/\partial B)|_{B=0} = 0$, while the second derivatives are given in Eq. (E20).

The ensemble average ZFS will be dominated by the diagonal terms in the double sum over families of trajectories, and its nonrelativistic component resulting from (94) is given by

$$ \frac{\chi^{(2)-(nr)}}{|\chi_L|} = 36\pi k_F a \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{1}{\eta} \cos^4 \varphi_{\nu\eta} \sin^3 \varphi_{\nu\eta} \\ \times \left[ 1 + \frac{3}{\cos^2 \varphi_{\nu\eta}} \left( \frac{1}{k_F a} \right)^2 \right] \\ \times R^2(L_{\nu\eta}/L_T) \sin^2(\theta_{\nu\eta}(k_F)). \quad (95) $$

$\chi^{(2)-(nr)}$ exhibits a paramagnetic character and the expected period-halving in $k_F a$ with respect to the ZFS of an individual nanoparticle. Its leading-order contribution in $k_F a$ [first term in the square brackets of Eq. (95)] represents the diagonal approximation to the average orbital response calculated in Ref. [29] and reproduced by the red line of Fig. 1 (once the nonrelativistic bulk ZFS $\chi_b^{(nr)}$ is included). The other component [second term in the square brackets of Eq. (95)] stands from the contribution of the Pauli susceptibility. However, it is not significant within our calculation scheme, as it is of an order in $k_F a$ which is not taken into account in the semiclassical approximation of the orbital effects. The remaining contribution to the ZFS resulting from $\Delta F^{(2)}$ in Eq. (94) is the correction $\Delta\chi^{(2)-(so)}$ associated to the SOC, that can be expressed as

$$ \frac{\Delta\chi^{(2)-(so)}}{|\chi_L|} = -\frac{9\pi}{16} \left(\frac{v_F}{c}\right)^2 \sum_{\substack{\nu>0 \\ \eta>2\nu}} \sin^4(2\varphi_{\nu\eta}) \\ \times R^2(L_{\nu\eta}/L_T) \sin(2\theta_{\nu\eta}(k_F)). \quad (96) $$

$\Delta\chi^{(2)-(so)}$ also presents the period-halving in $k_F a$ characteristic of an average response, but it does not have a definite sign. Like in the case of the SOC contribution of an individual nanoparticle discussed in Appendix E, $\Delta\chi^{(2)-(so)}$ is much smaller than its nonrelativistic counterpart (95), due to the prefactor $(v_F/c)^2$ and to one less power of $k_F a$. In a nanoparticle ensemble with an important size dispersion, the oscillating term $\sin^2(\theta_{\nu\eta}(k_F))$ in Eq. (95) averages to $1/2$

and we obtain the result $\chi_{\text{ens}}^d$ plotted as a blue line in Fig. 1 (once the nonrelativistic bulk ZFS $\chi_b^{(nr)}$ is included), while the contribution $\Delta\chi^{(2)-(so)}$ vanishes, and thus does not provide an SOC correction to $\chi_{\text{ens}}^d$.

## V. HALF-SPHERICAL NANOPARTICLES

As discussed at the end of Sec. IV G, the high symmetry of the spherical potential translates into the smallness of the SOC contribution to the ZFS. Thus, a reduction of these symmetries appears as a way to boost the relative importance of the SOC. A first step in the process of progressive destruction of symmetries is to consider a half-spherical confining potential with a magnetic field directed along the axial axis. Such a geometry has the advantage that its nonrelativistic eigenstates at zero magnetic field can be identified as a subset of those of the sphere [87], and many of the analytical developments performed for the case of the sphere can be readily adapted for the case in which the magnetic field is applied along the symmetry axis. Moreover, it is in metallic nanoparticles with the approximate shape of a half-sphere (HS) that the smallest g-factors (as low as 0.3) have been reported [11].

The kinetic and Darwin corrections are expected to change minimally when trading the spherical geometry by the HS, while the angular magneto-electric effect was shown, for the sphere, to be typically one order of magnitude smaller than the one arising from the SOC. Therefore, in this section we concentrate on the effect of the SOC for the ZFS for the reduced symmetry case of half-spherical nanoparticles. Since the total angular momentum is no longer a conserved quantity, the use of the coupled basis is of no aid to treat the SOC, and thus we present our calculations in the product basis, appealing to a numerical diagonalization once the SOC is included in the Hamiltonian.

Similar considerations to those discussed in Sec. IV A, concerning the passage from the ideal model of electrons confined in a sphere to the case of realistic nanoparticles, also apply for the HS geometries studied in this section.

### A. Nonrelativistic susceptibility

We consider noninteracting electrons, described by a Hamiltonian $\mathcal{H}^{(0)}$ of the form (20a), confined by a potential with the shape of a HS, which, with the choice of coordinates
<figure><img src="image_5.png" /><figcaption>FIG. 4. Half-spherical geometry with the coordinate axes used in the text. In red: Classical periodic orbit with the same topological indices $(\nu, \eta) = (1, 5)$ than in the case of Fig. 8 for the sphere, but now contained in two planes (whose intersections with the HS are indicated by the two blue semi-circles).</figcaption></figure>

of Fig. 4, writes

$$V(r, \theta) = V_0 [\Theta(r-a) + \Theta(\theta - \pi/2)] . \quad (97)$$

For finite $V_0$ such a potential does not lead to a separable Schrödinger equation, and it has the unphysical feature of taking two different values ($V_0$ and $2V_0$) in the classically forbidden region. However, these shortcomings are no longer found in the limit $V_0 \to \infty$ of hard walls, where the corresponding eigenstates of $\mathcal{H}^{(0)}$ are characterized by the set of quantum numbers $\{\lambda\} = \{n, l, m_z, m_s\}$ with the condition $l + m_z$ odd, and the associated spinors $\Psi_{n,l,m_z,m_s}^{(\text{HS})}(\mathbf{r})$ have the form (29) with the orbital wave function given by

$$\psi_{n,l,m_z}^{(\text{HS})}(\mathbf{r}) = \sqrt{2} R_{n,l}(r) Y_l^{m_z}(\vartheta) . \quad (98)$$

The radial wave function $R_{n,l}(r)$ is given in Eq. (32). The condition $l + m_z$ being odd arises since only the eigenstates of the sphere which are odd with respect to the reflection off the $z=0$ plane can be eigenstates of the HS, together with the property $Y_l^{m_z}(\pi - \theta, \varphi) = (-1)^{l+m_z} Y_l^{m_z}(\theta, \varphi)$ fulfilled by the spherical harmonics. We remark that the previous restriction excludes the isotropic $l=0$ states. The corresponding eigenenergies $E_{n,l}^{(0)}$ are given by Eq. (31).

In the nonrelativistic case, the application of a magnetic field $\mathbf{B} = B\hat{e}_z$ leads to the same energy corrections of Eqs. (33)–(36). The ZFS of the half-sphere follows from Eq. (6) by summing over $m_z$ and $m_s$ with the restriction of $l+m_z$ being odd. Including the spin-dependent component, we have

$$\frac{\chi^{(\text{nr})} - (\text{HS})}{|\chi_L|} = - \frac{6\pi}{k_F a} \sum_{l=1}^{\infty} \sum_{n=1}^{\infty} l \left\{ \begin{aligned} & (l^2 + 2) E_0 f'_{\mu}(E_{n,l}^{(0)}) \\ & + \frac{(l^2 - 3l/4 - 1)}{(l - 1/2)(l + 3/2)} \mathcal{R}_{n,l} f_{\bar{\mu}}(E_{n,l}^{(0)}) \end{aligned} \right\} . \quad (99)$$

The numerical evaluation of the above expression for the nonrelativistic ZFS of the HS leads to a result (not shown) that is indistinguishable from that of the sphere (see the black solid line in Fig. 1) for the presented $k_F a$-interval, while very small differences appear for smaller values of $k_F a$. For energies corresponding to the first eigenvalues there is a phase shift of half-period in the $k_F a$ oscillations between the ZFS of the sphere and the half-sphere [82], since the lowest-energy state response is always diamagnetic and the HS spectrum starts at the first excited state of the sphere. But this phase shift does not prevail for larger values of $k_F a$.

The concordance of the ZFS of the sphere and the HS for not too small values of $k_F a$ is understandable from a semiclassical viewpoint. Indeed, in Fig. 4 we see that the periodic orbits of the HS can be put in correspondence with those of the sphere by simply symmetrizing one of the two containing planes with respect to the $x-y$ plane. Since the two containing planes are associated with the same angle $\beta$ that defines their orientation with respect to the $z$-axis, Eq. (E14) gives account of the oscillating part of the density of states of the HS (up to a factor of $1/2$ arising from the double counting of the periodic orbits of the sphere). In the absence of SOC the integrations over the angles $\alpha$ and $\gamma$ follow like in the case of the sphere, leading to Eqs. (E17) and (E20a), from where the expression (A19) for the oscillating part of the ZFS $\chi^{(\text{orb})-\text{osc}}$ follows (once a factor of 2 arising from the reduced volume of the HS is included).

## B. Spin-orbit coupling for a HS confining potential

In the hard-wall limit of the potential (97) defining a half-spherical box, the spin-orbit coupling (24a) is given by the Hamiltonian

$$\mathcal{H}^{(\text{so})-(\text{HS})} = \mathcal{H}^{(\text{dome})} + \mathcal{H}^{(\text{floor})}, \quad (100)$$

where

$$\mathcal{H}^{(\text{dome})} = \frac{1}{2m^2c^2} \frac{V_0}{r} \delta(r-a) \mathbf{S} \cdot \mathbf{L}, \quad (101a)$$

$$\mathcal{H}^{(\text{floor})} = \frac{1}{2m^2c^2} \frac{\hbar V_0}{ir} \delta\left(\theta - \frac{\pi}{2}\right) \times \mathbf{S} \cdot \left( -\hat{e}_{\varphi} \partial_r + \hat{e}_r \frac{1}{r \sin\theta} \partial_{\varphi} \right), \quad (101b)$$

represent, respectively, the effect of the dome and the floor of the confining potential.

Writing $\mathbf{S} \cdot \mathbf{L} = (S_+L_- + S_-L_+)/2 + S_zL_z$, we see that the contribution $S_zL_z$ leads to spin-conserving matrix elements of $\mathcal{H}^{(\text{dome})}$ which are diagonal in the indices $m_z$ and $m_s$. That is,

$$\langle \Psi_{n',l',m_z,m_s}^{(\text{HS})} | \mathcal{H}^{(\text{dome})} | \Psi_{n,l,m_z,m_s}^{(\text{HS})} \rangle = (-1)^{(n-n')} m_z m_s \delta_{l',l} \zeta_{n',l} \zeta_{n,l} \frac{2E_0^2}{mc^2}, \quad (102)$$

where we have used that $\operatorname{sgn}[j'_l(\zeta_{n,l})] = (-1)^n$, independently of the value of $l$. The remaining terms of $\mathbf{S} \cdot \mathbf{L}$ result in spin-flip matrix elements of $\mathcal{H}^{(\text{dome})}$ which are nondiagonal in the indices $m_z$ and $m_s$, i.e.,
$$ \langle \Psi_{n', l', m'_z, m'_s}^{(\text{HS})} | \mathcal{H}^{(\text{dome})} | \Psi_{n, l, m_z, m_s = \pm 1/2}^{(\text{HS})} \rangle = (-1)^{(n-n')} \delta_{m'_z, m_z \pm 1} \delta_{m'_s, -m_s} \zeta_{n', l'} \zeta_{n, l} \sqrt{(l \mp m'_z + 1)(l \pm m'_z)} \mathcal{I}_{l', l, m'_z}^{(\text{dome})} \frac{2E_0^2}{mc^2}, \quad (103) $$

where we have defined

$$ \mathcal{I}_{l', l, m'_z}^{(\text{dome})} = \int_{\text{dome}} d\vartheta (Y_{l'}^{m'_z}(\vartheta))^* Y_{l}^{m'_z}(\vartheta). \quad (104) $$

Since $l + m_z$ and $l' + m'_z$ are both odd, the two spherical harmonics in the previous equation have different parities. Therefore, the integral over the dome does not trivially vanish in the case that interests us, where $m'_z = m_z \pm 1$, and then $l$ and $l'$ have a different parity. Equation (F14), in Appendix F, provides the result of the integral (104).

The Hamiltonian component (101b) writes

$$ \mathcal{H}^{(\text{floor})} = \frac{1}{4m^2c^2} \frac{\hbar^2 V_0}{i r} \delta\left(\theta - \frac{\pi}{2}\right) \left( (\sigma_x \sin \varphi - \sigma_y \cos \varphi) \partial_r + \frac{1}{r} (\sigma_x \cos \varphi + \sigma_y \sin \varphi + \sigma_z \cot \theta) \partial_\varphi \right), \quad (105) $$

and thus, the condition $\delta(\theta - \pi/2)$ results in vanishing spin-conserving matrix elements, while the nondiagonal, spin-flip matrix elements are given, in the limit of large $V_0$, by

$$ \langle \Psi_{n', l', m'_z, m'_s}^{(\text{HS})} | \mathcal{H}^{(\text{floor})} | \Psi_{n, l, m_z, m_s = \pm 1/2}^{(\text{HS})} \rangle = \frac{\hbar^2 V_0}{4m^2c^2} \delta_{m'_s, -m_s} \int_0^a dr \, r \int_0^{2\pi} d\varphi \exp(\pm i\varphi) \psi_{m'_z}^{(V_0, n', l')} (r, \theta = \pi/2, \varphi) \\ \times \left\{ \mp \partial_r + \frac{1}{ir} \partial_\varphi \right\} \psi_{m_z}^{(V_0, n, l)} (r, \theta = \pi/2, \varphi). \quad (106) $$

where $\psi_{m_z}^{(V_0, n, l)}(r, \theta, \varphi)$ and $\psi_{m'_z}^{(V_0, n', l')}(r, \theta, \varphi)$ converge, respectively, to the orbital wave functions $\psi_{n, l, m_z}(r, \theta, \varphi)$ and $\psi_{n', l', m'_z}(r, \theta, \varphi)$ when $V_0 \to \infty$. Thus, the limiting condition (F8) and the form (32) of the radial wave function for the hard-wall case allow us to write

$$ \langle \Psi_{n', l', m'_z, m'_s}^{(\text{HS})} | \mathcal{H}^{(\text{floor})} | \Psi_{n, l, m_z, m_s = \pm 1/2}^{(\text{HS})} \rangle = \pm \frac{(-1)^{(l+l'+m_z+m'_z)/2}}{2^{l+l'}} \delta_{m'_z, m_z \pm 1} \delta_{m'_s, -m_s} \frac{4E_0^2}{mc^2} \mathcal{J}_{n', l', n, l, m_z} \\ \times \frac{\sqrt{(2l+1)(2l'+1)}}{|j_{l'+1}(\zeta_{n', l'}) j_{l+1}(\zeta_{n,l})|} \frac{\sqrt{(l+m_z)!(l-m_z)!}}{\left(\frac{l+m_z-1}{2}\right)! \left(\frac{l-m_z-1}{2}\right)!} \frac{\sqrt{(l'+m'_z)!(l'-m'_z)!}}{\left(\frac{l'+m'_z-1}{2}\right)! \left(\frac{l'-m'_z-1}{2}\right)!}, \quad (107) $$

where $l$ and $l'$ have different parities and we have expressed the integral over the radial coordinate through

$$ \mathcal{J}_{n', l', n, l, m_z} = \int_0^1 \frac{d\zeta}{\zeta} j_{l'}(\zeta_{n', l'}\zeta) \left( \frac{l \mp m_z - 1}{\zeta} j_l(\zeta_{n, l}\zeta) - \zeta_{n, l} j_{l+1}(\zeta_{n, l}\zeta) \right). \quad (108) $$

Recalling that $j_l(\zeta) \sim \zeta^l/(2l+1)!!$ for small values of $\zeta$, and that $l=0$ is not allowed for the HS, we verify that the integral in Eq. (108) is divergence-free.

The axial symmetry of the HS translates in the conservation of the z-component of the total angular momentum. The subspaces with different $m_j = m_z + m_s$ are not coupled by $\mathcal{H}^{(\text{so})} - \mathcal{H}^{(\text{HS})}$, and within each of these subspaces the coupling between two eigenstates $\Psi_{n, l, m_z, m_s}^{(\text{HS})}$ and $\Psi_{n', l', m'_z, m'_s}^{(\text{HS})}$ of $\mathcal{H}^{(0)}$ only occurs in the following cases: (i) $\{l', m'_z, m'_s\} = \{l, m_z, m_s\}$; (ii) $m_s = 1/2$ and $\{m'_z, m'_s\} = \{m_z + 1, -1/2\}$; (iii) $m_s = -1/2$ and $\{m'_z, m'_s\} = \{m_z - 1, 1/2\}$ (with $l$ and $l'$ having different parity in the last two cases). These restrictions can be graphically represented in the plane $l-m_z$ (see Fig. 5), where only the subspaces with positive $m_j$

are considered. The diagonalization within each of these subspaces (those with negative $m_j$ follow from Kramer's degeneracy) leads to the eigenstates of $\mathcal{H}^{(0)} + \mathcal{H}^{(\text{so})} - \mathcal{H}^{(\text{HS})}$

$$ \Psi_{m_j, p}(\mathbf{r}) = \sum_{\substack{n, l, m_s = \pm 1/2 \\ l > |m_j - m_s| \\ l + m_j - m_s \neq 2}} C_{n, l, m_s}^{m_j, p} \Psi_{n, l, m_j - m_s, m_s}^{(\text{HS})}(\mathbf{r}), \quad (109) $$

labeled by the half-integer $m_j$ and the positive integer index $p$, with the associated eigenenergies $E_{m_j, p}$.
<figure><img src="image_6.png" /><figcaption>FIG. 5. The $l-m_z$ plane for positive $m_j = m_z + m_s$ for the half-spherical nanoparticle. Each arrow pointing up (down) represents the orientation of the spin $m_s = 1/2$ ($m_s = -1/2$) for the allowed states $l + m_z$ odd and $l > m_z$. At zero-field, the Hamiltonian of the half-spherical nanoparticle only couples states with the same $m_j$.</figcaption></figure>

### C. Perturbative treatment of the magnetic field

Under the application of a weak magnetic field $B$, the energy $E_{m_j,p}$ picks up a correction $\delta E_{m_j,p}$. According to Eq. (6), the latter determines the ZFS of the HS through the parameters

$$ \left. \frac{\partial \delta E_{m_j,p}}{\partial B} \right|_{B=0} = \mu_B \sum_{n,l,m_s} \left| C_{n,l,m_s}^{m_j,p} \right|^2 (m_j + m_s), \qquad (110a) $$

$$ \left. \frac{\partial^2 \delta E_{m_j,p}}{\partial B^2} \right|_{B=0} = \frac{\mu_B^2}{2E_0} \sum_{n,l,m_s} \sum_{n'=1}^{+1} \left( C_{n',l+2i,m_s}^{m_j,p} \right)^* \times C_{n,l,m_s}^{m_j,p} \mathcal{R}_{n',l+2i,n,l} \mathcal{Y}_{l,i}^{m_j-m_s} \times \Theta(l+2i-|m_j-m_s|), \quad (110b) $$

where the sums over $\{n, l, m_s\}$ have the same restrictions as in Eq. (109), and the Heaviside function prevents the consideration of unphysical terms having $l' = l - 2$. The angular matrix elements are calculated over the whole sphere, where the diagonal ones $\mathcal{Y}_{l,0}^{m_z} = \mathcal{Y}_l^{m_z}$ are given by Eq. (38b), while the nonvanishing off-diagonal ($l' = l \pm 2$) ones can be expressed, for $i = \pm 1$, as

$$ \mathcal{Y}_{l,i}^{m_z} = - \frac{1}{4(l + 1/2 + i)} \times \sqrt{\frac{[(l + 1 + i)^2 - m_z^2] [(l + i)^2 - m_z^2]}{(l + 3/2 + i)(l - 1/2 + i)}}. \quad (111) $$

The radial matrix elements are

$$ \mathcal{R}_{n',l',n,l} = \frac{1}{a^2} \int_0^a dr \, r^4 R_{n',l'}(r) R_{n,l}(r). \quad (112) $$

The diagonal ones $\{n', l'\} = \{n, l\}$ are given by Eq. (38a), while the off-diagonal ones can be obtained by numerical integration or by recurrence formulas, as shown in Appendix F3.

### D. Spin-orbit correction for the half-sphere

In Fig. 6 we present the numerical evaluation of the SOC correction $\Delta\chi^{(\text{so})-(\text{HS})}$ to the nonrelativistic ZFS $\chi^{(\text{nr})-(\text{HS})}$, according to Eqs. (6) and (110), as a function of $k_F a$ (red dots). These results are indistinguishable from those corresponding to the SOC correction for the sphere presented in Fig. 3(b) and reproduced by the black line of Fig. 6. The symmetry reduction in going from the sphere to the HS is thus not enough to yield an enhanced SOC correction in the last case. Upon increasing $a$, the SOC mixing of energy levels of the HS with different $(n, l)$ is favored. But such a tendency is countered by the relative weakening of the SOC matrix elements, as discussed at the end of Sec. IV G.

The ZFS correction $\Delta\chi^{(\text{so})-(\text{HS})}$, which is very close to the corresponding correction in the sphere, can be very well approximated by only considering the energy shifts arising from only the diagonal (in $m_z$, $m_s$, and $n$) matrix elements (102), as well as the corresponding diagonal SOC matrix elements for the sphere in the product basis (not shown). These accor-dances indicate that for the sphere, as well as for the HS, for the considered $k_F a$ values, the SOC can be accounted for in first-order perturbation theory.

## VI. CONCLUSION

Motivated by measurements of the magnetization in ensembles of noble-metal nanoparticles [18, 41] and the g-factor of an individual nanoparticle [10, 11], together with the theoretical proposals pointing to the key role played by the spin-orbit coupling in these experimental results [14, 15, 19], we considered in this work the relevance of such a relativistic effect on the magnetic response of confined electrons. In particular, we attempted to quantify the spin-orbit effect by going beyond the previous phenomenological approaches that assigned arbitrary values to its strength. We focused on the extrinsic SOC originating from the confining potential, treating model systems, and then discussing their applicability to the case of metallic nanoparticles of different shapes. The relevance of the SOC in the magnetic response was gauged against other weakly-relativistic corrections and finite-size effects inherent to the mesoscopic regime.

For spherical geometries, the inclusion of the SOC and other relativistic effects could be readily done, at the quantum level, by working in the coupled basis of the total angular momentum, within a perturbative treatment of the magnetic field. An equivalent formulation of the SOC correction was derived from a semiclassical approach including spin and orbital effects.

Analogously to the nonrelativistic zero-field susceptibility characterizing electrons confined in a spherical geometry [29], the weakly-relativistic corrections present oscillations as a function of $k_F a$ between paramagnetic and diamagnetic values, which are typically much larger than the background that sets the bulk response. The oscillations corresponding to the SOC contribution are out-of-phase with respect to those of the nonrelativistic ZFS. The typical values of the SOC contribu-
<figure><img src="image_7.png" /><figcaption>FIG. 6. Red dots: SOC correction $\Delta\chi^{(\text{so})-(\text{HS})}$ to the nonrelativistic ZFS $\chi^{(\text{nr})-(\text{HS})}$ in a half-sphere, from Eqs. (6) and (110), as a function of $k_F a$. Black line: reproduction of the result of Fig. 3(b) obtained for the sphere using Eq. (92b).</figcaption></figure>

tion are more than an order of magnitude larger than the ones due to the magneto-electric coupling, and more than an order of magnitude smaller than those arising from the kinetic energy correction. The latter, being the dominant weakly-relativistic correction to the ZFS, remains much smaller than the nonrelativistic ZFS. These small values of the weakly-relativistic corrections stem from their $(v_F/c)^2$ dependence, while the reduced effect of the SOC correction is associated with the high spatial symmetry of the spherical geometry.

In order to study the impact of a symmetry reduction on the magnetic response, the formalism developed to treat the SOC in the sphere was adapted to the case of a half-sphere. The nonrelativistic ZFS was shown to be the same for the sphere and the half-sphere in the semiclassical limit of large $k_F a$, as can be readily understood from semiclassical arguments. The inclusion of the SOC in the HS leads to corrections of the ZFS which are very close to those of the sphere, and therefore much smaller than the correction arising from the kinetic energy shift and the typical values of the nonrelativistic ZFS. The symmetry reduction when going from the spherical geometry to that of the HS is not enough to render relevant the SOC effects due to the electron confinement. Even if the HS confinement favors the SOC mixing of the unperturbed states in comparison with the case of the sphere, such an effect is offset by the generic suppression of the SOC matrix elements with the parameter $(v_F/c)^2$ and the size $a$. The SOC is thus weak enough to remain perturbative, and therefore it does not induce the transition in the statistical properties of the spec-

tra necessary to reverse the sign of the magnetic response of a nanoparticle ensemble.

Further symmetry reduction would lead to more quasi-degeneracies in the $B=0$ spectrum, enhancing the effect of the SOC. However, the cases of the quarter- and eighth-sphere, which can be worked out similarly to that of the HS, do not show a significant increase of the SOC contribution to the ZFS [88]. The case of cylindrical nanoparticles is important, as it is the geometry considered in Ref. [18]. At least in the most symmetric case where the magnetic field is oriented along the cylinder axis, the SOC due to the confinement seems to have minor importance [89]. Moving forward from integrable to chaotic geometries will increase the importance of the SOC mixing, but at the same time will be accompanied by a reduction of the typical values attained by the nonrelativistic ZFS [53, 90].

Our theoretical results show that the spin-orbit coupling generated by the potential discontinuity at the border of a metallic nanoparticle remains far away from the strong-coupling regime [19] where the change of the statistical properties of the spectrum may lead to a diamagnetic response for the case of disordered or classically chaotic systems. The nanorods of Ref. [18] are considerably larger than the nanoparticles of Refs. [22, 26, 30, 31] which might result in a stronger extrinsic SOC effect arising from the impurity potentials. However, the high-quality single-crystalline gold employed is not expected to result in important extrinsic effects, nor in an underlying chaotic classical dynamics of the con-
duction electrons. While the SOC contribution to the zero-field susceptibility of an individual nanoparticle is suppressed by a factor $(v_F/c)^2$ with respect to its nonrelativistic counterpart, the SOC contribution to the average response for ensembles of nanoparticles with an important size dispersion vanishes. As discussed in the introduction, the other mechanisms of SOC (host ionic lattice and impurities) are likewise expected to be ineffective towards driving the transition between statistical ensembles. We thus conclude that spin-orbit effects are by themselves not responsible for the large (with respect to the bulk) diamagnetic susceptibility observed in Refs. [18, 22, 26, 30, 31].

The magnetic response of micro- and nanostructures is characterized by the magnetic susceptibility in the singly-connected case and by the persistent current in multiply-connected geometries. Our work in the first problem contributes to the open discussion about the possible origin of the diamagnetic response of a large ensemble of metallic rings [62]. While the typical values of the persistent current measured in the second-generation experiments [64, 91, 92] agree with the theoretical predictions [56–59], the observed average diamagnetic persistent current forces us to quantify the role of potentially relevant physical effects, such as the spin-orbit coupling in constrained geometries or the attractive superconducting-like electron-electron interaction, the first of which constituted the object of the present study.

## ACKNOWLEDGMENTS

We thank Hélène Bouchiat, Gaëtan Percebois, and Pablo Tamborenea for stimulating discussions. We acknowledge financial support from the French National Research Agency (ANR) through Grant No. ANR-14-CE26-0005 Q-MetaMat. This work of the Interdisciplinary Thematic Institute QMat, as part of the ITI 2021–2028 program of the University of Strasbourg, CNRS, and Inserm, was supported by IdEx Unistra (ANR 10 IDEX 0002), and by SFRI STRAT'US project (ANR 20 SFRI 0012) and ANR-17-EURE-0024 under the framework of the French Investments for the Future Program.

## Appendix A: Semiclassical evaluation of the zero-field orbital susceptibility

In this appendix, we use the semiclassical approximation for the $l$-fixed DOS $\varrho_l(E)$ in order to calculate the contributions of Eq. (43) to the nonrelativistic ZFS. We also show the equivalence between this procedure with the use of the 3D semiclassical approximation employed in Ref. [29].

### 1. One-dimensional semiclassical approach

Similarly to the decomposition of the DOS into a smooth (Weyl) and an oscillatory (in energy) part, presented in Eq. (7) for the field-dependent case, the semiclassical approximation

to the radial problem allows to write the field-independent $l$-fixed spinless DOS as [70]

$$ \varrho_l(E) = \bar{\varrho}_l(E) + \varrho_l^{\text{osc}}(E), \quad (A1) $$

with

$$ \bar{\varrho}_l(E) = \frac{\tau_l(E)}{2\pi\hbar}, \quad (A2a) $$

$$ \varrho_l^{\text{osc}}(E) = \frac{\tau_l(E)}{\pi\hbar} \sum_{\eta=1}^{\infty} \cos\left(\eta \left[ \frac{S_l(E)}{\hbar} - \frac{3\pi}{2} \right]\right). \quad (A2b) $$

The sum is over the repetitions $\eta$ of the periodic orbit in the effective ($l$-dependent) 1D radial potential

$$ V_l^{\text{eff}}(r) = V(r) + \frac{\hbar^2}{2m} \frac{(l+1/2)^2}{r^2}, \quad (A3) $$

having an energy $E$, a period

$$ \tau_l(E) = \frac{\hbar}{E} \sqrt{\frac{E}{E_0} - \left(l + \frac{1}{2}\right)^2}, \quad (A4) $$

and a classical action

$$ S_l(E) = 2\hbar \left[ \sqrt{\frac{E}{E_0} - \left(l + \frac{1}{2}\right)^2} - \left(l + \frac{1}{2}\right) \arccos \left( \frac{l + 1/2}{\sqrt{E/E_0}} \right) \right]. \quad (A5) $$

For a given $E$, the maximum allowed $l$ is given by

$$ l_{\max} = \left\lfloor \sqrt{\frac{E}{E_0}} - \frac{1}{2} \right\rfloor, \quad (A6) $$

where the floor function $\lfloor \zeta \rfloor$ yields the integer part of $\zeta$.

The radial problem, being a one-dimensional configuration, results in a $\bar{\varrho}_l(E)$ with a $1/\sqrt{E}$ behavior for $E/E_0 \gg (l+1/2)^2$ and the same $\hbar$ dependence than $\varrho_l^{\text{osc}}(E)$. The decomposition (A1), when used in Eq. (45), reproduces [70, 93] the well-known zero-field semiclassical DOS of the sphere [94, 95], with a smooth part $g^{(3D)}(E)\mathcal{V}$ and an oscillating one, that like Eq. (8), has the form

$$ \varrho^{\text{osc}}(E, 0) = \sum_{\nu=1}^{\infty} \sum_{\eta=2\nu+1}^{\infty} \varrho_{\nu\eta}(E), \quad (A7) $$

where the spin degeneracy factor is included, and

$$ \varrho_{\nu\eta}(E) = \frac{4}{E_0} \sqrt{\frac{ka}{\pi}} \frac{(-1)^\nu}{\sqrt{\eta}} \cos \varphi_{\nu\eta} \sin^{3/2} \varphi_{\nu\eta} \cos \left( \theta_{\nu\eta}(k) \right). \quad (A8) $$

The sum is over families of classical periodic orbits lying on the equatorial plane of the sphere, labeled by the topological indexes $(\nu, \eta)$, with $\nu$ the number of turns around the center
(i.e., the winding number) and $\eta$ the number of specular reflections at the boundary (i.e., the number of bounces).<fnref n="7" /> We note $k = \sqrt{2mE}/\hbar$ and define

$$ \varphi_{\nu\eta} = \frac{\pi\nu}{\eta}, \qquad (A9) $$

corresponding to half the angle spanned between two consecutive bounces (see Fig. 7). Furthermore, we have defined the $k$-dependent phase

$$ \theta_{\nu\eta}(k) = kL_{\nu\eta} + \frac{\pi}{4} - \frac{3\pi}{2}\eta, \qquad (A10) $$

with $L_{\nu\eta} = 2\eta a \sin \varphi_{\nu\eta}$ the trajectory length.

In Ref. [95] it is shown that the effect of a magnetic field weak enough to produce a cyclotron radius much larger than $a$ is to affect each term of the sum (A7) representing the oscillating part of the density of states by a field-dependent orbital modulation factor for the $(\nu, \eta)$ family, given by

$$ M_{\nu\eta}^{(\text{orb})}(B) = j_0(2\pi\phi_{\nu\eta}/\phi_0). \qquad (A11) $$

We note $j_0(\zeta) = \sin \zeta / \zeta$ the zeroth-order spherical Bessel function of the first kind, while

$$ \phi_{\nu\eta} = B A_{\nu\eta} \qquad (A12) $$

is the magnetic flux enclosed by the orbits $(\nu, \eta)$ covering the area $A_{\nu\eta} = \frac{1}{2}\eta a^2 \sin(2\varphi_{\nu\eta})$, and $\phi_0 = hc/e$ is the flux quantum.

## 2. Equivalence between 1D and 3D semiclassical

The semiclassical evaluation of $\chi^{(\text{orb})}$ can be addressed by trading in Eqs. (43) the sums over the principal quantum number $n$ by energy-integrals and the use of the Poisson summation formula for the sum over $l$, resulting in

$$ \frac{\chi^{(\text{para})}}{|\chi_L|} = - \frac{6\pi E_0}{k_F a} \int_0^\infty dE f'_{\mu_0}(E) \sum_{\nu=-\infty}^{+\infty} \int_{-1/2}^{\sqrt{E/E_0}-1/2} dl \times \exp(2\pi i \nu l) l \left(l+\frac{1}{2}\right) (l+1) \varrho_l(E), \qquad (A13a) $$

$$ \frac{\chi^{(\text{dia})}}{|\chi_L|} = - \frac{6\pi}{k_F a} \int_0^\infty dE f_{\mu_0}(E) \sum_{\nu=-\infty}^{+\infty} \int_{-1/2}^{\sqrt{E/E_0}-1/2} dl \times \exp(2\pi i \nu l) \left(l + \frac{1}{2}\right) \mathcal{R}_l(E) \varrho_l(E), \qquad (A13b) $$

where $l$ is now understood as a continuous variable. Following Eq. (38a), we have defined $\mathcal{R}_l(E) = (1/3)[1 + (2E_0/E)(l + 3/2)(l - 1/2)]$.

In the leading order in $k_F a \gg 1$, the smooth part of $\chi^{(\text{orb})}$ is obtained by using $\bar{\varphi}_l(E)$ in Eqs. (A13) and only keeping the $\nu = 0$ term of the sum,

$$ \frac{\bar{\chi}^{(\text{para})}}{|\chi_L|} = - \frac{3}{k_F a} \int_0^\infty dE f'_{\mu_0}(E) \left(\frac{E}{E_0}\right)^{1/2} \times \int_0^1 d\zeta \zeta \left(\frac{E}{E_0}\zeta^2 - \frac{1}{4}\right) \sqrt{1 - \zeta^2}, \qquad (A14a) $$

$$ \frac{\bar{\chi}^{(\text{dia})}}{|\chi_L|} = - \frac{1}{k_F a E_0} \int_0^\infty dE f_{\mu_0}(E) \left(\frac{E}{E_0}\right)^{1/2} \times \int_0^1 d\zeta \zeta \left(1 + 2\zeta^2 - \frac{2}{E/E_0}\right) \sqrt{1 - \zeta^2}. \qquad (A14b) $$

Performing the integration over the variable $\zeta = (l + 1/2)\sqrt{E_0/E}$ we have

$$ \frac{\bar{\chi}^{(\text{para})}}{|\chi_L|} = - \frac{1}{k_F a} \int_0^\infty dE f'_{\mu_0}(E) \left(\frac{E}{E_0}\right)^{1/2} \left(\frac{2E}{5E_0} - \frac{1}{4}\right), \qquad (A15a) $$

$$ \frac{\bar{\chi}^{(\text{dia})}}{|\chi_L|} = - \frac{1}{k_F a E_0} \int_0^\infty dE f_{\mu_0}(E) \left(\frac{E}{E_0}\right)^{1/2} \left(\frac{3}{5} - \frac{2E_0}{3E}\right). \qquad (A15b) $$

An integration by parts in Eq. (A15b) leads to

$$ \frac{\bar{\chi}^{(\text{dia})}}{|\chi_L|} = \frac{1}{k_F a} \int_0^\infty dE f'_{\mu_0}(E) \left(\frac{E}{E_0}\right)^{1/2} \left(\frac{2E}{5E_0} - \frac{4}{3}\right). \qquad (A16) $$

Thus, the leading-order term in $k_F a$ of $\bar{\chi}^{(\text{para})}$ and $\bar{\chi}^{(\text{dia})}$ cancel each other, and $\bar{\chi}^{(\text{orb})}$ is then given by next-order contributions.<fnref n="8" /> However, such terms are not captured by the expressions of Eqs. (A15), which are only valid in the leading order in $k_F a$, since they result from $E$ and $l$ integrations in which the form (A4) was used beyond its regime of validity of $E/E_0 \gg (l + 1/2)^2$. In Ref. [53] it is shown that the magnetic field dependence of $\bar{\varphi}(E, B)$, given by the so-called zero-length trajectories, results in $\bar{\chi}^{(\text{orb})} = \chi_L$, while the numerical evaluation of Eqs. (42)–(43) presented in Fig. 1 approaches, in the limit of large radius $a$ where the role of the confinement potential becomes irrelevant, the bulk result given by Eqs. (1)–(2) (see the inset in Fig. 1).

The fact that $\bar{\chi}^{(\text{orb})}/|\chi_L|$ is of order $(k_F a)^0$, points to the importance of $\chi^{\text{osc}}$, which is obtained by using $\varrho_l^{\text{osc}}(E)$ in Eq. (A13). The rapidly oscillating (in $E$) phases $2\pi\nu \pm$

<footnote n="7">⁷ We have not included the contribution of the diametral (self-retracing) orbits with $\eta = 2\nu$ having only a twofold degeneracy, and thus associated with a contribution to the DOS which is of higher order in $\hbar$ than the ones of Eq. (A7). Moreover, diametral trajectories not enclosing a magnetic flux do not contribute to the ZFS, and as we will see in the sequel, the SOC is ineffective in their cases.</footnote>

<footnote n="8">⁸ For the degenerate case, in the leading order in $k_F a$, we have that $\bar{\chi}^{(\text{para})} = |\bar{\chi}^{(\text{dia})}| = (2/5) (k_F a)^2$. Notwithstanding, we stress that the separation between $\chi^{(\text{para})}$ and $\chi^{(\text{dia})}$ is only for computational purposes, and lacks physical reality.</footnote>
$\eta S_l(E)/\hbar$ allow us to perform a stationary-phase (sp) integration over $l$ with the condition $l_{sp} = \sqrt{E/E_0} \cos(\pi\nu/\eta) - 1/2$, yielding

$$ \frac{\chi^{(\text{para})-\text{osc}}}{|\chi_L|} = - \frac{6\sqrt{\pi}}{k_F a} \int_0^\infty dE f'_{\bar{\mu}_0}(E) \left(\frac{E}{E_0}\right)^{1/4} \\ \times \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \sin^{3/2} \varphi_{\nu\eta} \cos \varphi_{\nu\eta} \\ \times \left[ \frac{E}{E_0} \cos^2 \varphi_{\nu\eta} - \frac{1}{4} \right] \cos \left(\theta_{\nu\eta}(k)\right), \qquad \text{(A17a)} $$

$$ \frac{\chi^{(\text{dia})-\text{osc}}}{|\chi_L|} = - \frac{2\sqrt{\pi}}{k_F a E_0} \int_0^\infty dE f_{\bar{\mu}_0}(E) \left(\frac{E}{E_0}\right)^{1/4} \\ \times \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \sin^{3/2} \varphi_{\nu\eta} \cos \varphi_{\nu\eta} \\ \times \left[ 1 + 2 \cos^2 \varphi_{\nu\eta} - \frac{2E_0}{E} \right] \cos \left(\theta_{\nu\eta}(k)\right). \qquad \text{(A17b)} $$

The couple of integers $(\nu, \eta)$ of the one-dimensional problem, introduced respectively in Eqs. (A13) and (A2b), take now the role of the topological indexes labeling the families of classical periodic orbits in the sphere. The restriction of $\nu \neq 0$ appears since these contributions, considered in Eqs. (A14), lead to $\bar{\chi}^{(\text{orb})}$. Only positive values of $\nu$ are kept, as the sign of $\nu$ is associated with the orientation in which the periodic orbit is traveled. The condition $\eta > 2\nu$ appears as a restriction for the stationary-phase value $l_{sp}$ to be within the integration interval. The special case of $\eta = 2\nu$, where $l_{sp}$ is at the lower limit of the integration interval, corresponds to the diametral orbits (see footnote 7). The angle $\varphi_{\nu\eta}$ and the $k$-dependent phase $\theta_{\nu\eta}$ have been defined in Eqs. (A9) and (A10), respectively. The stationary-phase procedure yielding (A17) is analogous to that allowing to link $\varrho_l^{\text{osc}}(E)$ and $\varrho^{\text{osc}}(E, 0)$ [93]. Since we work in leading-order in $k_F a$, we can neglect the last terms in the square brackets of Eqs. (A17).

In the low-temperature limit, we use $f_{\bar{\mu}_0}(E) = \Theta(\bar{\mu}_0 - E)$ and, with the help of Fresnel integrals, to leading order in $k_F a$, we find

$$ \frac{\chi^{(\text{dia})-\text{osc}}}{|\chi_L|} = -2\sqrt{\pi}(k_F a)^{1/2} \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\eta^{3/2}} \sin^{1/2} \varphi_{\nu\eta} \cos \varphi_{\nu\eta} \\ \times (1 + \cos^2 \varphi_{\nu\eta}) \sin \left(\theta_{\nu\eta}(k_F)\right). \qquad \text{(A18)} $$

Thus, $\chi^{(\text{dia})-\text{osc}} \sim (k_F a)^{1/2}$ is of lower order in $k_F a$ than $\chi^{(\text{para})-\text{osc}} \sim (k_F a)^{3/2}$, and in the semiclassical limit we have $\chi^{(\text{orb})-\text{osc}} = \chi^{(\text{para})-\text{osc}}$.

The integration of $f'_{\bar{\mu}_0}(E)$ multiplied by a rapidly oscillating function of $E$ results in a general expression [53], that

when applied to Eq. (A17a) yields

$$ \frac{\chi^{(\text{orb})-\text{osc}}}{|\chi_L|} = 6\sqrt{\pi}(k_F a)^{3/2} \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \sin^{3/2} \varphi_{\nu\eta} \cos^3 \varphi_{\nu\eta} \\ \times R(L_{\nu\eta}/L_T) \cos \left(\theta_{\nu\eta}(k_F)\right), \qquad \text{(A19)} $$

which is the result obtained in Ref. [29] starting from the three-dimensional DOS (A7) with the modulation factor (A11). The thermal factor $R$, given by Eq. (13), depends on the ratio $L/L_T$, where $L$ is the trajectory length and $L_T = \hbar v_F / \pi k_B T$ the thermal length.

When $\chi_b^{(\text{nr})}$ is added to the semiclassical evaluation of Eq. (A19), the result is almost indistinguishable from the non-relativistic ZFS shown in Fig. 1 (in black), obtained from the quantum perturbative Eqs. (42) and (43).

## Appendix B: Nonrelativistic susceptibility using the coupled basis

As discussed in Secs. IVD and IV E, the treatment of the SOC is greatly simplified by working in the coupled basis of total angular momentum within the subspace decomposition (61). It is therefore useful to recast the nonrelativistic ZFS as the sum of two components arising from the p (product) and e (entangled) states associated, respectively, with the subspaces $S_{n,l,(d/u)}^P$ and $S_{n,l,m_j}^e$, writing

$$ \chi^{(\text{nr})} = \chi_{\text{p}}^{(\text{nr})} + \chi_{\text{e}}^{(\text{nr})}. \qquad \text{(B1)} $$

The p-states (d/u) are characterized, respectively, by the quantum numbers $\{n, l, m_z = -l, m_s = -1/2\}$ and $\{n, l, m_z = l, m_s = +1/2\}$ of the product basis, that can be expressed in a compact way as $\{n, l, m_z = \mp l, m_s = (1/2)(m_z/l)\}$. According to Eq. (33), at finite magnetic field, these states are associated with the energies

$$ E_{n,l}^{(\text{nr})} (d/u) = E_{n,l}^{(0)} \mp \mu_B B(l+1) + E_{n,l,\mp l}^{(\text{dia})}. \qquad \text{(B2)} $$

Applying Eq. (6) and summing over $m_z = \mp l$ we have

$$ \frac{\chi_{\text{p}}^{(\text{nr})}}{|\chi_{\text{L}}|} = - \frac{9\pi E_0}{k_F a} \sum_{l=0}^{\infty} \sum_{n=1}^{\infty} \left\{ (l+1)^2 f'_{\bar{\mu}_0} (E_{n,l}^{(0)}) + \frac{\mathcal{R}_{n,l}}{2E_0} \frac{l+1}{l+3/2} f_{\bar{\mu}_0} (E_{n,l}^{(0)}) \right\}. \qquad \text{(B3)} $$

In the nonrelativistic case, the subspace $S_{n,l,m_j}^e$ can be characterized by the product-basis states with quantum numbers $\{n, l, m_z = \mp l, m_s = -(1/2)(m_z/l)\}$ and $\{n, l, m_z \in [-l+1, l-1], m_s = \pm 1/2\}$. For each couple $(n, l)$, there are two states of the first kind and $4l-2$ of the second one, whose energies at finite magnetic field follow from Eq. (33). Once the sum over $m_z$ is performed in Eq. (6), they yield

$$ \frac{\chi_e^{(\text{nr})}}{|\chi_L|} = - \frac{3\pi E_0}{k_F a} \sum_{l=1}^{\infty} \sum_{n=1}^{\infty} l \left\{ (2l^2 + 1) f'_{\bar{\mu}_0} (E_{n,l}^{(0)}) + \frac{\mathcal{R}_{n,l}}{E_0} \frac{2l + 5/2}{l + 3/2} f_{\bar{\mu}_0} (E_{n,l}^{(0)}) \right\}. \qquad \text{(B4)} $$
Obviously, the addition of Eqs. (B3) and (B4) results in the expression of $\chi^{(nr)}$ presented in Sec. IV B [cf. Eq. (40)].

## Appendix C: Darwin correction

In this appendix, we consider the Darwin correction for a spherically symmetric confinement, given by the Hamiltonian contribution $\mathcal{H}^{(D)}$ of Eq. (28). In order to address the corresponding matrix elements, we first tackle a technical issue concerning the eigenvalues and eigenvectors of a spherical box confined by a finite-height potential.

### 1. Finite-box spherical potential

The spherical symmetry of the potential (26) allows to write the eigenstates as in Eqs. (29)–(30), where the radial wave function of the bound states is given by

$$R_{n,l}(r) = \sqrt{\frac{2}{a^3}} \frac{1}{|c_{n,l}|^{1/2}} \begin{cases} j_l(k_{n,l} r), & r \le a, \\ \frac{j_l(k_{n,l} a)}{k_l(\kappa_{n,l} a)} k_l(\kappa_{n,l} r), & r > a. \end{cases} \quad (C1)$$

We follow the standard convention of using $j_l(\zeta)$ for the spherical Bessel function of the first kind and order $l$. Similarly, we note $k_l(\zeta)$ the modified spherical Bessel function of the second kind and order $l$, obtained from the imaginary-argument spherical Hankel function $h_l^{(1)}(\zeta)$, i.e., $k_l(\zeta) = -i^l h_l^{(1)}(i\zeta)$. The condition of a bound state implies that its energy $E_{n,l}^{(0)} = E_0(k_{n,l} a)^2$ is smaller than $V_0$. We have defined $E_0 = \hbar^2/2ma^2$ and $\kappa_{n,l} = \sqrt{2m(V_0 - E_{n,l}^{(0)})/\hbar}$. The normalization is settled by the constant

$$c_{n,l} = \left( \frac{j_l(k_{n,l} a)}{k_l(\kappa_{n,l} a)} \right)^2 k_{l-1}(\kappa_{n,l} a) k_{l+1}(\kappa_{n,l} a) - j_{l-1}(k_{n,l} a) j_{l+1}(k_{n,l} a). \quad (C2)$$

The allowed $E_{n,l}^{(0)}$ and $k_{n,l}$ result from the solutions of the quantization condition

$$\kappa_{n,l} a \frac{\mathbf{k}'_l(\kappa_{n,l} a)}{\mathbf{k}_l(\kappa_{n,l} a)} = k_{n,l} a \frac{j'_l(k_{n,l} a)}{j_l(k_{n,l} a)}, \quad (C3)$$

which can be recast as

$$\kappa_{n,l} a \frac{\mathbf{k}_{l+1}(\kappa_{n,l} a)}{\mathbf{k}_l(\kappa_{n,l} a)} = k_{n,l} a \frac{j_{l+1}(k_{n,l} a)}{j_l(k_{n,l} a)}, \quad (C4)$$

once we employ the useful recurrence relationship

$$q'_l(\zeta) = \frac{l}{\zeta} q_l(\zeta) - q_{l+1}(\zeta), \quad (C5)$$

valid for $q_l(\zeta) = j_l(\zeta)$, as well as for $q_l(\zeta) = k_l(\zeta)$.

In the limiting case of a hard-wall potential ($V_0 \to \infty$), the support of the wave function is $r < a$, and the previous expressions result in

$$k_{n,l} a \to \zeta_{n,l}, \quad (C6a)$$

$$c_{n,l} \to [j_{l+1}(\zeta_{n,l})]^2, \quad (C6b)$$

where $\zeta_{n,l}$ stands for the nth root of $j_l(\zeta)$, and thus the radial wave function (C1) takes the simpler form of Eq. (32).

In the case where $V_0$ is large but remains finite, the second-order expansion of (C4) around $\zeta_{n,l}$ in the small parameter $u = \sqrt{E_0/V_0}$ allows to write the corrections to Eq. (C6) and to other important parameters as

$$k_{n,l} a \approx \zeta_{n,l} (1 - u + u^2), \quad (C7a)$$

$$c_{n,l} \approx [j_{l+1}(\zeta_{n,l})]^2 (1 + 3u + 3u^2), \quad (C7b)$$

$$j_l(k_{n,l} a) \approx j_{l+1}(\zeta_{n,l}) \zeta_{n,l} u, \quad (C7c)$$

$$j_l(k_{n,l} a) j'_l(k_{n,l} a) \approx -[j_{l+1}(\zeta_{n,l})]^2 \zeta_{n,l} (u + 2u^2), \quad (C7d)$$

$$\kappa_{n,l} a \approx u^{-1} - \frac{1}{2} \zeta_{n,l}^2 u + \zeta_{n,l}^2 u^2. \quad (C7e)$$

It is important to remark that Eq. (C7c) is valid up to quadratic order in $u$, and this limiting condition leads to the expression (59) of the SOC radial matrix element, which does not depend on $V_0$.

### 2. Matrix elements for the Darwin correction

According to Eq. (28), the diagonal Darwin energy correction is

$$E_{n,l}^{(D)} = \frac{E_0}{4mc^2} I_{n,l}^{(D)}, \quad (C8)$$

with the radial matrix element

$$I_{n,l}^{(D)} = a^2 \int_0^\infty dr [R_{n,l}(r)]^2 \frac{d}{dr} \left( r^2 \frac{dV}{dr} \right). \quad (C9)$$

For the potential (26) we have $V'(r) = V_0 \delta(r-a)$ and therefore the integral (C9) can be trivially performed, yielding

$$I_{n,l}^{(D)} = - V_0 a^4 \left. \frac{d}{dr} [R_{n,l}(r)] \right|_{r=a}, \quad (C10)$$

that, using the form (C1) of the radial wave function, results in

$$I_{n,l}^{(D)} = -4V_0 \frac{k_{n,l} a}{|c_{n,l}|} j_l(k_{n,l} a) j'_l(k_{n,l} a). \quad (C11)$$

According to Eq. (C7c), the radial matrix element $I_{n,l}^{(D)}$ diverges as $\sqrt{V_0}$ in the large $V_0$-limit. Such a divergence is unphysical since an infinite $V_0$ would imply an infinite electric field, for which the weakly-relativistic approach leading to Eq. (14) would not be valid. The subtleties related with a relativistic particle hitting a steep wall are extensively discussed
in the literature [83, 84]. The appearance of a divergence in our weakly-relativistic approach calls for a systematic expansion in the small parameter $u = \sqrt{E_0/V_0}$, as presented in the first part of this appendix. In addition, we need to verify that the physical constants of the problem are such that we work in the regime of validity of the perturbative approach. For the case of gold nanoparticles, the values of $E_F = 5.5$ eV and $W = 4.3$ eV [80] result in $u \simeq 0.75 (k_F a)^{-1}$. Therefore, in the semiclassical limit of $k_F a \gg 1$, to which our study is restricted, $u$ is indeed a small parameter.

Using the second-order expressions of Eq. (C7), the radial matrix element (C11) can be approximated by

$$I_{n,l}^{(D)} = 4V_0 \zeta_{n,l}^2 u (1 - 2u), \quad (C12)$$

and thus, the Darwin energy correction is

$$E_{n,l}^{(D)} = \frac{E_0^2 \zeta_{n,l}^2}{mc^2} (u^{-1} - 2). \quad (C13)$$

Notice that the forms (C13) and (60), respectively, of the Darwin and spin-orbit energy corrections, are both valid up to terms of order $(V_0)^0$. We remark that, while $\mathcal{H}^{(D)}$ couples states with different $n$, the resulting second-order corrections in $v_F/c$ can be neglected.

Similarly to the case of the kinetic correction treated in Sec. IV C, the Darwin energy (C13) induces, at $B = 0$, a

renormalization of the chemical potential

$$\Delta\mu^{(D)} = (u^{-1} - 2) \frac{\bar{\mu}_0 E_0}{mc^2}. \quad (C14)$$

## Appendix D: Matrix elements for the angular magneto-electric coupling

In this appendix we calculate the matrix elements of the angular magneto-electric Hamiltonian $\mathcal{H}^{(ame)}$ restricted to the subspace $S_{n,l,m_j}^e$, as discussed in Sec. IV E. According to Eq. (27b) the diagonal matrix element is

$$\begin{aligned} \mathcal{E}_{n,j,m_j,(\pm)}^{(ame)} &= \langle \Phi_{n,j,m_j}^{(\pm)} | \mathcal{H}^{(ame)} | \Phi_{n,j,m_j}^{(\pm)} \rangle \\ &= \frac{\mu_B B}{4mc^2} I_{n,j \pm 1/2}^{(\text{so})} \mathcal{I}_{j,m_j}^{(\text{d-ame})(\pm)}. \end{aligned} \quad (D1)$$

The radial matrix element coincides with that arising from the spin-orbit coupling (59) since, in the hard-wall limit,

$$\int_0^a dr \, r^3 [R_{n,l}(r)]^2 V'(r) = I_{n,l}^{(\text{so})}. \quad (D2)$$

Using the standard notation $\sigma_{\pm} = \sigma_x \pm i\sigma_y$, the angular matrix element is

$$\mathcal{I}_{j,m_j}^{(\text{d-ame})(\pm)} = \langle \Upsilon_{j,m_j}^{(\pm)} | \sin^2 \theta \sigma_z - \frac{\sin \theta \cos \theta}{2} (e^{-i\varphi}\sigma_+ + e^{i\varphi}\sigma_-) | \Upsilon_{j,m_j}^{(\pm)} \rangle. \quad (D3)$$

The spin-conserving component of $\mathcal{I}_{j,m_j}^{(\text{d-ame})(+)}$ is given by

$$\langle \Upsilon_{j,m_j}^{(+)} | \sin^2 \theta \sigma_z | \Upsilon_{j,m_j}^{(+)} \rangle = -\frac{m_j}{2} \left( \frac{j(j+2)+j+1+m_j^2}{j(j+1)(j+2)} \right). \quad (D4)$$

Since the two spin-flip components coincide, we only need

$$\begin{aligned} \langle \Upsilon_{j,m_j}^{(+)} | \sin \theta \cos \theta e^{-i\varphi} \sigma_+ | \Upsilon_{j,m_j}^{(+)} \rangle &= (-1)^{(m_j+1/2)} \sqrt{\frac{2\pi}{15}} \frac{\sqrt{(j+1)^2 - m_j^2}}{j+1} \int d\vartheta Y_{j+1/2}^{-(m_j+1/2)}(\vartheta) Y_2^1(\vartheta) Y_{j+1/2}^{m_j-1/2}(\vartheta) \\ &= (-1)^{(m_j+1/2)} \sqrt{\frac{2}{3}} \sqrt{(j+1)^2 - m_j^2} \begin{pmatrix} j+1/2 & j+1/2 & 2 \\ 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} j+1/2 & j-1/2 & 2 \\ -(m_j+1/2) & m_j-1/2 & 1 \end{pmatrix}, \end{aligned} \quad (D5)$$

where we have used the integration formula of three spherical harmonics in terms of Wigner-3j symbols. The first of the 3j symbols can be trivially calculated, while the use of Regge and permutation symmetries for the second one leads to

$$\langle \Upsilon_{j,m_j}^{(+)} | \sin \theta \cos \theta e^{-i\varphi} \sigma_+ | \Upsilon_{j,m_j}^{(+)} \rangle = \frac{m_j}{2} \frac{(j+1)^2 - m_j^2}{j(j+1)(j+2)}. \quad (D6)$$

Proceeding analogously with the other basis vector $\Upsilon_{j,m_j}^{(-)}$, while combining the spin-conserving and spin-flipping components, we have

$$\mathcal{I}_{j,m_j}^{(\text{d-ame})(\pm)} = \mp \frac{m_j(j + 1/2)}{j(j + 1)}. \quad (D7)$$

The off-diagonal matrix element of $\mathcal{H}^{(ame)}$ restricted to the subspace $S_{n,l,m_j}^e$ is

$$\langle \Phi_{n,j,m_j}^{(-)} | \mathcal{H}^{(ame)} | \Phi_{n,j,m_j}^{(+)} \rangle = \frac{\mu_B B}{4mc^2} I_{n,j \pm 1/2}^{(\text{so})} \mathcal{I}_{j,m_j}^{(\text{od-ame})(\pm)}, \quad (D8)$$
with

$$ \mathcal{I}_{j,m_j}^{(\text{od-ame})} = \langle \Upsilon_{j+1,m_j}^{(-)} | \sin^2 \theta \sigma_z - \frac{\sin \theta \cos \theta}{2} (e^{-i\varphi}\sigma_+ + e^{i\varphi}\sigma_-) | \Upsilon_{j,m_j}^{(+)} \rangle. \quad (\text{D9}) $$

The spin-conserving component of $\mathcal{I}_{j,m_j}^{(\text{od-ame})(+)}$ is given by

$$ \langle \Upsilon_{j+1,m_j}^{(-)} | \sin^2 \theta \sigma_z | \Upsilon_{j,m_j}^{(+)} \rangle = -\frac{1}{2} \frac{\sqrt{(j+1)^2 - m_j^2}}{(j+1)} \left( 1 + \frac{m_j^2}{j(j+2)} \right). \quad (\text{D10}) $$

Since here too the two spin-flip components coincide, we only need

$$ \begin{aligned} \langle \Upsilon_{j+1,m_j}^{(-)} | \sin \theta \cos \theta e^{-i\varphi} \sigma_+ | \Upsilon_{j,m_j}^{(+)} \rangle &= -(-1)^{(m_j+1)/2} \sqrt{\frac{8\pi}{15}} \frac{m_j}{j+1} \int d\vartheta Y_{j+1/2}^{-(m_j+1/2)}(\vartheta) Y_2^1(\vartheta) Y_{j+1/2}^{m_j-1/2}(\vartheta) \\ &= -\frac{m_j^2}{2} \frac{\sqrt{(j+1)^2 - m_j^2}}{j(j+1)(j+2)}, \end{aligned} \quad (\text{D11}) $$

and thus

$$ \mathcal{I}_{j,m_j}^{(\text{od-ame})} = -\frac{1}{2} \frac{\sqrt{(j+1)^2 - m_j^2}}{(j+1)}. \quad (\text{D12}) $$

which leads to the off-diagonal matrix element (78). The simple structure of the angular matrix elements (D7) and (D12) is a consequence of the Wigner–Eckart theorem applied to the diagonal and off-diagonal matrix elements of $\mathcal{H}^{(\text{ame})}$, given by Eqs. (D1) and (D8), respectively.

## Appendix E: Semiclassical description of the spin-orbit coupling

In this appendix we apply the semiclassical description of the spin dynamics in order to calculate the effect of the spin-orbit coupling on the ZFS, showing that the corresponding correction constitutes the semiclassical limit of the quantum perturbative results obtained in Sec. IV F.

### 1. Semiclassics without orbital effects

A trace formula describing the semiclassical DOS for a system with SOC has been developed from the Dirac equation [96], or alternatively, from a coherent-state path-integral description of the spin variables [3, 97–99]. In the so-called weak-coupling limit, in which the SOC is weak enough to provide a negligible perturbation of the orbital motion, and assuming that the applied magnetic field only acts on the spin variables without affecting the orbital part of the wave functions, the oscillating part of the DOS can be written as [100]

$$ \varrho^{\text{osc}}(E, B) = \sum_{\xi} \mathcal{M}_{\xi}^{(s)}(E, B) A_{\xi}(E) \cos \left( \frac{S_{\xi}(E)}{\hbar} - \frac{\pi}{2} \mu_{\xi} \right). \quad (\text{E1}) $$

The sum is over the periodic orbits $\xi$ of the uncoupled system, while $A_{\xi}(E)$, $S_{\xi}(E)$, and $\mu_{\xi}$ are, respectively, the corresponding stability determinant, classical action, and Maslov index. No assumption is made on whether the periodic orbits are isolated or not. The only difference of Eq. (E1) with the standard Berry–Tabor or Gutzwiller DOS is the appearance of the spin modulation factor, that for a spin 1/2 particle takes the form

$$ \mathcal{M}_{\xi}^{(s)}(E, B) = 2 \cos \left( \frac{1}{2} \mathcal{K}_{\xi}(E) \right), \quad (\text{E2}) $$

where $\mathcal{K}_{\xi}(E)$ is the phase accumulated by the precession of the spin under the effect of the Zeeman (16b) and the SOC (24a) terms of the Hamiltonian when the particle travels along the periodic orbit $\xi$. Writing

$$ \mathcal{H}^{(Z)} + \mathcal{H}^{(\text{so})} = \frac{\hbar}{2} \boldsymbol{\sigma} \cdot \mathbf{C}(\mathbf{r}, \mathbf{p}) \quad (\text{E3}) $$

with

$$ \mathbf{C}(\mathbf{r}, \mathbf{p}) = \frac{g_0 \mu_B}{\hbar} \mathbf{B} + \frac{1}{2m^2 c^2} \nabla V \times \mathbf{p}, \quad (\text{E4}) $$

we have, for the case in which $\mathbf{C}$ remains oriented along the $z$-axis, i.e., $\mathbf{C}(\mathbf{r}(t), \mathbf{p}(t)) = C(\mathbf{r}(t), \mathbf{p}(t)) \hat{e}_z$, a simple expression for the accumulated phase [100]

$$ \phi_{\xi}(E) = \int_{0}^{\tau_{\xi}} dt C(\mathbf{r}(t), \mathbf{p}(t)), \quad (\text{E5}) $$

where $\tau_{\xi}$ is the period of the classical orbit.

In the absence of magnetic field and SOC, $\mathcal{M}_{\xi}^{(s)} = 2$ for all $\xi$, and therefore the spin modulation factor merely represents spin degeneracy. Without SOC and under a magnetic field $\mathbf{B} = B \hat{e}_z$, we have

$$ \mathcal{K}_{\xi}(E) = \frac{g_0 \mu_B}{\hbar} B \tau_{\xi}. \quad (\text{E6}) $$
<figure><img src="image_8.png" /><figcaption>FIG. 7. Classical periodic orbit (1, 3) (in red) with $\phi_{\nu\eta}$ the half angle spanned between two collisions and the centripetal vector $\Delta\mathbf{p} = \mathbf{p}_f - \mathbf{p}_i$ being the impulse acquired in the collision.</figcaption></figure>

Using Eqs. (E1) and (12), the resulting ZFS in this particular case can be expressed as

$$ \chi^{(Z)-osc} = -\frac{1}{V} \frac{\partial^2}{\partial B^2} \left[ \sum_{\xi} R(\tau_{\xi}/\tau_T) \left(\frac{\hbar}{\tau_{\xi}}\right)^2 \times \cos\left(\frac{g_0 \mu_B B}{\hbar} \tau_{\xi}\right) \varrho_{\xi}^{osc}(E_F, 0) \right] \Bigg|_{B=0} \quad (E7) $$

where the thermal factor $R$ has been defined in Eq. (13) and $\tau_T = \hbar/\pi k_B T$. In the limit of $T \to 0$, Eq. (E7) simplifies to $\chi^{(Z)-osc} = (\mu_B^2/V)\varrho^{osc}(E_F, 0)$, yielding the small component of the Pauli susceptibility associated with the oscillatory part of the DOS [see the discussion after Eq. (46)].

We address now the accumulated phase, considering the effect of the SOC in the absence of a magnetic field, and focusing ourselves on the case of the sphere. Since the impulse at each bounce is given by the change in momentum $\Delta\mathbf{p}$, for a collision occurring at a time $t_c$ we have

$$ \int_{t_c}^{t_c^+} dt \nabla V(\mathbf{r}(t)) \times \mathbf{p}(t) = -\Delta\mathbf{p} \times \mathbf{p}, \quad (E8) $$

where $\mathbf{p}$ can be taken indistinctly as the momentum before or after the collision, as the angular momentum $\mathbf{L}$ is conserved over the collision (see Fig. 7). Specializing ourselves in the case of a periodic orbit in the equatorial plane defined by the topological indices $(\nu, \eta)$, we have

$$ -\Delta\mathbf{p} \times \mathbf{p} = p^2 \sin(2\varphi_{\nu\eta}) \hat{\mathbf{e}}_z, \quad (E9) $$

where the angle $\phi_{\nu\eta}$, defined in Eq. (A9), corresponds to half the angle spanned between two consecutive bounces (see Fig. 7). Taking into account the $\eta$ bounces of the periodic orbit,

$$ \lambda K_{\nu\eta}(E) = \pm \frac{1}{2} \left(\frac{v}{c}\right)^2 \eta \sin(2\phi_{\nu\eta}), \quad (E10) $$

where the $\pm$ corresponds to the counterclockwise (clockwise) orientation of the trajectory with respect to $\hat{e}_z$ and $v = p/m$. Since the spin modulation factor of Eq. (E2) is obtained as a trace over the spin Hilbert space [100], a change of the spin basis can be implemented without altering the results, and thus the expression (E10) is valid for any periodic orbit of the family defined by the indices $(\nu, \eta)$, not necessarily on the equatorial $z = 0$ plane. We notice that the SOC is irrelevant for the highly symmetric diametral orbits, having $\phi_{\nu\eta} = 0$. Furthermore, we remark that, in analogy with the quantum SOC matrix element (60), the accumulated phase $\lambda K_{\nu\eta}(E)$ is independent of $V_0$, and that it scales as $(v/c)^2$.

The individual trajectories of the family $(\nu, \eta)$ can be labeled by these two topological indices, together with the three Euler angles $(\alpha, \beta, \gamma)$, where $(\alpha, \beta)$ are the polar angles defining the unit vector $\hat{e}_\gamma$ perpendicular to the plane of the motion, and $\gamma$ is the rotation angle in that plane. We choose the above defined angles in such a way that $\beta$ is the angle formed by $\hat{e}_\gamma$ and $\hat{e}_z$ (see Fig. 8).

Comparison between Eqs. (E6) and (E10) shows that the kicked precessions occurring at the $\eta$ bounces of the periodic orbit labeled by $(\nu, \eta, \alpha, \beta, \gamma, (\pm))$ would have on the accumulated phase $\lambda K_{\nu\eta\alpha\beta\gamma(\pm)}(E)$ the same effect than that of a fictitious magnetic field

$$ \mathbf{B}_{\nu\eta\alpha\beta\gamma}^{f(\pm)} = \pm \frac{\hbar}{\mu_B} \frac{v^3}{4ac^2} \cos \varphi_{\nu\eta} \hat{\mathbf{e}}_{\gamma}, \quad (E11) $$

where the $\pm$ corresponds to the orientation of the trajectory with respect to $\hat{e}_\gamma$.

<figure><img src="image_9.png" /><figcaption>FIG. 8. In red: Classical periodic orbit $(\nu, \eta, \alpha, \beta, \gamma, (+))$ belonging to the degenerate family with topological indices $(\nu, \eta) = (1, 5)$, labeled by the Euler angles $(\alpha, \beta, \gamma)$, with a counterclockwise orientation with respect to $\hat{e}_\gamma$. This unit vector, associated with the polar angles $(\alpha, \beta)$, is perpendicular to the plane of motion (whose intersection with the sphere is indicated by the blue circle). We note $\gamma$ the rotation angle setting the trajectory in its plane.</figcaption></figure>
## 2. Semiclassics with orbital effects

We now include both terms of the vector C of Eq. (E4), comprising the effects of a magnetic field and the SOC. Towards simplifying the calculation of the corresponding modulation factor, we consider an effective field

$$ \mathbf{B}_{\nu\eta\alpha\beta\gamma}^{\text{eff}(\pm)} = B\hat{\mathbf{e}}_z + \mathbf{B}_{\nu\eta\alpha\beta\gamma}^{\text{f}(\pm)}, \qquad (\text{E12}) $$

which allows to obtain the spin modulation factor with the help of Eqs. (E2) and (E6) when trading in the last equation B by the magnitude of $\mathbf{B}_{\nu\eta\alpha\beta\gamma}^{\text{eff}(\pm)}$,

$$ \mathbf{B}_{\nu\eta\beta}^{\text{eff}(\pm)} = \left[ B^2 \pm \frac{\hbar}{\mu_B} \frac{v^3 B}{2ac^2} \cos \varphi_{\nu\eta} \cos \beta + \left( \frac{\hbar}{\mu_B} \frac{v^3}{4ac^2} \cos \varphi_{\nu\eta} \right)^2 \right]^{1/2}, \qquad (\text{E13}) $$

which is independent on the angles $\alpha$ and $\gamma$.

As indicated above, Eq. (E1) only considers the effect of the magnetic field on the spin variables. In order to lift up such a restriction and include the perturbation of the classical action of the periodic orbits by a weak magnetic field, we must go one step backwards from Eq. (E1) and write the oscillating part of the DOS as a sum over individual periodic orbits $(\nu, \eta, \alpha, \beta, \gamma(\pm))$, i.e.,

$$ \varrho^{\text{osc}}(E, B) = -\frac{1}{\pi} \text{Im} \left\{ \sum_{\substack{\nu>0 \\ \eta>2\nu}} \sum_{(\pm)} \frac{1}{4\pi^2} \int_0^{2\pi} d\alpha \int_0^{2\pi} d\gamma \int_0^{\pi/2} d\beta \sin\beta \mathcal{M}_{\nu\eta\alpha\beta\gamma}^{(\pm)}(E, B) \mathcal{B}_{\nu\eta}(E) \exp(i\theta_{\nu\eta}(k)) \right\}, \qquad (\text{E14}) $$

where [95, 101]

$$ \mathcal{B}_{\nu\eta}(E) = -i \frac{\sqrt{\pi ka}}{E_0} \frac{(-1)^\nu}{\sqrt{\eta}} \cos \varphi_{\nu\eta} \sin^{3/2} \varphi_{\nu\eta}, \qquad (\text{E15}) $$

and the k-dependent phase $\theta_{\nu\eta}(k)$ has been defined in Eq. (A10). The integral over $\beta$ stops at $\pi/2$ since the angles between $\pi/2$ and $\pi$ correspond to trajectories within the integration interval, but traveled in the clockwise direction. It is easy to recover from Eq. (E14) the particular cases of Eqs. (E1) and (A8), for which the orbit-dependent modulation factor $\mathcal{M}_{\nu\eta\alpha\beta\gamma}^{(\pm)}(E, B)$ leads, upon summing over the individual orbits of the family $(\nu, \eta)$, to the form (E2) and (A11), respectively. In the general case (E4), having a magnetic field and SOC, $\mathcal{M}_{\nu\eta\alpha\beta\gamma}^{(\pm)}(E, B)$ factors out in a spin and an orbital part for the individual trajectories, even if such separation is not possible at the level of families of periodic orbits.

Considering the contribution from the counterclockwise and the clockwise orientations, performing the trivial integrals over the angles $\alpha$ and $\gamma$, and taking into account the perturbation on the classical action of the periodic orbits by the effect of a weak magnetic field, we have

$$ \varrho^{\text{osc}}(E, B) = -\frac{1}{\pi} \text{Im} \left\{ \sum_{\substack{\nu>0 \\ \eta>2\nu}} \int_0^{\pi/2} d\beta \sin\beta \mathcal{B}_{\nu\eta}(E) \left[ \mathcal{M}_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(+)}) \exp\left(i\left[\theta_{\nu\eta}(k) - \frac{2\pi\phi_{\nu\eta}}{\phi_0} \cos\beta\right]\right) \right. \right. \\ \left. \left. + \mathcal{M}_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(-)}) \exp\left(i\left[\theta_{\nu\eta}(k) + \frac{2\pi\phi_{\nu\eta}}{\phi_0} \cos\beta\right]\right) \right] \right\}. \qquad (\text{E16}) $$

The k-dependent phase $\theta_{\nu\eta}$ has been defined in Eq. (A10), while $\pm \phi_{\nu\eta} \cos\beta$ is the flux of the magnetic field $B \hat{\mathbf{e}}_z$ piercing the trajectories $(\nu, \eta, \alpha, \beta, \gamma, (\pm))$, and $\phi_{\nu\eta}$ is given by Eq. (A12). Towards the calculation of the ZFS, it is convenient to recast Eq. (E16) as

$$ \varrho^{\text{osc}}(E, B) = \frac{1}{E_0} \sqrt{\frac{ka}{\pi}} \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \cos \varphi_{\nu\eta} \sin^{3/2} \varphi_{\nu\eta} \left[ I_{\nu\eta,c}(E, B) \cos(\theta_{\nu\eta}(k)) + I_{\nu\eta,s}(E, B) \sin(\theta_{\nu\eta}(k)) \right] \qquad (\text{E17}) $$

with

$$ I_{\nu\eta,c}(E, B) = \int_0^{\pi/2} d\beta \sin\beta \cos\left(\frac{2\pi\phi_{\nu\eta}}{\phi_0} \cos\beta\right) \left[ \mathcal{M}_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(+)}) + \mathcal{M}_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(-)}) \right], \qquad (\text{E18a}) $$

$$ I_{\nu\eta,s}(E, B) = \int_0^{\pi/2} d\beta \sin\beta \sin\left(\frac{2\pi\phi_{\nu\eta}}{\phi_0} \cos\beta\right) \left[ \mathcal{M}_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(+)}) - \mathcal{M}_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(-)}) \right]. \qquad (\text{E18b}) $$
In order to obtain $\chi^{\text{osc}}$ from Eq. (12), we need the second derivative of $\rho^{\text{osc}}$ with respect to $B$, which will be determined by

$$ \left. \frac{\partial^2 I_{\nu\eta,c}}{\partial B^2} \right|_{B=0} = \int_0^{\pi/2} d\beta \sin\beta \left[ -\left(\frac{e}{\hbar c} A_{\nu\eta}\right)^2 \cos^2\beta \left( M_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(+)}) + M_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(-)}) \right) \right]_{B=0} \\ + \left. \frac{\partial^2}{\partial B^2} \left( M_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(+)}) + M_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(-)}) \right) \right|_{B=0}, \qquad (\text{E19a}) $$

$$ \left. \frac{\partial^2 I_{\nu\eta,s}}{\partial B^2} \right|_{B=0} = \frac{2e}{\hbar c} A_{\nu\eta} \int_0^{\pi/2} d\beta \sin\beta \cos\beta \left. \frac{\partial}{\partial B} \left( M_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(+)}) - M_{\nu\eta}^{(s)}(E, B_{\nu\eta\beta}^{\text{eff}(-)}) \right) \right|_{B=0}. \qquad (\text{E19b}) $$

Performing the $\beta$ integrals and neglecting the terms of order $(v/c)^4$ we have

$$ \left. \frac{\partial^2 I_{\nu\eta,c}}{\partial B^2} \right|_{B=0} = -\frac{4}{3} \left(\frac{e}{\hbar c} A_{\nu\eta}\right)^2 - 4 \left(\frac{\mu_B}{\hbar} \frac{L_{\nu\eta}}{v}\right)^2, \qquad (\text{E20a}) $$

$$ \left. \frac{\partial^2 I_{\nu\eta,s}}{\partial B^2} \right|_{B=0} = -\frac{1}{3} \frac{e^2}{\hbar c^4} \frac{v}{ma} A_{\nu\eta} L_{\nu\eta}^2 \cos \varphi_{\nu\eta}. \qquad (\text{E20b}) $$

Using Eqs. (E17) and (12), we verify that the first term on the right-hand side of Eq. (E20a) leads to the expression (A19) of $\chi^{\text{(orb)}-\text{osc}}$, the oscillating part for the orbital component of the ZFS, obtained in Ref. [29] and rederived in Appendix A. Similarly, the second term on the right-hand side of Eq. (E20a) leads to the expression (E7) of the oscillating part of the Pauli spin susceptibility. While Eq. (E20a) simply yields the oscillating part of the nonrelativistic ZFS, its relativistic counterpart follows from Eq. (E20b) and can be written as

$$ \frac{\chi^{(\text{so})-\text{osc}}}{|\chi_L|} = 6 \left(\frac{v_F}{c}\right)^2 \left(\pi k_F a\right)^{1/2} \sum_{\substack{\nu>0 \\ \eta>2\nu}} (-1)^\nu \sqrt{\eta} \sin^{5/2} \varphi_{\nu\eta} \\ \times \cos^3 \varphi_{\nu\eta} R(L_{\nu\eta}/L_T) \sin(\theta_{\nu\eta}(k_F)). \qquad (\text{E21}) $$

Comparing $\chi^{(\text{so})-\text{osc}}$ with $\chi^{(\text{orb})-\text{osc}}$ [cf. Eq. (A19)], we see that the former is smaller than the latter by the factor $(v_F/c)^2$ and due to one less power of $k_F a$. In addition, these two terms are “out of phase”, concerning the $k_F a$ oscillations.

### 3. Equivalence between the semiclassical and quantum perturbative approaches

In the remaining part of this appendix, we show that the previously obtained form (E21) of the spin-orbit correction to the ZFS can also be derived by applying the one-dimensional semiclassical approximation to the quantum perturbative result (92b). For that purpose, we follow the same steps as in Appendix A, addressing the quantum result with the help of the Poisson summation rule and only keeping the highest-order terms in $k_F a$ [see Eqs. (A13)]. The smooth part of $\chi^{(\text{so})}$ resulting from the $\nu = 0$ term has been evaluated in Sec. IV G [see Eq. (93b)]. Using for the oscillating part the same stationary-phase condition than in Appendix A, we have,

to leading order in $k_F a$,

$$ \frac{\chi^{(\text{so})-\text{osc}}}{|\chi_L|} = - \frac{12\sqrt{\pi}}{k_F a} \frac{E_0^2}{mc^2} \int_0^\infty dE f'_{\bar{\mu}_0}(E) \left(\frac{E}{E_0}\right)^{9/4} \\ \times \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \sin^{3/2} \varphi_{\nu\eta} \cos^3 \varphi_{\nu\eta} \cos(\theta_{\nu\eta}(k)). \qquad (\text{E22}) $$

An integration by parts leaves us with the integral of $f'_{\bar{\mu}_0}(E)$ multiplied by a rapidly oscillating function of $E$, that leads to [53]

$$ \frac{\chi^{(\text{so})-\text{osc}}}{|\chi_L|} = \frac{12\sqrt{\pi}}{k_F a} \frac{E_0^2}{mc^2} \sum_{\substack{\nu>0 \\ \eta>2\nu}} \frac{(-1)^\nu}{\sqrt{\eta}} \sin^{3/2} \varphi_{\nu\eta} \\ \times \cos^3 \varphi_{\nu\eta} R(\tau_{\nu\eta}/\tau_T) \\ \times \frac{d}{dE} \left( \left( \left(\frac{E}{E_0}\right)^{9/4} \cos(\theta_{\nu\eta}(k)) \right) \right|_{E=E_F}, \qquad (\text{E23}) $$

which to leading order in $k_F a$ coincides with the semiclassical result (E21). Such an agreement demonstrates that the semiclassical treatment of the spin-orbit interaction [100] is valid not only in the large-spin limit, but also for a finite spin, including spin 1/2.

## Appendix F: Matrix elements for the case of the half-sphere

In this appendix we prove a general result concerning the limiting values of the wave function and its derivative at a potential discontinuity, in one- and three-dimensional geometries, which is crucial in order to obtain the spin-orbit matrix elements (102) and (106). We also use nontrivial recurrence relations allowing to calculate the angular matrix element (103) and the radial one, Eq. (112).

### 1. Wave-function values close to a potential barrier

We consider a one-dimensional Schrödinger equation with a potential barrier of height $V_0$ at $x = a$ (see Fig. 9). We will
<figure><img src="image_10.png" /><figcaption>FIG. 9. Sketch of a one-dimensional sharp potential barrier of height $V_0$ at $x = a$.</figcaption></figure>

prove that its solutions $\psi(x)$ verify

$$ \lim_{V_0 \to \infty} \frac{\sqrt{2mV_0}}{\hbar} \psi(a) = -\psi'(a^{-}). \quad (F1) $$

Without loss of generality, we take $V(x) = 0$ for $x \in (a - \epsilon, a)$ and $V(x) = V_0$ for $x \in (a, a + \epsilon)$, where $\epsilon$ is a small, but finite length. Integrating the Schrödinger equation in the small interval $(a - \epsilon, a + \epsilon)$ yields

$$ -\frac{\hbar^2}{2m} [\psi'(a + \epsilon) - \psi'(a - \epsilon)] + V_0 \int_a^{a+\epsilon} dx \psi(x) = E \int_{a-\epsilon}^{a+\epsilon} dx \psi(x). \quad (F2) $$

Assuming that the barrier extends to $x = +\infty$, we have $\psi(x) = \psi(a) \exp(-\kappa x)$ in the interval $(a, a + \epsilon)$, with $\kappa = \sqrt{2m(V_0 - E)/\hbar}$. Taking the limit of large $V_0$ while keeping a finite value of $\epsilon$ we obtain

$$ \psi'(a + \epsilon) = -\kappa \psi(a) \exp(-\kappa \epsilon) \\ \approx 0, \quad (F3a) $$

$$ V_0 \int_a^{a+\epsilon} dx \psi(x) \approx \frac{V_0}{\kappa} \psi(a). \quad (F3b) $$

Thus, in the leading order in $V_0$, Eq. (F2) can be recast as

$$ \frac{\hbar^2}{2m} \psi'(a - \epsilon) + \sqrt{\frac{V_0 \hbar^2}{2m}} \psi(a) = 0. \quad (F4) $$

Since $\psi'(x)$ is a continuous function in the open interval $(a - \epsilon, a)$, we promptly obtain the announced result (F1). The particular case of the radial wave function is contained in Eq. (C7) and has been used in obtaining Eq. (59).

The calculation of the spin-orbit matrix elements in the case of the HS necessitates the generalization of (F1) to a two-dimensional potential barrier. Since the potential (97) defining the confinement in the HS with a finite height $V_0$ has axial symmetry, $m_z$ is a good quantum number, and the eigenfunctions have the form $\psi_{m_z}(r, \theta, \varphi) = \mathbf{f}_{m_z}(r, \theta) \exp(im_z\varphi)$, with $\mathbf{f}_{m_z}$ verifying

$$ -\frac{\hbar^2}{2m} \left( \frac{1}{r^2} \partial_r r^2 \partial_r + \frac{1}{r^2 \sin \theta} \partial_\theta \sin \theta \partial_\theta - \frac{m_z^2}{r^2 \sin^2 \theta} \right) \mathbf{f}_{m_z}(r, \theta) + V(r, \theta) \mathbf{f}_{m_z}(r, \theta) = E \mathbf{f}_{m_z}(r, \theta). \quad (F5) $$

Marching on the footprints of the derivation for the one-dimensional case, we integrate the previous equation in the small angular interval $(\pi/2 - \epsilon, \pi/2 + \epsilon)$, obtaining

$$ -\frac{\hbar^2}{2m} \left( \left. \frac{2\epsilon}{r^2} \partial_r r^2 \partial_r \mathbf{f}_{m_z}(r, \theta) \right|_{\theta=\pi/2} + \left. \frac{1}{r^2} \partial_\theta \mathbf{f}_{m_z}(r, \theta) \right|_{\theta=\pi/2-\epsilon} \left. - \frac{2\epsilon m_z^2}{r^2} \mathbf{f}_{m_z}(r, \theta) \right|_{\theta=\pi/2} \right) + V_0 \int_{\pi/2}^{\pi/2+\epsilon} d\theta \mathbf{f}_{m_z}(r, \theta) \\ = 2\epsilon E \mathbf{f}_{m_z}(r, \theta) \bigg|_{\theta=\pi/2}. \quad (F6) $$

In the leading order in $V_0$ we have $\mathbf{f}_{m_z}(r, \theta) = \mathbf{f}_{m_z}(r, \theta=\pi/2) \exp(-\sqrt{2mV_0}r^2[\theta - \pi/2]/\hbar)$ in the interval $(\pi/2, \pi/2 + \epsilon)$. Thus, taking the limit of large $V_0$ and then that of $\epsilon \to 0$ we obtain

$$ \frac{r\sqrt{2mV_0}}{\hbar} \mathbf{f}_{m_z}(r, \theta=\pi/2) \approx -\partial_\theta \mathbf{f}_{m_z}(r, \theta) \bigg|_{\theta=\pi/2-}, \quad (F7) $$

which generalizes Eq. (F1) to the case of a two-dimensional potential barrier with axial symmetry.

Taking the hard-wall limit of $V_0 \to \infty$ for the potential (97), the solutions $\psi_{m_z}^{(V_0, n, l)}(r, \theta, \varphi)$ converge towards the orbital wave functions (98), with the condition of $l + m_z$ being odd. Thus, we can write

$$ \frac{\sqrt{mV_0}}{\hbar} \psi_{m_z}^{(V_0, n, l)}(r, \theta=\pi/2, \varphi) \approx -\sqrt{l(l+1)-m_z(m_z+1)} \frac{R_{nl}(r)}{r} e^{-i\varphi} Y_l^{m_z+1}(\theta=\pi/2, \varphi), \quad (F8) $$
where we have used

$$ \frac{\partial Y_l^{m_z}(\theta, \varphi)}{\partial \theta} = m_z \cot(\theta) Y_l^{m_z}(\theta, \varphi) + \sqrt{(l-m_z)(l+m_z+1)} e^{-i\varphi} Y_l^{m_z+1}(\theta, \varphi), \quad (F9) $$

taking $\theta = \pi/2$.

## 2. Angular integral of the matrix element (103)

The angular integral (104) appearing in the nondiagonal matrix element (103) of $\mathcal{H}^{(\text{so})-(\text{HS})}$ can be expressed as

$$ \mathcal{I}_{l',l,m'_z}^{(\text{dome})} = \frac{1}{2} \sqrt{(2l+1)(2l'+1)} \sqrt{\frac{(l-m)!(l'-m'_z)!}{(l+m'_z)!(l'+m'_z)!}} \mathcal{Y}_{l',l}^{m'_z}, \quad (F10) $$

where we have used the definition of the spherical harmonic $Y_l^{m_z}(\vartheta)$ in terms of the associated Legendre function $P_l^{m_z}(\cos\theta)$ [with the standard convention of assigning the Condon-Shortley phase $(-1)^{m_z}$ to the latter], and introduced

$$ \mathcal{Y}_{l',l}^{m'_z} = \int_0^1 dx P_{l'}^{m'_z}(x) P_l^{m'_z}(x). \quad (F11) $$

In the case that interests us, with $l' + m'_z$ odd and $l + m'_z$ even, we can prove, directly from the differential equations fulfilled by the associated Legendre functions $P_l^{m'_z}$ and $P_{l'}^{m'_z}$ with $l \neq l'$, the following useful identity:

$$ \mathcal{Y}_{l',l}^{m'_z} = -\frac{l' + 1 - m'_z}{l'(l' + 1) - l(l + 1)} P_l^{m'_z}(0) P_{l'+1}^{m'_z}(0), \quad (F12) $$

which, together with the relationship

$$ P_l^{m'_z}(0) = \frac{(-1)^{(l+m'_z)/2}}{2^l} \frac{(l+m'_z)!}{\left(\frac{l+m'_z}{2}\right)! \left(\frac{l-m'_z}{2}\right)!}, \quad (F13) $$

valid for even $l + m'_z$, lead to

$$ \mathcal{I}_{l',l,m}^{(\text{dome})} = \frac{(-1)^{m+(l+l'+1)/2}}{2^{l+l'}} \frac{\sqrt{(2l+1)(2l'+1)}}{l(l+1)-l'(l'+1)} \\ \times \frac{\sqrt{(l+m)!(l-m)!}}{\left(\frac{l+m}{2}\right)!\left(\frac{l-m}{2}\right)!} \frac{\sqrt{(l'+m)!(l'-m)!}}{\left(\frac{l'+m-1}{2}\right)!\left(\frac{l'-m-1}{2}\right)!}, \quad (F14) $$

allowing to give a closed expression to the nondiagonal matrix element (103) of $\mathcal{H}^{(\text{dome})}$.

## 3. Definite integrals of two spherical Bessel functions

The radial matrix elements (112), appearing in the perturbative treatment of the magnetic field for the HS of Sec. VC, as well as in the diagonalization to obtain the finite-field spectrum of the sphere [29], can be expressed as

$$ \mathcal{R}_{n',l',n,l} = \frac{2 L_{l',l}^{(4)}(\zeta_{n',l'}, \zeta_{n,l})}{|j_{l'+1}(\zeta_{n',l'}) j_{l+1}(\zeta_{n,l})|}, \quad (F15) $$

with

$$ L_{l',l}^{(q)}(\alpha, \beta) = \int_0^1 d\zeta \, \zeta^q j_{l'}(\alpha\zeta) j_l(\beta\zeta), \quad (F16) $$

and where $\alpha$ and $\beta$ are such that $j_{l-2}(\alpha) = 0$ and $j_l(\beta) = 0$, respectively.

The particular case of $\mathcal{R}_{n',l,n,l}$ appearing in the first term of Eq. (110b) requires a simpler term

$$ K_l^{(q)}(\alpha, \beta) = L_{l,l}^{(q)}(\alpha, \beta), \quad (F17) $$

invoking the definite integral of two spherical Bessel functions with the same order. Using repetitively the recursive formulas developed for the indefinite integrals of spherical Bessel functions of the same order [102], we have

$$ K_l^{(4)}(\alpha, \beta) = \left(\frac{\alpha^2 + \beta^2}{2\alpha\beta}\right)^l K_0^{(4)}(\alpha, \beta) - \sum_{d=1}^{l} \frac{(\alpha^2 + \beta^2)^{d-1}}{(2\alpha\beta)^d} \left\{ 2j_{l-d}(\alpha) j_{l-d}(\beta) \right. \\ \left. - \alpha \left( \frac{2[2(l-d)+3]}{\alpha^2 - \beta^2} - 1 \right) j_{l+1-d}(\alpha) j_{l-d}(\beta) + \beta \left( \frac{2[2(l-d)+3]}{\alpha^2 - \beta^2} + 1 \right) j_{l-d}(\alpha) j_{l+1-d}(\beta) \right\} \quad (F18) $$
with

$$ K_0^{(4)}(\alpha, \beta) = \frac{1}{2\alpha\beta} \left\{ \frac{1}{(\alpha - \beta)^3} \left[ 2(\alpha - \beta) \cos(\alpha - \beta) + ((\alpha - \beta)^2 - 2) \sin(\alpha - \beta) \right] - \frac{1}{(\alpha + \beta)^3} \left[ 2(\alpha + \beta) \cos(\alpha + \beta) + ((\alpha + \beta)^2 - 2) \sin(\alpha + \beta) \right] \right\}. \quad (F19) $$

The recursive formulas developed for the indefinite integrals of spherical Bessel functions of different orders [102] allow us to write

$$ L_{l-2,l}^{(4)}(\alpha, \beta) = \frac{2l-1}{2} \frac{\alpha}{\beta} K_{l-1}^{(4)}(\alpha, \beta) - \frac{2l+1}{2} K_{l-2}^{(4)}(\alpha, \beta), \quad (F20) $$

which can be directly evaluated from Eq. (F18), and also expressed as

$$ L_{l-2,l}^{(4)}(\alpha, \beta) = \frac{1}{4} \left[ (2l-1) \left(\frac{\alpha}{\beta}\right)^2 - 2l-3 \right] \left[ \left(\frac{\alpha^2 + \beta^2}{2\alpha\beta}\right)^{l-2} K_0^{(4)}(\alpha, \beta) - \sum_{d=1}^{l-2} \frac{(\alpha^2 + \beta^2)^{d-1}}{(2\alpha\beta)^d} \left\{ 2j_{l-2-d}(\alpha) j_{l-2-d}(\beta) - \alpha \left( \frac{2[2(l-d)-1]}{\alpha^2 - \beta^2} - 1 \right) j_{l-1-d}(\alpha) j_{l-2-d}(\beta) + \beta \left( \frac{2[2(l-d)-1]}{\alpha^2 - \beta^2} + 1 \right) j_{l-2-d}(\alpha) j_{l-1-d}(\beta) \right\} \right. \\ \left. + \frac{2l-1}{4} \frac{\alpha}{\beta^2} \left( \frac{2[2l-1]}{\alpha^2 - \beta^2} - 1 \right) j_{l-1}(\alpha) j_{l-2}(\beta). \quad (F21) $$

The integral $L_{l+2,l}^{(4)}(\alpha, \beta)$ can be obtained from the expression of $L_{l-2,l}^{(4)}(\alpha, \beta)$ above by implementing the shift of $l$ to $l + 2$ while exchanging the role of $\alpha$ and $\beta$.

[1] G. Bergmann, Phys. Rep. 107, 1 (1984).

[2] J. B. Miller, D. M. Zumbühl, C. M. Marcus, Y. B. Lyanda-Geller, D. Goldhaber-Gordon, K. Campman, and A. C. Gossard, Phys. Rev. Lett. 90, 076807 (2003).

[3] O. Zaitsev, D. Frustaglia, and K. Richter, Phys. Rev. B 72, 155325 (2005).

[4] N. Nagaosa, J. Sinova, S. Onoda, A. H. MacDonald, and N. P. Ong, Rev. Mod. Phys. 82, 1539 (2010).

[5] I. Žutić, J. Fabian, and S. Das Sarma, Rev. Mod. Phys. 76, 323 (2004).

[6] H.-A. Engel, E. Rashba, and B. Halperin, in Handbook of Magnetism and Advanced Magnetic Materials, Vol. 5, H. Krohnmüller and S. Parkin (eds.) (John Wiley & Sons Ltd, Chichester, 2007).

[7] G. A. Intronati, P. I. Tamborenea, D. Weinmann, and R. A. Jalabert, Phys. Rev. Lett. 108, 016601 (2012).

[8] N. Taniguchi, A. Hashimoto, B. D. Simons, and B. L. Altshuler, Europhys. Lett. 27, 335 (1994).

[9] Y. Alhassid, Rev. Mod. Phys. 72, 895 (2000).

[10] D. G. Salinas, S. Guéron, D. C. Ralph, C. T. Black, and M. Tinkham, Phys. Rev. B 60, 6137 (1999).

[11] D. Davidović and M. Tinkham, Phys. Rev. Lett. 83, 1644 (1999).

[12] A. Kawabata, J. Phys. Soc. Jpn. 29, 902 (1970).

[13] W. P. Halperin, Rev. Mod. Phys. 58, 533 (1986).

[14] P. W. Brouwer, X. Waintal, and B. I. Halperin, Phys. Rev. Lett. 85, 369 (2000).

[15] K. A. Matveev, L. I. Glazman, and A. I. Larkin, Phys. Rev. Lett. 85, 2789 (2000).

[16] H. Bouchiat and G. Montambaux, J. Phys. (France) 50, 2695 (1989).

[17] Y. Imry, in Coherence Effects in Condensed Matter Systems, ed. by B. Kramer (Plenum, New York, 1991).

[18] P. G. van Rhee, P. Zijlstra, T. G. A. Verhagen, J. Aarts, M. I. Katsnelson, J. C. Maan, M. Orrit, and P. C. M. Christianen, Phys. Rev. Lett. 111, 127202 (2013).

[19] B. Murzaliev, M. Titov, and M. I. Katsnelson, Phys. Rev. B 100, 075426 (2019).

[20] H. Hori, T. Teranishi, Y. Nakae, Y. Seino, M. Miyake, and S. Yamada, Phys. Lett. A 263, 406 (1999).

[21] Y. Nakae, Y. Seino, T. Teranishi, M. Miyake, S. Yamada, and H. Hori, Physica B 284-288, 1758 (2000).

[22] H. Hori, Y. Yamamoto, T. Iwamoto, T. Miura, T. Teranishi, and M. Miyake, Phys. Rev. B 69, 174411 (2004).

[23] Y. Yamamoto, T. Miura, M. Suzuki, N. Kawamura, H. Miyagawa, T. Nakamura, K. Kobayashi, T. Teranishi, and H. Hori, Phys. Rev. Lett. 93, 116801 (2004).

[24] Y. Yamamoto and H. Hori, Rev. Adv. Mater. Sci. 12, 23 (2006).

[25] E. Guerrero, M. A. Muñoz-Márquez, E. Fernández-Pinel, P. Crespo, A. Hernando, and A. Ferández, J. Nanopart. Res. 10, 179 (2008).

[26] E. Guerrero, M. A. Muñoz-Márquez, M. A. García, P. Crespo, E. Fernández-Pinel, A. Hernando, and A. Fernández, Nanotechnology 19, 175701 (2008).

[27] J. Bartolomé, F. Bartolomé, L. M. García, A. I. Figueroa, A. Repollés, M. J. Martínez-Pérez, F. Luis, C. Magén, S. Selenska-Pobell, F. Pobell, T. Reitz, R. Schönemann,
T. Herrmannsdörfer, M. Merroun, A. Geissler, F. Wilhelm, and A. Rogalev, Phys. Rev. Lett. 109, 247203 (2012).

[28] M. Agrachev, S. Antonello, T. Dainese, M. Ruzzi, A. Zoleo, E. Aprà, N. Govind, A. Fortunelli, L. Sementa, and F. Maran, ACS Omega 2, 2607 (2017).

[29] M. Gómez Viloria, G. Weick, D. Weinmann, and R. A. Jalabert, Phys. Rev. B 98, 195417 (2018).

[30] P. Crespo, R. Litrán, T. C. Rojas, M. Multigner, J. M. de la Fuente, J. C. Sánchez-López, M. A. García, A. Hernando, S. Penadés, and A. Fernández, Phys. Rev. Lett. 93, 087204 (2004).

[31] P. Dutta, S. Pal, M. S. Seehra, M. Anand, and C. B. Roberts, Appl. Phys. Lett. 90, 213102 (2007).

[32] Y. Imry, Phys. Rev. B 91, 104503 (2015).

[33] P. Crespo, M. A. García, E. F. Pinel, M. Multigner, D. Alcántara, J. M. de la Fuente, S. Penadés, and A. Hernando, Phys. Rev. Lett. 97, 177203 (2006).

[34] B. Donnio, P. García-Vázquez, J.-L. Gallani, D. Guillon, and E. Terazzi, Adv. Mater. 19, 3534 (2007).

[35] J. S. Garitaonandia, M. Insausti, E. Goikolea, M. Suzuki, J. D. Cashion, N. Kawamura, H. Ohsawa, I. G. de Muro, K. Suzuki, F. Plazaola, and T. Rojo, Nano Lett. 8, 661 (2008).

[36] J. de la Venta, V. Bouzas, A. Pucci, M. A. Laguna-Marco, D. Haskel, S. G. E. te Velthuis, A. Hoffmann, J. Lal, M. Bleuel, G. Ruggeri, C. de Julián Fernández, and M. A. García, J. Nanosci. Nanotechnol. 9, 6434 (2009).

[37] B. Donnio, A. Derory, E. Terazzi, M. Drillon, D. Guillon, and J.-L. Gallani, Soft Matter 6, 965 (2010).

[38] U. Maitra, B. Das, N. Kumar, A. Sundaresan, and C. N. R. Rao, ChemPhysChem 12, 2322 (2011).

[39] R. Gréget, G. L. Nealon, B. Vileno, P. Turek, C. Mény, F. Ott, A. Derory, E. Voirin, E. Rivière, A. Rogalev, F. Wilhelm, L. Joly, W. Knafo, G. Ballon, E. Terazzi, J.-P. Kappler, B. Donnio, and J.-L. Gallani, ChemPhysChem 13, 3092 (2012).

[40] G. Percebois, D. Weinmann, R. A. Jalabert, and G. Weick, arXiv:2107.04508.

[41] G. L. Nealon, B. Donnio, R. Gréget, J.-P. Kappler, E. Terazzi, and J.-L. Gallani, Nanoscale 4, 5244 (2012).

[42] T. Rangel, D. Kecik, P. E. Trevisanutto, G.-M. Rignanese, H. Van Swygenhoven, and V. Olevano, Phys. Rev. B 86, 125125 (2012).

[43] P. Monod and A. Janossy, J. Low Temp. Phys. 26, 311 (1977).

[44] M. Isasa, E. Villamor, L. E. Hueso, M. Gradhand, and F. Casanova, Phys. Rev. B 91, 024402 (2015).

[45] F. Schwabl, Advanced Quantum Mechanics (Springer-Verlag, Berlin, 1997).

[46] L. D. Landau and E. M. Lifshitz, Statistical Physics (Pergamon, Oxford, 1985).

[47] M. Suzuki, N. Kawamura, H. Miyagawa, J. S. Garitaonandia, Y. Yamamoto, and H. Hori, Phys. Rev. Lett. 108, 047201 (2012).

[48] A. Rukhadze and V. Silin, J. Exp. Theor. Phys. 11, 463 (1960).

[49] E. M. Chudnovsky, J. Phys. A: Math. Gen. 14, 2091 (1981).

[50] J. M. van Ruitenbeek and D. A. van Leeuwen, Phys. Rev. Lett. 67, 640 (1991).

[51] D. A. van Leeuwen, Magnetic Moments in Metacluster Molecules (PhD thesis, University of Leiden, unpublished, 1993).

[52] H.-J. van Leeuwen, J. Phys. Radium 2, 361 (1921).

[53] K. Richter, D. Ullmo, and R. A. Jalabert, Phys. Rep. 276, 1 (1996).

[54] D. Ullmo, K. Richter, and R. A. Jalabert, Phys. Rev. Lett. 74, 383 (1995).

[55] M. C. Gutzwiller, Chaos in Classical and Quantum Mechanics (Springer-Verlag, Berlin, 1990).

[56] A. Schmid, Phys. Rev. Lett. 66, 80 (1991).

[57] F. von Oppen and E. K. Riedel, Phys. Rev. Lett. 66, 84 (1991).

[58] B. L. Altshuler, Y. Gefen, and Y. Imry, Phys. Rev. Lett. 66, 88 (1991).

[59] V. Ambegaokar and U. Eckern, Phys. Rev. Lett. 65, 381 (1990).

[60] D. Ullmo, H. U. Baranger, K. Richter, F. von Oppen, and R. A. Jalabert, Phys. Rev. Lett. 80, 895 (1998).

[61] L. Lévy, D. Reich, L. Pfeiffer, and K. West, Physica B: Condens. Matter 189, 204 (1993).

[62] R. Deblock, R. Bel, B. Reulet, H. Bouchiat, and D. Mailly, Phys. Rev. Lett. 89, 206803 (2002).

[63] H. Bary-Soroker, O. Entin-Wohlman, and Y. Imry, Phys. Rev. Lett. 101, 057001 (2008).

[64] W. E. Shanks, Persistent Currents in Normal Metal Rings (PhD thesis, Yale University, 2011).

[65] H. Mathur and A. D. Stone, Phys. Rev. B 44, 10957(R) (1991).

[66] A. G. Aronov and Y. B. Lyanda-Geller, Phys. Rev. Lett. 70, 343 (1993).

[67] J. Splettstoesser, M. Governale, and U. Zülicke, Phys. Rev. B 68, 165341 (2003).

[68] L. L. Foldy and S. A. Wouthuysen, Phys. Rev. 78, 29 (1950).

[69] R. Mondal, M. Berritta, C. Paillard, S. Singh, B. Dkhil, P. M. Oppeneer, and L. Bellaiche, Phys. Rev. B 92, 100402 (2015).

[70] G. Weick, R. A. Molina, D. Weinmann, and R. A. Jalabert, Phys. Rev. B 72, 115410 (2005).

[71] G. Weick, G.-L. Ingold, R. A. Jalabert, and D. Weinmann, Phys. Rev. B 74, 165421 (2006).

[72] M. Brack, Rev. Mod. Phys. 65, 677 (1993).

[73] R. A. Molina, D. Weinmann, and R. A. Jalabert, Phys. Rev. B 65, 155427 (2002).

[74] W.-D. Schöne, W. Ekardt, and J. M. Pacheco, Phys. Rev. B 50, 11079 (1994).

[75] J. P. Bucher, P. Xia, and L. A. Bloomfield, Phys. Rev. B 42, 10858 (1990).

[76] N. Pavloff and M. Hansen, Z. Phys. D 24, 57 (1992).

[77] K. Richter, D. Ullmo, and R. A. Jalabert, Phys. Rev. B 54, R5219 (1996).

[78] K. Richter, D. Ullmo, and R. A. Jalabert, J. Math. Phys. 37, 5087 (1996).

[79] J. M. van Ruitenbeek and D. A. van Leeuwen, Mod. Phys. Lett. B 07, 1053 (1993).

[80] N. Ashcroft and N. Mermin, Solid State Physics (Sounders College, Philadelphia, 1976).

[81] F. von Oppen, Phys. Rev. B 50, 17151 (1994).

[82] S. Frauendorf, V. M. Kolomietz, A. G. Magner, and A. I. Sanzhur, Phys. Rev. B 58, 5622 (1998).

[83] W. Greiner, Relativistic Quantum Mechanics, Vol. 3 (Springer, 1990).

[84] P. Alberto, M. Fiolhais, and M. Oliveira, Eur. J. Phys. 19, 553 (1998).

[85] A. Mikhailov and V. Polikanov, J. Exp. Theor. Phys. 27, 95 (1968).

[86] C. Yannouleas and R. A. Broglia, Ann. Phys. (N.Y.) 217, 105 (1992).

[87] C.-H. Chang, Phys. Rev. E 67, 046201 (2003).

[88] M. Gómez Viloria, unpublished.

[89] M. Gómez Viloria, Effect of Spin-Orbit Coupling on the Magnetic Response of Confined Systems (PhD thesis, University of Strasbourg, 2021).

[90] K. Richter and B. Mehlig, Europhys. Lett. 41, 587 (1998).

[91] H. Bluhm, N. C. Koshnick, J. A. Bert, M. E. Huber, and K. A. Moler, Phys. Rev. Lett. 102, 136802 (2009).
[92] A. C. Bleszynski-Jayich, W. E. Shanks, B. Peaudecerf, E.GNOSSEL, F. von Oppen, L. Glazman, and J. G. E. Harris, *Science* **326**, 272 (2009).

[93] G. Weick, *Quantum Dissipation and Decoherence of Collective Excitations in Metallic Nanoparticles* (PhD thesis, Université Louis Pasteur and Universität Augsburg, 2006).

[94] R. Balian and C. Bloch, *Ann. Phys. (N.Y.)* **69**, 76 (1972).

[95] K. Tanaka, S. C. Creagh, and M. Brack, *Phys. Rev. B* **53**, 16050 (1996).

[96] J. Bolte and S. Keppeler, *Phys. Rev. Lett.* **81**, 1987 (1998).

[97] A. Alscher and H. Grabert, *J. Phys. A: Math. Gen.* **32**, 4907 (1999).

[98] M. Pletyukhov and O. Zaitsev, *J. Phys. A: Math. Gen.* **36**, 5181 (2003).

[99] A. D. Ribeiro, M. A. M. de Aguiar, and A. F. R. de Toledo Piza, *J. Math. Phys.* **48**, 112103 (2007).

[100] O. Zaitsev, *J. Phys. A: Math. Gen.* **35**, L721 (2002).

[101] J. B. Keller and S. Rubinow, *Ann. Phys. (N.Y.)* **9**, 24 (1960).

[102] J. Bloomfield, S. Face, and Z. Moss, arXiv:1703.06428.