On the Origins of Spontaneous Spherical Symmetry-Breaking in Open-Shell Atoms Through Polymer Self-Consistent Field Theory

Two features of atomic systems that emerge from quantum theory are the presence of highly inhomogeneous shell structure and the spontaneous breaking of spherical symmetry. These two attributes play an important role in determining atomic size and facilitating chemical bonding. An intuitive understanding of the mechanism behind them would be helpful in considerations of molecules and bulk materials. On one hand, quantum superposition arguments can be made in support of the point of view that all isolated atoms, even open shell atoms, are spherically symmetric [9, 13]. On the other hand, the concept of an “isolated” atom is purely theoretical, and operationally, it is nowadays accepted that open-shell atoms used as building blocks for molecular calculations are not generally treated as have spherically symmetric electron densities [11]. In particular, the recent density functional theory (DFT) study of open-shell atoms by Chowdhury and Perdew [8] has examined the implications of spherical symmetry-breaking for molecular bonding and atomization energies.

An alternative to the Kohn-Sham (KS) DFT used by Chowdhury and Perdew [8] is polymer self-consistent field theory (SCFT). Based on a quantum-classical isomorphism introduced by Feynman [12, 7], classical statistical mechanics is used to represent quantum particles as ring “polymers”, that is, extended non-point like objects, embedded in an extra thermal dimension [39]. SCFT has been shown to be equivalent to KS-DFT, assuming the Pauli-exclusion principle holds [39], and incorporates all quantum effects [41] in terms of classical quantities which enables explanations of quantum phenomena more in tune with classical intuition. The use of a Pauli potential and lack of explicit orbitals means that SCFT is related to orbital-free (OF) DFT, and as is typical in OF-DFT, the exclusion principle can be incorporated using approximations [40]. SCFT predicts molecular bonding [35], atomic shell structure [24], includes temperature dependence [39], and can be related to dynamical quantum mechanics [41]. SCFT also has foundational implications due to its classical ensemble formulation [41]. Like DFT, SCFT is based on an energy functional, and since the quantities are real and classical, it can be decomposed into thermodynamic components in order to explain the origins of predicted structures [29].

The purpose of this paper is to show that SCFT spontaneously predicts shell structure and spherical symmetry-breaking in isolated atoms, and to give the thermodynamic explanations for non-spherical ground-state structures. We partition the SCFT free energy functional for the electron density into components: the translational entropy contribution, polymer configurational entropy contribution (equivalent to the quantum kinetic energy), and the internal energy which includes electron-electron, electron-nucleus, and Pauli-exclusion interactions. We find that for the atoms we studied from hydrogen to neon, all thermodynamic contributions cause the free energy to increase for spherical symmetry-breaking except the electron-nucleus potential. In other words, the overall energy is reduced when the electrons can move closer to the nucleus by breaking spherical symmetry and this more than compensates for all other factors which would tend to increase the energy.

Neutral atoms in their ground-state were studied in a previous work [24] using a spherical-averaging approximation, which amounts to representing only the radial part of the electron density, as well as the adoption of partial atomic shell information. These approximations were useful to see how well the model could replicate atomic trends for a large range of elements, most importantly if it could produce the proper shell structure. We will generalize that work by restoring the full angular distribution of the electron density. We organize the paper as follows. In section II, the conceptual basis of the theory is outlined and a new derivation of the model equations is given. The fields in the model are introduced, the free energy components are given, and the section ends with a discussion about the numerical methods used. The atomic binding energies for the elements hydrogen to neon and the angular electron density contour plots for the elements boron to neon are presented in section III, along with tables of data for the free energy components and various density constraints that the model satisfies with a proof of the equivalence with KS-DFT. The Discussion in section IV demonstrates that the external potential is responsible for the spherical symmetry breaking. The paper concludes in section V and some future research directions are discussed.

The pioneering work of Feynman in developing the path integral formulation of quantum theory led to the deduction of the isomorphism between quantum mechanics and the statistical mechanics of ring-like molecules, using a Wick rotation of time $t=-i\hbar\beta$ to a parameter $\beta$ [12]. The Wick rotation essentially allows one to transform a dynamics problem into a statics problem in one higher dimension, where temporal variables are typically exchanged for spatial variables [4]. Part of the intuition behind why the Wick rotation takes quantum mechanics to the statistical mechanics of ring-like molecules stems from the nature of the trace operator in the partition function, which is a sum of the Boltzmann factors over configurations starting and ending at the same position. This cyclic aspect is crucial to the interpretation, especially after having transferred to the path integral picture, because the path integral allows one to see that each of the trajectories followed through the imaginary time configuration parameter space correspond to the different possible arrangements of quantum particles comprising a system of quantum particles [12, 7]. The correspondence with the mathematics of polyatomic fluids noted by Chandler [7] lends itself to label each of the quantum particles in the system as atoms comprising a molecule. The probabilistic uncertainty in the position and momenta of the quantum system are then manifested as the collection of system arrangements associated to a given configuration. The connection to polymer SCFT is then made clear through the insights of Matsen [28]. The periodicity of the imaginary time parameter $\beta$ also happens to be a part of one of the conditions to be a Kubo-Martin-Schwinger (KMS) state, which is a general notion of being in a thermal state [14]. Although the mathematics is essentially identical between quantum mechanical and polymeric systems, the interpretation of the fundamental constituents and the mechanisms that govern their behaviour has changed dramatically. In the quantum mechanical case, the time evolution of the system could trace out many different paths, forming a distribution of them that expresses the uncertainty in which path will be followed. In the statistical mechanical case, the system could explore many different energy configurations, forming a distribution of them where thermal fluctuations represent the uncertainty in which configuration will be chosen. A set of postulates has been proposed to bridge between finite-temperature quantum DFT and ring polymer SCFT [41].

