Fermionic Phonons: Exact Analytic Results and Quantum Statistical Mechanics for a One Dimensional Harmonic Crystal

TL;DR

This paper derives exact analytic solutions for the energy states of a one-dimensional harmonic crystal and explores their thermal properties, revealing the need for many energy levels for accurate modeling at moderate temperatures.

Contribution

It provides the first exact analytic expressions for energy eigenvalues and eigenfunctions of a 1D harmonic crystal, along with numerical analysis of thermal properties.

Findings

Large number of energy levels needed for accurate results

Analytic expressions for eigenvalues and eigenfunctions

Temperature-dependent energy and density profiles

Abstract

Analytic expressions for the energy eigenvalues and eigenfunctions of a one-dimensional harmonic crystal are obtained. The average energy and density profiles are obtained numerically as a function of temperature. A surprisingly large number of energy levels (eg.\ 5,000 for 4 particles) are required for reliable results at even moderate temperatures.

Equations

$$U(\mathbf{q})$$

$$K_{jk}=\left\{\begin{array}[]{ll}K,&j=k\\ 1,&j=k\pm 1,\\ 0,&\mbox{ otherwise,}\end{array}\right.$$

$$S_N(K-\mu )$$

$$=$$

$${\mu }_n=K+2\cos \frac{n\pi }{N+1},n=1,2,\dots ,N.$$

$$u_{n,j}=\sqrt{\frac{2}{N+1}}\sin \frac{jn\pi }{N+1},j=1,2,\dots ,N.$$

$$\underline{\underline{X}}\equiv\{{\bf u}_{1},{\bf u}_{2},\ldots,{\bf u}_{N}\}=\left(\begin{array}[]{cccc}u_{1,1}&u_{2,1}&\ldots&u_{N,1}\\ u_{1,2}&u_{2,2}&\ldots&u_{N,2}\\ \vdots&&\ddots&\vdots\\ u_{1,N}&u_{2,N}&\ldots&u_{N,N}\end{array}\right),$$

$$U(\mathbf{q})$$

$${\mathbf{q}}^{\prime }\equiv {\underline{\underline{X}}}^{\mathrm{T}}\mathbf{d}.$$

$$\mathcal{K}(\mathbf{p})=\frac{1}{2m}\mathbf{p}\cdot \mathbf{p}=\frac{1}{2m}{\mathbf{p}}^{\prime }\cdot {\mathbf{p}}^{\prime },$$

$$\mathcal{H}(\mathbf{\Gamma })=\mathcal{K}(\mathbf{p})+U(\mathbf{q})=\frac{1}{2m}{\mathbf{p}}^{\prime }\cdot {\mathbf{p}}^{\prime }-\frac{\lambda }{2}\underline{\underline{D}}:{\mathbf{q}}^{\prime }{\mathbf{q}}^{\prime }.$$

$$\hat{\mathcal{H}}$$

$$E_{n,l_n}=(l_n+0.5)\hslash {\omega }_n,l_n=0,1,2,\dots ,$$

$${\varphi }_{n,l_n}(Q_n)\equiv \frac{1}{\sqrt{2^{l_n}{l}_n!\sqrt{\pi }}}{e}^{-Q_n^2/2}{\mathrm{H}}_{l_n}(Q_n),$$

$${\varphi }_{\mathbf{l}}(\mathbf{Q})=\prod _{n=1}^N{\varphi }_{n,l_n}(Q_n),$$

$${\psi }^{\pm }(\mathbf{x})\equiv \frac{1}{\sqrt{{\chi }^{\pm }N!}}\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p\psi (\hat{\mathrm{P}}\mathbf{x}),$$

$${\chi }^{\pm }\equiv \sum _{\hat{\mathrm{P}}}(\pm 1{)}^p⟨\psi (\hat{\mathrm{P}}\mathbf{q})|\psi (\mathbf{q})⟩.$$

$$\psi (\mathbf{x})=\sum _{\mathbf{s},\mathbf{n}}⟨{\alpha }_{\mathbf{s}}{\varphi }_{\mathbf{n}}|\psi ⟩{\Phi }_{\mathbf{n},\mathbf{s}}(\mathbf{x}),{\Phi }_{\mathbf{n},\mathbf{s}}(\mathbf{x})\equiv {\alpha }_{\mathbf{s}}(\mathbf{\sigma }){\varphi }_{\mathbf{n}}(\mathbf{q}).$$

$${\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\mathbf{x})=\frac{1}{\sqrt{{\chi }_{\mathbf{n},\mathbf{s}}^{\pm }N!}}\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p{\alpha }_{\mathbf{s}}(\hat{\mathrm{P}}\mathbf{\sigma }){\varphi }_{\mathbf{n}}(\hat{\mathrm{P}}\mathbf{q}).$$

$$\Phi_{{\bf n},{\bf s}}^{\pm}({\bf x})=\left\{\begin{array}[]{l}\displaystyle\frac{1}{\sqrt{\tilde{\chi}^{+}_{{\bf n},{\bf s}}}}\left[\tilde{\alpha}^{+}_{\bf s}({\bm{\sigma}})\tilde{\phi}^{+}_{\bf n}({\bf q})+\tilde{\alpha}^{-}_{\bf s}({\bm{\sigma}})\tilde{\phi}^{-}_{\bf n}({\bf q})\right]\\ \displaystyle\rule{0.0pt}{17.07182pt}\frac{1}{\sqrt{\tilde{\chi}^{-}_{{\bf n},{\bf s}}}}\left[\tilde{\alpha}^{+}_{\bf s}({\bm{\sigma}})\tilde{\phi}^{-}_{\bf n}({\bf q})+\tilde{\alpha}^{-}_{\bf s}({\bm{\sigma}})\tilde{\phi}^{+}_{\bf n}({\bf q})\right],\end{array}\right.$$

$${\tilde{\varphi }}_{\mathbf{n}}^{\pm }(\mathbf{q})\equiv \frac{1}{\sqrt{N!}}\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p{\varphi }_{\mathbf{n}}(\hat{\mathrm{P}}\mathbf{q}),$$

$${\tilde{\alpha }}_{\mathbf{s}}^{\pm }(\mathbf{\sigma })\equiv \frac{1}{\sqrt{N!}}\sum _{\hat{\mathrm{P}}}(\pm 1{)}^p{\alpha }_{\mathbf{s}}(\hat{\mathrm{P}}\mathbf{\sigma }).$$

$${\chi }_{\mathbf{n}}^{\pm }\equiv ⟨{\tilde{\varphi }}_{\mathbf{n}}^{\pm }|{\tilde{\varphi }}_{\mathbf{n}}^{\pm }⟩,\text{mbox}and{\chi }_{\mathbf{s}}^{\pm }\equiv ⟨{\tilde{\alpha }}_{\mathbf{s}}^{\pm }|{\tilde{\alpha }}_{\mathbf{s}}^{\pm }⟩.$$

$${\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\hat{\mathrm{P}}\mathbf{x})=(\pm 1{)}^p{\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\mathbf{x}),$$

$${\tilde{\chi }}_{\mathbf{n},\mathbf{s}}^{\pm }$$

$$=$$

$${\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\mathbf{x})=\frac{1}{\sqrt{N!{\chi }_{\mathbf{n},\mathbf{s}}^{\pm }}}\sum _{P=1}^{N!}(\pm 1{)}^P{\alpha }_{\mathbf{s}}({\mathbf{\sigma }}_P){\varphi }_{\mathbf{n}}({\mathbf{q}}_P).$$

$${\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\mathbf{x})$$

$${\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\mathbf{x})=\frac{1}{\sqrt{2{\chi }_{\mathbf{n},\mathbf{s}}^{\pm }}}\left[{\alpha }_{\mathbf{s}}(\mathbf{\sigma }){\varphi }_{\mathbf{n}}(\mathbf{q})\pm {\alpha }_{\mathbf{s}}({\mathbf{\sigma }}^{\prime }){\varphi }_{\mathbf{n}}({\mathbf{q}}^{\prime })\right].$$

$${\Phi }_{\mathbf{n},\mathbf{s}}^{\pm }(\mathbf{x})$$

$${\tilde{\alpha }}_{\mathbf{s}}^{\pm }(\mathbf{\sigma })$$

Full text

I Introduction

Exact analytic results for model systems have proved useful in the development of the physical sciences. They provide benchmarks against which to test approximate techniques, and they give insight into the mechanisms for the physical behavior of more realistic systems. Although generally those models that are amenable to exact analytic solution are necessarily a simplification of reality, the results do have the great advantage of being unambiguous and free from doubts concerning approximations, convergence, numerical techniques etc. Further, the parameter space may be rapidly explored, allowing general conclusions to be drawn, and numerically sensitive regimes to be identified.

This paper treats a model of a one-dimensional crystal with nearest neighbor harmonic interactions. The potential energy is a quadratic form, with energy eigenvalues and eigenvectors being obtained explicitly. These allow the vibrational modes to be identified, and the exact quantum mechanical solution to be invoked.

The model is realistic in that it includes particle interactions. It is therefore more sophisticated than the quantum ideal gas, or the independent quantum harmonic oscillator that are the routine examples studied by quantum statistical mechanics. Because of the interactions, new insight is provided into phenomena such as collective motion, vibrational excitations, and phonons that is not available with the ideal systems.

An example of the utility of the present exact model calculations is that even for a small system (eg. 4 or 5 particles), 5,000 energy levels are required to obtain quantitatively accurate results at the moderately high temperatures where classical effects become noticeable, which is typical for terrestrial condensed matter. The computational advantage of the present analytic model can be quantified by comparison with the work of Hernando and Vaníček,Hernando13 who, for 4 or 5 interacting Lennard-Jones particles, obtained numerically 50 energy eigenvalues.

II Exact Analysis

II.1 Model

Following earlier work by the author,TDSM consider a one-dimensional harmonic crystal in which the particles are attached by linear springs to each other and to lattice sites. Let the coordinate of the $j$th particle be $q_{j}$, and let its lattice position (ie. in the lowest energy state) be $\overline{q}_{j}=j\Delta_{q}$. The lattice spacing is also the relaxed inter-particle spring length. There are fixed ‘wall’ particles at $q_{0}=0$ and $q_{N+1}=(N+1)\Delta_{q}$. Let $d_{j}\equiv q_{j}-\overline{q}_{j}$ be the displacement from the lattice position; for the wall particles, $d_{0}=d_{N+1}=0$. The system has over-all number density $\rho=\Delta_{q}^{-1}$.

In this model, there is an external harmonic potential of spring constant $\kappa$ acting on each particle centered at its lattice site. The inter-particle spring has strength $\lambda$ and relaxed length $\Delta_{q}$. With these the potential energy is

[TABLE]

Here $\underline{\underline{K}}^{(N)}$ is an $N\times N$ tridiagonal matrix with elements

[TABLE]

where $K\equiv-2-\kappa/\lambda$.

It should be mentioned explicitly that the lattice positions are in order, $\overline{q}_{j}<\overline{q}_{j+1}$, $j=1,2,\ldots,N$. However for the particle positions themselves there is no similar constraint on their order.

It is axiomatic in quantum mechanics that the Hamiltonian operator must be fully symmetric with respect to the permutation of identical particles. Messiah61 ; Merzbacher70 Otherwise the symmetry of the wave function would not be preserved during its evolution, which is to say that bosons would decay into fermions, and vice versa. The above potential energy is not symmetric (eg. in transposing particle positions $q_{j}$ and $q_{k}$, one replaces the respective one-body terms with $\kappa[q_{k}-\overline{q}_{j}]^{2}/2$ and $\kappa[q_{j}-\overline{q}_{k}]^{2}/2$, as well as the respective pair interactions with $\lambda[q_{j\pm 1}-q_{k}\mp\Delta_{q}]^{2}/2$ and $\lambda[q_{k\pm 1}-q_{j}\mp\Delta_{q}]^{2}/2$, which changes the value of the potential energy). This means that the particles are not identical and that the wave function is not restricted by symmetrization. The original version of this paper was vitiated by the author’s misunderstanding of this point.

II.2 Eigenvalues and Eigenvectors

The eigenvalues of the potential energy matrix are required, and these may be obtained from the characteristic equation,

[TABLE]

This is just the recursion relation for the Tchebyshev polynomials of the second kind, which are denoted $S_{n}(x)=U_{n}(x/2)$ by Abramowitz and Stegun, Eq. (AS22.7.6).Abramowitz70 They give $S_{N}(2\cos\theta)=U_{N}(\cos\theta)=\sin((N+1)\theta)/\sin\theta$, which evidently vanishes when $\theta_{n}=\pm n\pi/(N+1)$, $n=1,2,\ldots,N$. Hence the characteristic equation vanishes when $K-\mu_{n}=2\cos\theta_{n}$, or

[TABLE]

These are the eigenvalues, and since $K\equiv-2-\kappa/\lambda$, they are negative.

The corresponding eigenvectors, ${\bf u}_{n}$, have elements

[TABLE]

It may be shown that these form an orthonormal set.

II.3 Normal Modes

The matrix of eigenvectors,

[TABLE]

is orthogonal, $\underline{\underline{X}}^{\mathrm{T}}\underline{\underline{X}}=\underline{\underline{X}}\,\underline{\underline{X}}^{\mathrm{T}}=\mathrm{I}$. (It is also symmetric, $\underline{\underline{X}}^{\mathrm{T}}=\underline{\underline{X}}$.) With this the potential energy may be written

[TABLE]

Here $\underline{\underline{D}}\equiv\underline{\underline{X}}^{\mathrm{T}}\underline{\underline{K}}\,\underline{\underline{X}}$ is diagonal, $D_{n,n^{\prime}}=\mu_{n}\delta_{n,n^{\prime}}$, and the modes are defined as

[TABLE]

Suppose that the particles have mass $m$, and that the momentum of the $j$th particle is $p_{j}=m\dot{q}_{j}$. Hence the classical kinetic energy is

[TABLE]

where ${\bf p}^{\prime}\equiv\underline{\underline{X}}^{\mathrm{T}}{\bf p}$. Accordingly the classical Hamiltonian is

[TABLE]

Evidently, the normal modes represent independent harmonic oscillators.

II.4 Quantum Mechanics

The Hamiltonian operator in the normal mode representation is

[TABLE]

where $\hat{P}_{n}\equiv\hat{p}_{n}^{\prime}/\sqrt{m\hbar\omega_{n}}$, $Q_{n}\equiv q_{n}^{\prime}\sqrt{m\omega_{n}/\hbar}$, and $m\omega_{n}^{2}\equiv-\lambda\mu_{n}$. Since $\mu_{n}<0$ and $\lambda>0$, the frequencies are real.

For each mode, the energy eigenvalues areMessiah61 ; Merzbacher70

[TABLE]

and the corresponding energy eigenfunctions are the Hermite functions,Messiah61 ; Merzbacher70 ; WikiQHO

[TABLE]

where $\mathrm{H}_{l_{n}}(Q)$ is the Hermite polynomial of degree $l_{n}$. The eigenfunctions of the system have product form,

[TABLE]

and the energy in such an eigenstate is of course $E_{\bf l}=\sum_{n=1}^{N}E_{n,l_{n}}$. Since the mode amplitudes are a function of the positions, the energy eigenfunction will often instead be written $\phi_{\bf l}({\bf q})$.

III Results

III.1 Computational Details

The computational implementation of the eigenvalue and eigenvector calculations is obvious and need not be detailed here. There is perhaps one interesting algorithmic challenge for the quantum statistical mechanical aspects of the problem.

Because of the Maxwell-Boltzmann factor, it was useful to order the possible energy states beginning with the lowest. The $N$ mode frequencies $\omega_{n}$ are known, and one has to order the possible sets of quanta, ${\bf l}=\{l_{1},l_{2},\ldots,l_{N}\}$, in terms of their energy, $E_{{\bf l}(l)}\leq E_{{\bf l}(l+1)}$, $l=1,2,\ldots$. This was done by initially creating $l^{\mathrm{max}}$ states from the lowest frequency mode only, ${\bf l}(l)=\{l,0,0,\ldots,0\}$, $l=1,2,\ldots,l^{\mathrm{max}}$. Then, one quantum of the next higher frequency mode was added and the corresponding state inserted in order, bumping up the higher energy states, and discarding the highest. This continued until the insertion point reached $l^{\mathrm{max}}$, at which point the cycle was repeated for the next higher frequency.

III.2 Units

In order to make contact with a real physical system, in the results below parameters for a Lennard-Jones model of neon may be used. These are a mass $m=3.35\times 10^{-26}$ kg, a separation of zero force $r_{\mathrm{e}}=3.13\times 10^{-10}$ m, and a Lennard-Jones well-depth $\varepsilon=4.93\times 10^{-22}$ J. Sciver12

The Lennard-Jones pair potential, $u_{\mathrm{LJ}}(r)=\varepsilon[(r_{\mathrm{e}}/r)^{12}-2(r_{\mathrm{e}}/r)^{6}]$, allows one to define the Lennard-Jones frequency, $\omega_{\mathrm{LJ}}=\sqrt{72\varepsilon/(mr_{\mathrm{e}}^{2})}$, which equals $3.29\times 10^{12}$ Hz for neon. This is derived from the curvature of the Lennard-Jones potential at the zero force separation. The results below are mainly expressed in units of the Lennard-Jones zero force separation $r_{\mathrm{e}}$, and the Lennard-Jones frequency $\omega_{\mathrm{LJ}}$.

The thermal wave length, $\Lambda_{\mathrm{th}}\equiv\sqrt{2\pi\hbar^{2}/(mk_{\mathrm{B}}T)}$, can be used as a guide to the importance of symmetrization effects. (Here $\hbar$ is Planck’s constant, $k_{\mathrm{B}}$ is Boltzmann’s constant, and the convenient inverse temperature is $\beta\equiv 1/k_{\mathrm{B}}T$.) When the thermal wave length is comparable to, or larger than, the lattice spacing $\Delta_{q}\equiv\overline{q}_{j+1}-\overline{q}_{j}$, symmetrization effects should be measurable. They are also measurable when the spring constants are small.

III.3 Numerical Results

Figure 1 shows the positions of the particles in the eigenvectors of the potential energy matrix for the four particle harmonic crystal. For the lowest frequency mode, the particles are all displaced in the same direction (in phase), with the central two particles having twice the amplitude of the outer two. For the second lowest frequency, the middle pair of particles are out of phase. For the second highest frequency, both outer pairs of particles are out of phase. And finally, for the highest frequency, all three consecutive pairs of particles are out of phase.

Figure 2 shows the quantized energy levels for the harmonic crystal for several sets of parameters. Initially the energy of each level increases rapidly with level number. But for large level numbers, the rate of change of energy with energy level is sub-linear. Or to put it another way, the density of energy states increases with increasing energy. For $l\gtrsim 10$, the data is well fitted by $E_{l}\propto l^{1/N}$, as can be seen in the inset to the figure.

Figure 3 shows the average energy for a canonical equilibrium system as a function of inverse temperature. The data tests the dependence of the average on the number of levels used, $l_{\mathrm{max}}$. One can see that for $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 0.8$, the results for all $l_{\mathrm{max}}\geq 500$ are indistinguishable. For $\beta\hbar\omega_{\mathrm{LJ}}\lesssim 0.5$, there is a discernable difference between $l_{\mathrm{max}}=$ 2,000 and 5,000, and this difference increases with increasing temperature (decreasing inverse temperature).

The figure also shows the energy of the ground state, $E_{1}=\sum_{n=1}^{N}\hbar\omega_{n}/2$. It can be seen that for these parameters, excited states make negligible contribution when the temperature is lower than $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 2$.

The present crystal has $2N$ harmonic modes in the classical Hamiltonian ($N$ in the potential energy, and $N$ in the kinetic energy). Hence by the equipartition theorem, the average energy is $\langle{\cal H}\rangle_{\mathrm{cl}}=N/\beta$, which is the dotted curve in Fig. 3. One can see that this lies increasingly below the quantum results as the temperature decreases (inverse temperature increases), and it lies increasingly above the quantum results as the temperature increases. In the present case there is a region, $0.3\lesssim\beta\hbar\omega_{\mathrm{LJ}}\lesssim 0.4$, in which the quantum results for $l_{\mathrm{max}}=5,000$ coincide with the classical result. The data suggests what one knows to be true: the classical result must be the limiting result at high temperatures, but an increasing number of energy levels contribute to the average as the temperature increases. Hence for fixed $l_{\mathrm{max}}$, there is always a temperature above which the quantum results become inaccurate.

One can draw the important conclusion from this figure that starting from the exact quantum approach is a very inefficient way of obtain the classical result, which is the exact high temperature limit.

Figure 4 shows the average energy for several values of the interparticle spring constant $\lambda$. In the low temperature limit, $\beta\rightarrow\infty$, the average energy is the ground state energy, $E_{0}=\sum_{n=1}^{N}\hbar\omega_{n}/2$. It can be seen in Fig. 4 that for $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim$ 3–9 ($\lambda/m\omega_{\mathrm{LJ}}^{2}=$ 0.5–0.02) the system may be considered to be in the ground state.

Figure 5 shows the density profiles corresponding to the same cases as the preceding figure at the temperature $\beta\hbar\omega_{\mathrm{LJ}}=2$. The profiles are normalized to integrate to $N=4$. In general the peaks and troughs in the profiles indicate that the particles are mainly localized to their respective lattice positions. At the highest interparticle spring constant shown, the density is zero between lattice points, which means that there is little overlap between the particles. Conversely, at the lowest coupling shown, the density peaks are much broader and there is a high probability of finding a particle between the lattice positions.

IV Conclusion

This paper has obtained analytic expressions for the energy eigenvalues and eigenfunctions of a realistic system composed of interacting particles, namely a one dimensional harmonic crystal. Although the results are mainly practical and utilitarian in nature, there is one generic conceptual point that emerges: because wave function symmetrization does not apply to crystals (i.e. solids in which the particles may be distinguished by their attachment to fixed points in space), the consequent quantized vibrational modes, otherwise known as phonons, are not bosons. It would be more accurate to call them ‘non-bosonic phonons’ rather than ‘fermionic phonons’, as in the title to this paper. This is a general conceptual point that holds beyond the one-dimensional harmonic crystal treated here.

As mentioned the bulk of the paper is concerned with numerical results for the one-dimensional harmonic crystal. A number of quantitative conclusions can be drawn from these, such as the degree to which the average energy increases, and the structure in the density profile decreases, with increasing temperature.

The advantage of being able to easily explore the parameter space with the analytic model can be illustrated by comparison with the exact but numerical studies of Hernando and Vaníček,Hernando13 who, with praiseworthy effort, obtained the first 50 energy eigenvalues of a one-dimensional Lennard-Jones system. Of course their Lennard-Jones particles confined by a weak harmonic potential is not the same as the present particles interacting with linear springs and confined by fixed wall particles. Nevertheless, in so far as the separation between the Lennard-Jones particles is approximately $r_{\mathrm{e}}$,Hernando13 the Lennard-Jones frequency $\omega_{\mathrm{LJ}}=\sqrt{72\varepsilon/(mr_{\mathrm{e}}^{2})}$ seems the appropriate energy scale to use in qualitatively comparing the two systems.

One point of interest is that Hernando and VaníčekHernando13 obtained the average density profile for 4 particles for several temperatures, namely $\beta\hbar\omega_{\mathrm{LJ}}=$ 6.8, 11.9, and 47.5 in the present units. They found marked changes in the profiles, with the highest temperature case losing almost completely the density peaks present at the lower temperatures. This means that the excited states contribute significantly to their higher temperature results. The fact that Fig. 3 for the present harmonic crystal at $\rho=r_{\mathrm{e}}^{-1}$ shows no evidence of excited states for $\beta\hbar\omega_{\mathrm{LJ}}\gtrsim 2$, and significant change from the ground state only for $\beta\hbar\omega_{\mathrm{LJ}}\lesssim 0.5$, suggests that the inverse temperature for the harmonic crystal should be multiplied by a factor of 10–20 for the purposes of comparison with the Lennard-Jones system. Given that the 50 excited states obtained by Hernando and VaníčekHernando13 make a significant difference to their density profiles, it is questionable whether or not 50 energy eigenvalues are actually enough to accurately describe the system at those particular temperatures. In particular, the present Fig. 3 shows that 5,000 energy eigenvalues are necessary for the harmonic crystal for $\beta\hbar\omega_{\mathrm{LJ}}\lesssim 0.5$, which, if scaled by 10–20, would be comparable to the higher temperatures used in the earlier study. Comparing the two different systems is obviously of questionable validity, but it is intriguing that the mean field classical phase space calculations of the present authorAttard18c were in better agreement with the lower temperature results of Hernando and VaníčekHernando13 than with their higher temperature ones.

In any case the larger lesson from the present exact calculations is that in the terrestrial regime, where quantum effects are comparable to, or a perturbation on the classical result, it is necessary to obtain a prohibitively large number of energy eigenvalues for even quite a small system. One can conclude from this that it is better to treat terrestrial condensed matter systems as a quantum perturbation of a classical system, rather than as a fully quantum system. It is hoped quantitatively to confirm or refute this conclusion with explicit classical phase space calculations for the present model in a future publication.

Appendix A Symmetrization for Spin-Position Factorization

This appendix is independent of the one-dimensional harmonic crystal explored in the text.

The set of commuting dynamical variables for one particle $j$ may be taken to be ${\bf x}_{j}=\{{\bf q}_{j},\sigma_{j}\}$, where $\sigma_{j}\in\{-S,-S+1,\ldots,S\}$ is the $z$-component of the spin of particle $j$. (See Messiah (1961) §14.1 or Merzbacher (1970) §20.5.) Messiah61 ; Merzbacher70 Note that here $\sigma$ is not a spin operator or a Pauli spin matrix. Label the $2S+1$ spin eigenstates of particle $j$ by $s_{j}\in\{-S,-S+1,\ldots,S\}$, and the spin basis function by $\alpha_{s_{j}}(\sigma_{j})=\delta_{s_{j},\sigma_{j}}$. Note that this is not a spinor. For $N$ particles, ${\bm{\sigma}}\equiv\{\sigma_{1},\sigma_{2},\ldots,\sigma_{N}\}$, and similarly for ${\bf s}$ and ${\bf q}$, and the basis functions for spin space are $\alpha_{\bf s}({\bm{\sigma}})=\delta_{{\bf s},{\bm{\sigma}}}=\prod_{j=1}^{N}\delta_{s_{j},\sigma_{j}}$. The states ${\bf n}$ may be single or multi-particle states.

An unsymmetrized wave function $\psi({\bf x})$ in general has symmetrized form

[TABLE]

with the symmetrization factor being

[TABLE]

Alternatively, one can expand the wave function in terms of spin-position basis functions,

[TABLE]

Here and below the $\alpha_{\bf s}$ will be called the spin basis functions, and the $\phi_{\bf n}$ will be called the position basis functions. This nomenclature favors brevity over precision; better might be, for example, the basis functions for spin and position space, respectively. Instead of position one could use the momentum representation. It will often prove useful to choose the $\phi_{\bf n}$ to be energy eigenfunctions. The $\alpha_{{\bf s}}({\bm{\sigma}})$ and the $\phi_{{\bf n}}({\bf q})$ form a complete orthonormal set.

The symmetrization of any wave function in spin-position space can be accomplished by using symmetrized basis functions in its expansion. The latter are given by

[TABLE]

This is exact.

In place of this exact symmetrization, one can invoke an approximation that relies upon the factorization of the spin-position basis function into the sum of products of symmetrized position basis functions and symmetrized spin basis functions, namely

[TABLE]

The merits or otherwise of this approximation are discussed at the end of §A.1 below. Although the left hand side is normalized by the overall factor of $\surd\tilde{\chi}^{\pm}_{{\bf n},{\bf s}}$, the individual factors on the right hand side are not normalized. This is essential to the correct formulation of the ansatz and is emphasized by the tilde. The un-normalized symmetrized basis functions are

[TABLE]

and

[TABLE]

The $\alpha_{{\bf s}}({\bm{\sigma}})$ and the $\phi_{{\bf n}}({\bf q})$ are normalized. The $\surd N!$ here is an immaterial constant that is convenient but not essential. Respective symmetrization factors for use below may be defined in terms of these,

[TABLE]

It is essential that these symmetrization factors are not used to normalize the $\tilde{\phi}_{\bf n}^{\pm}({\bf q})$ and the $\tilde{\alpha}_{\bf s}^{\pm}({\bm{\sigma}})$ individually. The reason for this is that without individual normalization, all terms in the approximation for $\Phi_{{\bf n},{\bf s}}^{\pm}({\bf x})$ have equal weight when written as permutation sums. If the individual factors were normalized, then the symmetric factors would have a different weight to the antisymmetric factors, and since different products have two, one, or zero of each of these, individual terms in the permutation sum would have different weights.

This point also explains why the two products in each line of the approximation are simply added together symmetrically. In principle, the two products could be superposed with a relative phase factor and with a relative probability factor. The reason they aren’t is that ultimately the expression is meant to approximate a permutation sum in which all terms have equal weight (apart from the $(-1)^{p}$ for fermions). These two points will be taken up in the next subsubsection, and with the concrete example of two particles following that.

The approximation is sufficient to ensure the symmetrization of the spin-position basis functions,

[TABLE]

as can be confirmed by inspection.

The overall normalization factor for the symmetrized basis function $\Phi_{{\bf n},{\bf s}}^{\pm}({\bf x})$ for the approximation is

[TABLE]

It should be noted that in certain states the Fermi exclusion principle means that $\tilde{\alpha}_{\bf s}^{-}({\bm{\sigma}})$ or $\tilde{\phi}_{\bf n}^{-}({\bf q})$ vanish. In such states $\chi^{-}_{\bf s}$ and $\chi^{-}_{\bf n}$ also respectively vanish.

A.1 Comparison of Exact and Approximate Forms

One can label the $N!$ permutations of the dynamical variables ${\bm{\sigma}}_{P}$ and ${\bf q}_{P}$, $P=1,2,\ldots,N!$, in such a way that the permutation has the same parity as its label. The exact result is

[TABLE]

The approximate expression may be written

[TABLE]

In the final equality, the single sum is over terms where the permutation of the spins is the same as that of the positions. This is the same as the exact result. (Evidently, if this were the only contribution, then the normalization constants would be related as $\tilde{\chi}^{\pm}_{{\bf n},{\bf s}}=4\chi^{\pm}_{{\bf n},{\bf s}}/N!$.) The double sum, where the permutation of the spins differs from that of the particles, does not appear in the exact result and is unphysical in the sense that it disassociates the spin and position of each particle. It may be noted that only permutations of spin and position with the same parity make a non-zero contribution to this double sum. In the case of $N=2$ no such terms exist, so the double sum is zero, and the approximation is exact in this case (see also next).

Notice that the agreement of part of the approximation with the exact formulation depends upon using the un-normalized symmetrized basis functions, $\tilde{\alpha}^{\pm}_{\bf s}({\bm{\sigma}})$ and $\tilde{\phi}^{\pm}_{\bf n}({\bf q})$, and upon superposing the two terms without any phase factor or probability weight. As mentioned, the justification for this is that this procedure weights all terms in the permutation sum equally.

In the light of this analysis of the symmetrized basis functions, it is worth discussing whether or not the approximation Eq. (A.5) has any advantages over the exact form, Eq. (A.4). One could argue that the sum of products in the approximate form has a transparent interpretation that lends itself to the physical interpretation of symmetrization effects, and of physical phenomena such as Bose-Einstein condensation and Fermi exclusion. Indeed, the long-standing electronic orbital theory for optical spectra is predicated on this sum of products form, (which is exact for $N=2$; see also next). It could also be argued that there could be computational advantages to obtaining and storing the symmetrized spinless functions, and at a later stage combining them with the symmetrized spin functions for different values of $S$. Finally, it might be argued that exploring the properties of the symmetrized spinless functions leads directly to an understanding of the spatial (or momentum) localization of symmetrization effects, which is often missed in the formal treatment of symmetrization. Such localization shows that in many important terrestrial cases the nearest neighbor dimers give the dominant contribution, and for such pairs the approximation is exact.

A.2 Simple Example, $N=2$

For two particles, $N=2$, the comparison of the exact formulation for symmetrization with the approximate form can be performed rather directly. With ${\bf x}^{\prime}\equiv\hat{\mathrm{P}}_{12}{\bf x}$, the exact result is

[TABLE]

The approximation gives

[TABLE]

This is the same as the exact form. (Evidently in this case $\tilde{\chi}^{\pm}_{{\bf n},{\bf s}}=2\chi^{\pm}_{{\bf n},{\bf s}}$.)

One can illustrate the approximation further by making direct contact with conventional electronic orbital theory for the simple case of two fermions, $N=2$, $S=1/2$. In this case the symmetrized, un-normalized spin basis functions are

[TABLE]

and

[TABLE]

This may be re-written to show the same and different state occupancies explicitly,

[TABLE]

Here and throughout, the complementary Kronecker delta is $\overline{\delta}_{j,k}\equiv 1-\delta_{j,k}$.

For the basis functions in position one can consider single particle states, $\phi_{\bf n}({\bf q})=\phi_{{\bf n}_{1}}({\bf q}_{1})\phi_{{\bf n}_{2}}({\bf q}_{2})$. In the case of electronic orbitals, typically ${\bf n}_{j}=\{n_{j},l_{j},m_{j}\}$. Then

[TABLE]

and

[TABLE]

For use shortly, it follows that

[TABLE]

The overall normalization factor is

[TABLE]

For fermions, this vanishes if both are in the same state. Since the normalization factor appears in the denominator as the square root, the symmetrized spin-position basis function, $\Phi^{-}_{{\bf n},{\bf s}}({\bf x})$, also vanishes in this case.

Putting these together, the anti-symmetrized two fermion wave function is

[TABLE]

For $s_{1}=s_{2}$ this is

[TABLE]

This vanishes if ${\bf n}_{1}={\bf n}_{2}$, which is just the Fermi exclusion principle.

For ${\bf n}_{1}={\bf n}_{2}$ one has

[TABLE]

This vanishes if $s_{1}=s_{2}$, which is again the Fermi exclusion principle. The term in brackets is equivalent to the so-called singlet state.

For ${\bf n}_{1}\neq{\bf n}_{2}$ one has

[TABLE]

These terms are equivalent to the so-called triplet state (the first term represents the $++$ and $--$ states, and the second and third terms are a superposition of the $+-$ and $-+$ state).