In the static case considered here, there are two postulates needed, in addition to those from statistical mechanics, to describe quantum many-body systems of Fermions: 1) “pairs” (i.e. 1-2) of quantum particles are classically modelled as stochastic chains in a 4-dimensional thermal space and 2) the chains have excluded-volume with respect to 4-dimensional thermal space. The first postulate can be seen to be equivalent to the Heisenberg uncertainty principle [39, 40], and the introduction of the term “pairs” is to account for Fermion spin, while the second postulate is a (conjectured) statement of the Pauli-exclusion principle. Similar to ring polymer SCFT, the base entities in this theory are fundamentally distinguishable from each other, so the notion of particle “pairs” with higher-dimensional excluded-volume is the mechanism that emulates particle statistics from quantum statistical mechanics. To recover the standard electron density, field, and electron number that appear in DFT, we simply sum up all of the pairs, which we denote with a Greek index. The two postulates in combination with the usual approximations used in the study of many-body systems (i.e. Born-Oppenheimer approximation, point-particle nuclei) are then enough to derive the expression for the partition function of the system. Working in the canonical ensemble for simplicity, with $N$-body potential $U$, the $N$-body partition function $Q_{N}$ can be expressed as a path integral in configuration space as

$\displaystyle Q_{N}(\beta)=\prod^{N}_{i}\int\mathrm{d}\bm{r}_{i}\int\mathcal{D}[\bm{r}_{i}]e^{-\int^{\beta}_{0}\mathrm{d}\tau\left[\sum_{j}\frac{m}{2\hbar^{2}}\left|\frac{\mathrm{d}\bm{r}_{j}(\tau)}{\mathrm{d}\tau}\right|^{2}+U[\hat{n}](\bm{r}_{1}(\tau),\ldots,\bm{r}_{N}(\tau))\right]}$

(1)

where the parameter $\beta$ is the Lagrange multiplier that ensures the expectation value of the free energy remains constant, which will end up being the reciprocal thermal energy $\beta=1/k_{B}T$ (where $k_{B}$ is Boltzmann’s constant and $T$ is the temperature) due to the KMS condition [14]; $\bm{r}_{i}(s)$ is the parametrized curve representing the $i$th quantum particle as depicted in figure 1, with the parameter $s$ running from $0$ (a high classical temperature) to $\beta$ (a lower temperature); and the $N$-body potential $U$ is expressed as a functional of an electron density operator $\hat{n}(\bm{r})$, which is defined to be

$\displaystyle\hat{n}(\bm{r})=\sum_{\mu}\hat{n}_{\mu}(\bm{r})=\sum_{\mu}\sum^{N_{\mu}}_{i}\delta(\bm{r}-\bm{r}_{i})\,.$

(2)

$U[\hat{n}]$ can be re-expressed in terms of the fields $\mathcal{N}_{\mu}(\bm{r})$ using the functional Dirac delta $\delta[\mathcal{N}_{\mu}-\hat{n}_{\mu}]$, which can in-turn be expressed in terms of its functional Fourier transform representation with respect to conjugate fields $W_{\mu}(\bm{r})$, allowing eqn. 1 to be expressed as

	$\displaystyle Q_{N}(\beta)$	$\displaystyle=\prod_{\mu}\int\int\mathcal{D}[\mathcal{N}_{\mu}]\mathcal{D}[W_{\mu}]e^{-\beta U[\mathcal{N}]+\beta\int\mathrm{d}\bm{r}^{\prime}W_{\mu}(\bm{r}^{\prime})\mathcal{N}_{\mu}(\bm{r}^{\prime})}\prod^{N_{\mu}}_{i}\int\mathrm{d}\bm{r}_{i}$
		$\displaystyle\qquad\qquad\times\int\mathcal{D}[\bm{r}_{i}]e^{-\int^{\beta}_{0}\mathrm{d}\tau\left[\frac{m}{2\hbar^{2}}\left|\frac{\mathrm{d}\bm{r}_{i}(\tau)}{\mathrm{d}\tau}\right|^{2}+\int\mathrm{d}\bm{r}^{\prime}W_{\mu}(\bm{r}^{\prime})\hat{n}_{\mu}(\bm{r}^{\prime})\right]}\,.$		(3)

The $\hat{n}$ operator implicitly carries the $N$ coordinate dependencies with it, so inserting eqn. 2 into eqn. 3 allows one to see, after some manipulations, that the argument of the exponential in the configuration path integral is now completely separable into $N$ one-body terms, producing a product of $N$ separable path integrals. The configuration integrals can then be evaluated one at a time with the result being $N$ identical terms. The $N$-body partition function $Q_{N}(\beta)$ can then finally be expressed as

$\displaystyle Q_{N}(\beta)=\prod_{\mu}\int\int\mathcal{D}[\mathcal{N}_{\mu}]\mathcal{D}[W_{\mu}]\,e^{-\beta F[\mathcal{N}_{\mu},W_{\mu}]}$

(4)

where

$\displaystyle F[\mathcal{N}_{\mu},W_{\mu}]=-\frac{1}{\beta}N_{\mu}\ln(Q_{\mu}[W](\beta))+U[\mathcal{N}_{\mu}]-\int\mathrm{d}\bm{r}^{\prime}W_{\mu}(\bm{r}^{\prime})\mathcal{N}_{\mu}(\bm{r}^{\prime})$

(5)

and $Q_{\mu}[W](\beta)$ is a single-pair partition function that we have defined according to the expression

$\displaystyle Q_{\mu}[W](\beta)=\int\mathrm{d}\bm{r}\,q_{\mu}(\bm{r},\bm{r},\beta)$

(6)

where $q_{\mu}(\bm{r},\bm{r},\beta)$ represents the propagation of a single pair from initial position $\bm{r}$ at $s=0$ to final position $\bm{r}$ at $s=\beta$. The single-pair propagator $q_{\mu}(\bm{r},\bm{r}^{\prime},s)$ can be expressed as a path integral:

$\displaystyle q_{\mu}(\bm{r},\bm{r}^{\prime},s)=\mathcal{A}\int\mathcal{D}[\bm{r}]e^{-\int^{s}_{0}\mathrm{d}\tau\left[\frac{m}{2\hbar^{2}}\left|\frac{\mathrm{d}\bm{r}(\tau)}{\mathrm{d}\tau}\right|^{2}+W_{\mu}(\bm{r}(\tau))\right]}$

(7)

where $\mathcal{A}$ is a formally infinite normalization constant coming from the kinetic degrees of freedom, whose value we shall not be concerned with since it will not appear in any of the quantities of interest. As a propagator, it can be shown that $q_{\mu}(\bm{r},\bm{r}^{\prime},s)$ equivalently satisfies the modified diffusion equation

$\displaystyle\frac{\partial q_{\mu}(\bm{r},\bm{r}^{\prime},s)}{\partial s}=-H^{\text{eff}}_{\mu}q_{\mu}(\bm{r},\bm{r}^{\prime},s)=\frac{\hbar^{2}}{2m}\nabla^{2}q_{\mu}(\bm{r},\bm{r}^{\prime},s)-W_{\mu}(\bm{r})q_{\mu}(\bm{r},\bm{r}^{\prime},s)$

(8)

with initial condition $q_{\mu}(\bm{r},\bm{r}^{\prime},0)=\delta(\bm{r}-\bm{r}^{\prime})$, which makes it possible to evaluate $Q_{\mu}[W](\beta)$ and thus $F[\mathcal{N},W]$ exactly. It is worth pointing out that the Hamiltonian $H_{\text{eff}}$ above is the same Hamiltonian appearing in the Kohn-Sham equation from KS-DFT, a fact we use to prove the equivalence of the two theories in appendix C. Although it is possible to analytically evaluate the functional eqn. 5, we are not so fortunate with eqn. 4, whose path integrals are too unwieldy to perform exact calculations with to get the free energy [30]. However, eqn. 5 plays the same role that the action does in the real-time quantum mechanical path integral, so a solution can be sought which extremizes $F[\mathcal{N},W]$ by setting its first variation to zero and then approximating the integrand with the extremum of $F$. The free energy can then be calculated from $F[n,w]$, where $n$ and $w$ are the mean-fields for which the functional $F$ has a saddle point [28, 30, 22]. We can further justify the preservation of exactness in the model from varying eqn. 5, since any neglected higher-order contributions can be packaged into the unknown functional $U$, whose approximations occupy a large portion of current DFT research [19, 43, 20]. The total kinetic energy functional $K[n,w]$ can also be calculated from the expectation value of the kinetic term in the many-body Hamiltonian, after some manipulation, to be

$\displaystyle K[n,w]$

$\displaystyle=-\sum^{N}_{j}\bigg{\langle}\frac{\hbar^{2}}{2m}\nabla^{2}_{j}\bigg{\rangle}=-\frac{\hbar^{2}}{2m}\sum_{\mu}\sum^{N}_{j}\frac{1}{Q_{\mu}[w](\beta)}\int\mathrm{d}\bm{r}_{j}\nabla^{2}_{\bm{r}^{\prime}_{j}}q_{\mu}(\bm{r}_{j},\bm{r}^{\prime}_{j},\beta)\bigg{|}_{\bm{r}^{\prime}_{j}=\bm{r}_{j}}\,.$

(9)

Proceeding from the variational principle outlined above, using the path integral form of the single-pair propagator, the mean-field density $n_{\mu}(\bm{r},\beta)$ and field $w_{\mu}(\bm{r},\beta)$ corresponding to each pair are found to be

$\displaystyle w_{\mu}(\bm{r},\beta)=\frac{\delta U[n]}{\delta n_{\mu}(\bm{r},\beta)}\ \ \text{and}\ \ n_{\mu}(\bm{r},\beta)=\frac{N_{\mu}}{Q_{\mu}[w](\beta)}q_{\mu}(\bm{r},\bm{r},\beta)\,.$

(10)

The potential $U[n]$ is the only remaining quantity yet to be specified, which will give us the expressions for the fields $w_{\mu}(\bm{r},\beta)$ experienced by each pair in the system, and finally the expression for the free energy $F[n,w]$. In our model, the electron pairs experience four fields in the vicinity of the atomic nucleus: the Coulomb field between the nucleus and the electron pairs $w_{\mu}^{e-n}(\bm{r},\beta)$, the Coulomb field between the electrons $w_{\mu}^{e-e}(\bm{r},\beta)$, and the exchange field between electron pairs $w_{\mu}^{x}(\bm{r},\beta)$ representing two separate fields. The first of these two fields is the electron self-interaction field $w_{\mu}^{sic}(\bm{r},\beta)$, which corrects for the interaction of the electron with its own field that is not accounted for in the electron-electron Coulomb field $w_{\mu}^{e-e}(\bm{r},\beta)$; the self-interaction correction introduced in prior work [24] is employed in this work as well. The second of these two fields, as is commonplace in all OF-DFT approaches, is the Pauli-exclusion field $w_{\mu}^{P}(\bm{r},\beta)$, which accounts for the repulsion felt by electron pairs with the same configuration attempting to occupy the same location at the same (imaginary) time, as stipulated by the Pauli-exclusion principle. As will be discussed later, the Pauli-exclusion field used in this work accounts for some exchange effects but does not account for correlations.

The electron-nucleus potential $U_{e-n}[n]$ takes the form

$\displaystyle U_{e-n}[n]=-\int\int\mathrm{d}\bm{r}\,\mathrm{d}\bm{r}^{\prime}\,n(\bm{r},\beta)\frac{\rho_{\text{nuc}}(\bm{r}^{\prime})}{\left|\bm{r}-\bm{r}^{\prime}\right|}\,,$

(11)

where $\rho_{\text{nuc}}(\bm{r})$ is the nuclear density, which we take to be $\rho_{\text{nuc}}(\bm{r})=N\delta(\bm{r})$. The electron-nucleus field for each pair $w^{e-n}_{\mu}(\bm{r},\beta)$ is then found to be

$\displaystyle w^{e-n}_{\mu}(\bm{r},\beta)=-\int\mathrm{d}\bm{r}^{\prime}\,\frac{\rho_{\text{nuc}}(\bm{r}^{\prime})}{\left|\bm{r}-\bm{r}^{\prime}\right|}=-\frac{N}{|\bm{r}|}\,.$

(12)

The potential due to electron-electron Coulomb-type interactions $U_{e-e}[n]$ is similarly given by

$\displaystyle U_{e-e}[n]=\frac{1}{2}\int\int\mathrm{d}\bm{r}\,\mathrm{d}\bm{r}^{\prime}\,n(\bm{r},\beta)\frac{n(\bm{r}^{\prime},\beta)}{\left|\bm{r}-\bm{r}^{\prime}\right|}\,,$

(13)

and the electron-electron field for each pair $w^{e-e}_{\mu}(\bm{r},\beta)$ is then found to be

$\displaystyle w^{e-e}_{\mu}(\bm{r},\beta)=\int\mathrm{d}\bm{r}^{\prime}\frac{n(\bm{r}^{\prime},\beta)}{|\bm{r}^{\prime}-\bm{r}|}\,.$

(14)

Both eqn. 12 and eqn. 14 indicate that each pair experiences the exact same electron-electron and electron-nucleus field. The total field experienced by each pair will however not be the same, thanks to the exchange effects introduced by the other two fields.

Following previous work [40], the closest classical analogue of the Pauli-exclusion principle is the notion of excluded-volume, which in polymer SCFT, is often implemented as a Dirac delta energy penalty for overlapping polymer segments. If we are to be truly faithful to the exclusion principle however, the energy penalty should be for overlapping polymer segments from different polymer contours representing pairs of quantum particles to account for spin. Since the position along the polymer contour is parametrized by a parameter $s$, the energy penalty due to overlapping polymer contours occurs only for contours at the same value of $s$. Recall from the quantum-classical correspondence that the parameter $s$ can be interpreted as an imaginary time, so the Pauli-exclusion repulsion is akin to a particle pair feeling the excluded-volume repulsion when at the same place and (imaginary) time as another pair [40]. This idea is difficult to implement in practice however, so we approximate it by projecting out the degrees of freedom from the $s$ parameter space, which effectively amounts to imposing the excluded-volume energy penalties for every value of $s$. The downside to this approximation is that it ignores the inter-contour correlations and will clearly overestimate the excluded-volume felt between the pairs [40]. The approximate Pauli-exclusion potential is then given by

$\displaystyle U_{P}[n]=\frac{1}{2g_{0}}\sum_{\begin{subarray}{c}\mu,\nu\\ \mu\neq\nu\end{subarray}}\int\int\mathrm{d}\bm{r}\mathrm{d}\bm{r}^{\prime}n_{\mu}(\bm{r},\beta)\delta(\bm{r}-\bm{r}^{\prime})n_{\nu}(\bm{r}^{\prime},\beta)=\frac{1}{2g_{0}}\sum_{\begin{subarray}{c}\mu,\nu\\ \mu\neq\nu\end{subarray}}\int\mathrm{d}\bm{r}n_{\mu}(\bm{r},\beta)n_{\nu}(\bm{r},\beta)$

(15)

where $g_{0}$ is a constant with the same units as a density of states [40]. In principle, since the excluded-volume interaction is independent of the system under study, $g_{0}$ should be a universal constant whose value can be determined by comparing the Pauli potential for a very simplistic system (e.g. a uniform gas with only excluded-volume interactions) to experimental results. However, because the Pauli potential is being approximated in this work, $g_{0}$ is taken to be arbitrary and we choose its value once for all calculations; the value chosen and how it was chosen will be discussed in the Results section III. The approximate form of the Pauli-exclusion field for each pair $w^{P}_{\mu}(\bm{r},\beta)$ is then calculated as

$\displaystyle w^{P}_{\mu}(\bm{r},\beta)=\frac{1}{g_{0}}\sum_{\begin{subarray}{c}\gamma\\ \gamma\neq\mu\end{subarray}}n_{\gamma}(\bm{r},\beta)\,.$

(16)

In previous work [24], some constraints on the exact form of the Pauli-exclusion field were given which allowed for the accuracy of the approximate expression eqn. 16 to be assessed. The relevant constraints were as follows:

	$\displaystyle w^{P}(\bm{r},\beta)\geq 0\ \ ,\ \ \raisebox{2.15277pt}{\scalebox{1.0}{$\displaystyle\lim_{|\bm{r}|\rightarrow\infty}\;$}}w^{P}(\bm{r},\beta)=0\ \ ,\ \ \int\mathrm{d}\bm{r}\,w^{P}(\bm{r},\beta)n(\bm{r},\beta)<\infty\ \ ,\ \ w^{P}(\bm{r},\beta)=0\ \ \text{for}\ \ N=2\ \ ,$
	$\displaystyle w^{P}[\lambda^{3}n](\lambda\bm{r},\beta)=\lambda^{2}w^{P}[n](\lambda\bm{r},\beta)$		(17)

where in the last criterion $\lambda$ is a scale factor and the functional dependence of the field on the density has been explicitly reinstated [25]. We found that all but the last of the constraints in eqn. 17 were satisfied by eqn. 16, with $w^{P}(\bm{r},\beta)$ overestimating the excluded-volume interactions by precisely the amount required to fulfill the last constraint in eqn. 17; a point also discussed in a previous work [40].

The electron self-interaction field $w_{\mu}^{\text{sic}}(\bm{r},\beta)$ used in this work was first introduced in reference [24], and is essentially a Fermi-Amaldi self-interaction correction applied to each particle pair; see reference [24] for further discussion. It has the form

$\displaystyle w_{\mu}^{\text{sic}}(\bm{r},\beta)=-\frac{1}{N_{\mu}}\int\mathrm{d}\bm{r}^{\prime}\frac{n_{\mu}(\bm{r}^{\prime},\beta)}{|\bm{r}-\bm{r}^{\prime}|}$

(18)

where the corresponding potential $U_{\text{sic}}[n]$ is simply

$\displaystyle U_{\text{sic}}[n]=-\sum_{\mu}\frac{1}{2N_{\mu}}\int\int\mathrm{d}\bm{r}\mathrm{d}\bm{r}^{\prime}\frac{n_{\mu}(\bm{r}^{\prime},\beta)n_{\mu}(\bm{r},\beta)}{|\bm{r}-\bm{r}^{\prime}|}\,.$

(19)

In this form, eqn. 18 directly preserves the desirable qualities of the original Fermi-Amaldi electron self-interaction correction for hydrogen and helium [3]. Furthermore, because eqn. 18 acts on electron pairs and $N_{\mu}=\{1,2\}$, then eqn. 18 effectively accounts for the self-interaction of every electron in the atomic system. Those familiar with the Hartree-Fock model will immediately recognize that our model with eqn. 18 is identical to the situation in the Hartree-Fock model, which models exchange effects exactly. However, because we are approximating the Pauli-exclusion field by projecting out the degrees of freedom from the $s$ parameter space, our model in its current implementation will not reproduce the precise binding energies predicted by Hartree-Fock theory. This is because our Pauli-exclusion field overestimates the excluded-volume felt by the polymer contours in the $s$-parameter space, hence electron pairs feel too much repulsion between each other and the electron shells will be too distant from their neighbours, raising the free energy substantially in some cases. The Hartree-Fock model on the other hand, is a wavefunction-based model, so the Pauli-exclusion effect is automatically encoded into the wavefunction due to the spin-statistics theorem [36]. Implementing the exact expression for the Pauli potential would then allow our model to coincide exactly with Hartree-Fock theory.

As was mentioned previously, electron pairs are taken to be ring polymers embedded in a 3+1-dimensional thermal space under the influence of a potential $U$. The ring polymers are confined to explore the thermal space according to this potential, which they can do using two different mechanisms: translation and configuration. The translational degrees of freedom refer to the three dimensional motion of the polymer as a whole, while the configurational degrees of freedom refer to how the polymer confirmations change while holding one point of the polymer fixed in space. Each of these mechanisms has an entropy associated to them, which we will denote as $S_{t}$ and $S_{c}$, respectively. The behaviour of the polymer can then be explained by looking at the competition between the degrees of freedom encoded in the entropies and those restricted by the potential $U$. Therefore, by rephrasing the free energy per pair $F[n_{\mu},w_{\mu}]$ in terms of these entropies, we can exactly describe the process in the polymer picture by which the electrons surrounding the atomic nucleus would break spherical-symmetry; and with the notion of pairs, we can pinpoint exactly which pairs affect this change. Following previous work 39 and reference 29, the free energy per pair $F[n_{\mu},w_{\mu}]$ can be re-expressed as

	$\displaystyle F[n_{\mu},w_{\mu}]$	$\displaystyle=-\frac{1}{\beta}\int\mathrm{d}\bm{r}^{\prime}n_{\mu}(\bm{r}^{\prime},\beta)\ln{\left(\frac{n_{\mu}(\bm{r}^{\prime},\beta)}{N_{\mu}}\right)}+U[n_{\mu}]$
		$\displaystyle\quad+\frac{1}{\beta}\int\mathrm{d}\bm{r}^{\prime}n_{\mu}(\bm{r}^{\prime},\beta)\left[\ln{\left(q_{\mu}(\bm{r}^{\prime},\bm{r}^{\prime},\beta)\right)}+\beta w_{\mu}(\bm{r}^{\prime},\beta)\right]$		(20)

where we can identify the last term with the free energy contribution coming from the configurational entropy $S_{c}[n_{\mu},w_{\mu}]$ ¹¹1This configurational entropy expression corrects a typo in reference 39 which is missing a factor of $-\beta$. and the first term with the translational entropy $S_{t}[n_{\mu}]$, which we write as

	$\displaystyle S_{c}[n_{\mu},w_{\mu}]$	$\displaystyle=-\int\mathrm{d}\bm{r}^{\prime}n_{\mu}(\bm{r}^{\prime},\beta)\left[\ln{\left(q_{\mu}(\bm{r}^{\prime},\bm{r}^{\prime},\beta)\right)}+\beta w_{\mu}(\bm{r}^{\prime},\beta)\right]$		(21)
	$\displaystyle S_{t}[n_{\mu}]$	$\displaystyle=\int\mathrm{d}\bm{r}^{\prime}n_{\mu}(\bm{r}^{\prime},\beta)\ln{\left(\frac{n_{\mu}(\bm{r}^{\prime},\beta)}{N_{\mu}}\right)}\,.$		(22)

We shall examine the intuitive polymeric interpretation of spherical symmetry-breaking in the Results III and Discussion IV sections.

To solve the set of self-consistent equations 8,6, and 10, the biggest obstacle is the solution to the modified diffusion equation eqn. 8, which needs to be solved a number of times corresponding to the total number of pairs $N_{p}$ in the system — per self-consistent iteration — for every value of both spatial positions $\bm{r}$ and $\bm{r}^{\prime}$. This double spatial dependence of $q_{\mu}(\bm{r},\bm{r}^{\prime},s)$ means that traditional real-space methods are impractical for computational efficiency [39, 30, 28]. Instead, what is usually done in the polymer SCFT community, is to use the spectral method: all spatially-dependent functions are decomposed in terms of a set of basis functions, for which the position dependence of each function can be integrated out and the resulting equations become matrix equations for the unknown expansion coefficients. The problem encountered earlier with real-space methods is then made to vanish and is replaced with solving a matrix equation $N_{p}$ times per self-consistent iteration [30, 39]. In this work we use the spectral method with non-orthogonal Gaussian basis sets outlined in reference 24 to solve the modified diffusion equation. The method, along with the numerical procedure, has been discussed in reference 24, however, only the case of spherical Gaussian basis functions was treated. To solve the angular problem, we will need the full angular Gaussian basis functions, which are given by the (normalized) expression

$\displaystyle f_{i}(\bm{r})=\left(2(2c_{pl})^{l+\frac{3}{2}}\left[\Gamma\left(l+\frac{3}{2}\right)\right]^{-1}\right)^{\frac{1}{2}}Z^{m}_{l}(\theta,\phi)r^{l}e^{-c_{pl}r^{2}}\,,$

(23)

where the index $i$ represents the tuple of indices $(p,l,m)$ and $Z^{l}_{m}(\theta,\phi)$ are the real spherical harmonics defined in appendix B. From here, the basis-specific quantities outlined in appendix B of the overlap matrix (eqn. 50), Laplace matrix (eqn. 52), and Gamma tensor (eqn. 53) must be recalculated, but the general non-orthogonal spectral equations listed in appendix A remain unchanged. Although the differences from the spherical case only appear in the basis-specific quantities, it should be pointed out that the addition of angularity increases the dimensions of the basis-specific quantities substantially, severely increasing the computational runtime for even a modest sized basis set, and also introduces a complicated quantity into the expression for the Gamma tensor known as the real Gaunt coefficients; all details can be found in appendix B and references 24, 39.

The Gaussian exponents $c_{pl}$ are chosen according to an “even-tempered”-type scheme outlined in reference 24. In this work $N_{b}=425$ Gaussian basis functions were used for the angular results since this number provided excellent resolution and converged far enough to the infinite basis set limit, although more basis functions are used in this case compared with reference 24 because we assigned numbers for each $l$ value of the real spherical harmonics, which in turn have $2l+1$ types of basis functions (for the number of $m$ values associated to each $l$). Convergence here is judged according to the spectral convergence criterion used in reference 24. The number 425 comes from assigning 150 basis functions to $l=0$, 50$\times$3 basis functions to $l=1$, and 25$\times$5 basis functions to $l=2$. Increasing $l$ values are assigned smaller numbers of basis functions because the corresponding basis functions become much more diffuse and start to represent smaller portions of the electron density profile. The angular results are also only presented for the first 10 elements, so 150 basis functions for the $l=0$ portion is more than enough to achieve good resolution; the approximations used in this work also limit our accuracy a lot more than the basis set truncation error does. The set of minimum exponents we chose was $c_{1,1,1}=(10^{-15},10^{-10},10^{-6})$ and the set of maximum exponents was $c_{150,50,25}=(10^{11},10^{5},10^{3})$, where each entry corresponds to $l$ values in increasing order, respectively. Anything higher than 425 would only add a small amount of resolution and would require a relatively large increase in computation time. For 425 basis functions, every element was able to satisfy a tolerance of at least $10^{-6}$, with some going as far as $10^{-8}$. We decided to keep the same value for $g_{0}$ that was used in reference 24 of $g_{0}=0.1$ for the arbitrary constant $g_{0}$ associated to the Pauli-exclusion field eqn. 16, since this will make comparison with previous results easier. The full angular Gaussian basis set used in this work does not introduce or change any previously encountered numerical considerations addressed in reference 24.

The atomic binding energies corresponding to $g_{0}=0.1$ for the elements hydrogen to neon, strictly enforcing the maximum occupancy of 2 electrons per pair, are shown in table LABEL:tab2 for the spectral expansion of the density with angular basis functions, and table LABEL:tab3 for the restriction to spherically-symmetric basis functions; atomic units are used unless otherwise specified. The binding energies predicted by our model are contrasted with those predicted by Hartree-Fock theory, since in our neglect of correlations and use of an exact self-interaction correction, Hartree-Fock theory should be considered as “exact”. One difference with the results from reference 24 can be seen by looking at the percentage deviation with Hartree-Fock for the two tables: The model with angular dependence is much closer to Hartree-Fock for the first 6 elements, but rapidly worsens due to an overabundance of Pauli-exclusion repulsion felt between electron pairs. In the angular case considered here, the atomic shell configurations are not prebuilt into the code, so the overabundance of the Pauli-exclusion force causes the pairs after carbon to spread too far apart. However, the observed shell structure does arise solely from the information provided by the electron pair configurations and still somewhat resembles what we expect (figures 2-7)²²2The density contour plots of the first four elements are not shown as they are all spherical and have identical pair to shell structure. Reference 24 can be consulted for plots of hydrogen to beryllium using the model from this work.. Moreover, we do see spontaneous spherical symmetry-breaking, which is first predicted to occur at carbon as opposed to boron. The shapes of the pair densities after boron do not match naive expectations; this is addressed later in this section. The magnitude of the symmetry-breaking effect can also be modified through the value of $g_{0}$ (i.e. a smaller value produces more noticeable deviations). In particular, if $g_{0}$ is given pair dependence, then a simple arithmetic sequence of increasing values chosen from the neighbourhood around $0.1$ such that the smallest value is assigned to the first pair, is sufficient for boron to break spherical symmetry. Although nature predicts spherical symmetry-breaking to first occur at boron, this does not mean the polymer excluded-volume picture is invalid: symmetry-breaking has a subtle effect on the binding energies, as can be seen by comparing tables LABEL:tab3 and LABEL:tab2, so the approximation used for the Pauli-exclusion field may be too coarse to allow such a prediction, in which case the exact field would need to be implemented to sufficiently test this. The fact that spontaneous spherical symmetry-breaking does occur, and only one element off from where it is supposed to be, is a very encouraging result.

There is no difference in the binding energies between the spherical and angular results for the first 4 elements since these elements are known to have spherical ground-state distributions and minimal Pauli-exclusion repulsion. The lack of a binding energy difference between the spherical and non-angular cases of boron is attributed to the approximate Pauli potential used in this work, which predicts symmetry-breaking to occur at carbon instead of boron. Carbon is the first element where any difference in the binding energy between the spherical and non-angular cases can be seen: the percent difference of the SCFT model with the prediction from Hartree-Fock theory in the angular case is approximately an order of magnitude smaller than in the spherical case. The agreement for carbon is due in part to the cancellation of certain errors as opposed to a genuine agreement with Hartree-Fock theory, but the other trends in tables LABEL:tab3 and LABEL:tab2 suggest that the angular case does yield an improvement in the agreement with Hartree-Fock theory. The rest of the elements from tables LABEL:tab3 and LABEL:tab2 display very minor changes in the binding energy, which agrees nicely with the results of Chowdhury and Perdew [8], who report that the effect of symmetry-breaking has a small impact on the binding energy. These results suggest that the spherical-averaging approximation used in reference 24 performs quite well in most scenarios, and that it is physically reasonable to use it for isolated atoms provided the aim is not to investigate delicate features that arise from angularity.

In figures 2-7, we see that the first pair density always remains spherically-symmetric, which makes sense as this pair corresponds to the innermost electrons in the atom whose density profile is dominated by the spherical 1s contribution. In fact, because the magnitude of the density for the first pair is so much larger than the rest, the total density profile only marginally deviates from a spherical distribution. The pair densities beyond the first pair do not correspond to the orbital picture we get from other DFT approaches, and the non-spherical pair densities only resemble a single lobe in contrast to the multiple lobes expected from orbital pictures such as Hartree-Fock theory. However, pair densities do not correspond to the squared modulus of individual orbitals from KS-DFT, rather, they correspond to sums of squared moduli of orbitals, as shown in appendix C. The pair density profiles corresponding to non-spherical pairs somewhat resemble the situation in polymer macro-phase separation [28], where one pair occupies one of the lobes in one region and the other occupies the partner lobe across from it, together forming a hybrid 2s-2p-like structure. This macro-phase behaviour is unsurprising given that the system is being modelled as a classical polymeric system with higher-dimensional excluded-volume, a model which is equivalent to the wavefunction picture through the theorems of DFT [41]. The total densities of the atoms seen in figures 3-7 have approximately the same profiles as the densities predicted by quantum mechanics in the wavefunction picture, but due to the inexact Pauli-exclusion field used in this work, we do not expect to produce identical profiles.

In light of the macro-phase-like behaviour seen in figures 3-7, there are a number of constraints that the true electron density must satisfy in order to be considered physically acceptable, which we can use as benchmarks to assess the density profiles predicted by the model used in this work. The two main constraints on the electron density are that it must be positive for all positions $\bm{r}$ and that its integral over all space must yield the electron number $N$ (or the pair electron number if we are dealing with an individual pair density). Two further constraints are

$\displaystyle 1\geq\frac{3\pi}{4K}\left[\frac{\pi}{2}\int\mathrm{d}\bm{r}\,n^{3}(\bm{r},\beta)\right]^{\frac{1}{3}}\ \ \text{and}\ \ 1\geq\frac{1}{2K}\int\mathrm{d}\bm{r}\,\left(\nabla\sqrt{n(\bm{r},\beta)}\right)^{2}$

(24)

where $K$ represents the kinetic energy of the system [26], which is given by eqn. 43 in the model. The first constraint of eqns. 24 is the requirement that the electron density be contained in the function space $L^{3}$, while the second is the requirement that the kinetic energy associated with the electron density be bounded below by the von Weizsacker kinetic energy [26].

The pair density profiles corresponding to figures 2-7 clearly demonstrate that the pair densities, and thus the total density, are non-negative for all positions $\bm{r}$, since we know the density goes to zero for large $r$. Likewise, the expected electron numbers corresponding to each pair were obtained from the pair density integrals to within numerical accuracy (e.g. basis set truncation), adding up to the desired total electron number in every case. Table LABEL:tab5 displays the right-hand side values of eqns. 24 for the elements hydrogen to neon, where we can see that the density predicted by the model always satisfies these two inequalities. In particular, one can notice that the right-hand side values of the second inequality in eqns. 24 for hydrogen and helium are almost exactly 1, which is the statement that the von Weizsacker kinetic energy functional is exact for one and two electron systems.

Tables LABEL:tab4 and LABEL:tab6 list the pairwise components of the potential energy terms, the energy contributions from the configurational and translational entropies, and the total free energies for the elements carbon and fluorine, respectively. In both elements we can see that the lowering of the free energy due to spherical symmetry-breaking is produced from pairs 2 and 3, which both adopt opposing lobe shapes, in contrast to the purely spherical distributions in the spherically-averaged case, that distributes their free energy contribution uniformly amongst the two. It is this feature in particular that accounts for the free energy difference between the spherical and angular cases. Looking more closely at tables LABEL:tab4 and LABEL:tab6, the cause of this feature is the fact that the third pair density can occupy a region closer to the atomic nucleus by violating spherical symmetry, which is evidenced by the much lower $U_{e-n}$ value in both elements for the angular case. The electron-electron plus self-interaction correction potential, the Pauli potential, and both entropic contributions to the free energy for pair 3 are all worse in the angular case, suggesting that the move towards the nucleus more than compensates for the interaction penalties with other pairs. In order for pair 3 to accomplish the move from a spherical distribution to a lobe distribution, pair 2 must also transform to a mirroring lobe distribution so that its original spherical shape does not overlap as much with the new lobe shape of pair 3; this effect is propagated with the other pairs (except the first), converting them into non-spherical distributions as well.

Despite the failure of the Pauli-exclusion field to satisfy the coordinate scaling relation eqn. 17, the scaling argument presented in reference 24 shows that the picture of higher-dimensional excluded-volume interactions between pairs of threads recovers the Thomas-Fermi quantum kinetic energy term and the Dirac exchange term in the uniform density limit. The analysis from section III also demonstrates that the pair densities making up the total density satisfy all constraints necessary to guarantee a physically acceptable electron density. Together with the proof that pair densities do not necessarily correspond to individual squared moduli of orbitals, and the formal equivalence of the polymer-thread picture with quantum DFT through the quantum-classical isomorphism [12], the macro-phase-like behaviour exhibited in figures 3-7 show that spontaneous shell structure and spherical-symmetry breaking are robust predictions. The symmetry-breaking arises from the energetic benefit of electrons distributing asymmetrically closer to the nucleus. It is also clear from adding up the individual pair density profiles for any given atom, that the total density profiles only deviate slightly from spherical symmetry, which is consistent with the findings of Chowdhury and Perdew [8] that asymmetries in the electron density have a small effect. Together, these two observations highlight an important distinction: the ability of the pairs to individually break spherical-symmetry allows the atom to lower its binding energy in all cases, but this does not necessarily mean that the total electron density also breaks spherical symmetry.

One should consider whether the macro-phase-like behaviour is simply an artefact of the specific approximation for the Pauli potential being used in this work, as it is reasonable to speculate that the exaggerated repulsion produced by the approximate Pauli-exclusion field causes the pairs to clear their local neighbourhood. That is, are we only observing isolated atoms to be spherically-asymmetric because we are using an approximate model? If we implemented the exact Pauli potential, would total electron densities always be found to be spherically symmetric, consistent with the arguments of references 9, 13?. This seems unlikely, since this would require an unjustifiable perfect balance between the Pauli potential and other factors. Other frustrated systems induce spontaneous symmetry breaking, for example in true polymeric systems, SCFT is used to predict the micro-phases of block copolymers [30, 29, 28]. Returning to the present model, the shell structure for carbon displayed in figure 3 demonstrates non-spherical structure yet maintains nearly identical shell structure to that predicted by Hartree-Fock theory. As was mentioned earlier, there is probably some cancellation of errors happening within carbon due in part to competing Pauli pair repulsions, but it seems unlikely based on the trends from the other atoms combined with the difference in sensitivity between the binding energies and the density profiles, that the magnitude of this effect would be large enough to account for the macro-phase-like structure while also producing the minute differences with Hartree-Fock theory that are observed.

The ring polymer SCFT formulation of quantum mechanics predicts the spontaneous emergence of atomic shell structure and spherical symmetry-breaking in isolated neutral atoms hydrogen to neon Using postulated pair structure of the model and ideas of higher-dimensional excluded-volume in cooperation with an exact self-interaction correction, the model shows excellent agreement with Hartree-Fock theory for the atomic binding energies and density profiles of the first six elements, providing exact matches for the elements hydrogen and helium. However, due to the approximation made on the Pauli-exclusion field, the predicted shell structure starts to deviate significantly past the element neon and the symmetry-breaking is first predicted to occur at carbon instead of boron. Consistent with Chowdhury and Perdew [8], the symmetry-breaking effect is found to have a very small impact on the binding energies, which suggests that the spherical-averaging approximation is physically reasonable when investigating atomic systems. The pair density contour plots also display behaviour similar to polymer macro-phase separation, where individual electron pairs occupy single lobe structures that together form a dumbbell shape analogous to the 2p orbital shape. It is further shown that the predicted densities satisfy known constraints and still produce the same total electronic density profile that is predicted by quantum mechanics.

There are a number of future directions to consider, now that the basic engine from reference 24 has been constructed. One possible direction could be to extend the work of reference 35 in modelling systems of diatomic molecules to arbitrary formations of molecules or even solid-state lattices, since the Gaussian methodology is easily adapted to any number of complex geometries. The initial groundwork involved in this direction would involve switching to contracted Gaussian basis sets [16, 18, 15], since they allow for many fewer basis functions to be used while still maintaining roughly the same level of precision and resolution; the price tag associated to the contracted sets comes in the form of additional minimization routines that either minimize the spectral representation of the free energy eqn. 42 with respect to even-tempered parameters [33, 16] or fit the Gaussians to a Slater-type function [15, 16]. Both of the these methods typically require derivative information to perform the minimization, which is undesirable because the derivatives may not be well-defined or could possibly lead to numerical instabilities yielding false minima. A method that uses only the Nelder-Mead algorithm was originally developed for this work in anticipation of investigating molecular systems, which uses the spectral coefficients of the density and the Gaussian exponents from the uncontracted result to solve for the contraction coefficients, and then minimizes the sum of squared deviations between the two to solve for the exponents of the contracted set. After implementing the contracted sets and updating the current computational engine, one would need to generalize the centres of the Gaussian basis sets to arbitrary positions and modify the structure of the computation to accommodate multiple atoms. The shifting of the centres of the Gaussians from the origin to arbitrary positions is not so easily done with the spherical harmonic representation used in this work and might be better facilitated using a Cartesian representation of the Gaussians [16], which only entails re-deriving the basis specific matrices and the spectral components of the electron-nucleus field. Fortunately, the expressions for these quantities have already been derived, although the expressions are much more complicated. In the case of solid-state lattices, pseudo-potentials and other modifications would need to be introduced as well [20, 43]. Further increases in accuracy could also be achieved by combining the molecular dynamics framework of Car and Parrinello 5 with the model, to better account for the nuclear degrees of freedom.

Another possible direction is to implement an exact Pauli-exclusion field, so that the comparison of the present model with Hartree-Fock theory can be completed, and the equivalence of higher-dimensional excluded-volume with the Pauli-exclusion principle can be verified. This direction would be an important test of a foundational aspect of the theory [41], and would provide further evidence in support of the symmetry breaking mechanism presented here, in which electrons lower the free energy by breaking spherical symmetry to form electron distributions that approach closer to the atomic nucleus.

Other possible directions include adding correlation terms and relativistic corrections, or clarifying the mechanism for electron spin in the model. Correlation fields are stipulated to be the only thing missing from this model that prevents it from completely agreeing with the predictions of non-relativistic quantum mechanics for neutral atoms not in the presence of any other fundamental fields. Therefore, finding a mechanism for correlations within the polymer picture would also be of foundational importance to the model. An investigation of quantum correlations would also naturally lead to the topic of quantum entanglement, which further connects with information-theoretic approaches to DFT [31, 1] that may be useful in learning more about properties of correlations in many-body systems from the perspective of DFT. Investigating the mechanism that represents electron spin would also complement the study of correlations. Lastly, relativistic corrections including fine structure [27, 2] and finite-size nuclear centres [42] could be added to the model in order to study heavier elements, which may yield useful information on how relativistic effects manifest themselves in the polymer picture.

This research was financially supported by the Natural Sciences and Engineering Research Council of Canada (NSERC).