Exact solution of the many-body Schrödinger equation for few electrons

\label{sec:mbse}

In one-dimensional systems, the fully interacting Hamiltonian for \(N\) electrons has the form \[\label{eq:1dham} \hat{H}=\sum_{j=1}^N \left(-\frac{d^2}{dx_j^2}+v_{\rm ext}(x_j)\right)+\sum_{j<k}^N v_{\rm int}(x_j, x_k),\] where the interaction potential \(v_{\rm int}(x_j, x_k)\) is usually Coulombic, though the following discussion also applies for other types of interaction, including more than two-body ones. In 1D one often uses the soft Coulomb interaction \(1 / \sqrt{(x_j-x_k)^2+1}\), where a softening parameter (usually set to one) is introduced in order to avoid the divergence at \(x_j=x_k\), which is non-integrable in 1D.

Mathematically, the Hamiltonian (eq. \eqref{eq:1dham}) is equivalent to that of a single (and hence truly independent) electron in \(N\) dimensions, with the external potential \[v_{\rm ext}^{Nd}(x_1...x_N)=\sum_{j=1}^N v_{\rm ext}(x_j)+\sum_{j<k}^N v_{\rm int}(x_j, x_k).\] For small \(N\) it is numerically feasible to solve the \(N\)-dimensional Schrödinger equation \[\label{eq:SENd} \hat{H}\Psi_j(x_1...x_N)=E_j\Psi_j(x_1...x_N)\] which provides a spatial wave function for a single particle in \(N\) dimensions. This equivalence is not restricted to one-dimensional problems. One can generally map a problem of \(N\) electrons in \(d\) dimensions onto the problem of a single particle in \(Nd\) dimensions, or indeed a problem with multiple types of particles (e.g. electrons and protons) in \(d\) dimensions, in the same way.

What we exploit in Octopus is the basic machinery for solving the Schödinger equation in an arbitrary dimension, the spatial/grid bookkeeping, the ability to represent an arbitrary external potential, and the intrinsic parallelization. In order to keep our notation relatively simple, we will continue to discuss the case of an originally one-dimensional problem with \(N\) electrons. Grid-based solutions of the full Schrödinger equation are not new, and have been performed for many problems with either few electrons (in particular H\(_2\), D\(_2\) and H\(_2^+\)\cite{Ranitovic21012014,lein2002}) or model interactions \cite{luo2013}, including time-dependent cases \cite{ramsden2012}.

The time-dependent propagation of the Schrödinger equation can be carried out in the same spirit, since the Hamiltonian is given explicitly and each “single-particle orbital” represents a full state of the system. A laser or electric-field perturbation can also be applied, depending on the charge of each particle (given in the input), and taking care to apply the same effective field to each particle along the polarization direction of the field (in 1D, the diagonal of the hyper-cube).

Solving eq. (\ref{eq:SENd}) leaves the problem of constructing a wave function which satisfies the anti-symmetry properties of \(N\) electrons in one dimension. For fermions one needs to ensure that those spatial wave functions \(\Psi_j\) which are not the spatial part of a properly anti-symmetric wave function are removed as allowed solutions for the \(N\)-electron problem. A graphical representation of which wave functions are allowed is given by the Young diagrams (or tableaux) for permutation symmetries, where each electron is assigned a box, and the boxes are then stacked in columns and rows (for details see, for example, Ref. \cite{McWeeny}). Each box is labeled with a number from 1 to \(N\) such that the numbers increase from top to bottom and left to right. All possible decorated Young diagrams for three and four electrons are shown in Fig. \ref{fig:young}. Since there are two different spin states for electrons, our Young diagrams for the allowed spatial wave functions contain at most two columns. The diagram d) is not allowed for the wave function of three particles with spin \(1/2\), and diagrams k) to n) are not allowed for four particles. To connect a given wave function \(\Psi_j\) with a diagram one has to symmetrize the wave function according to the diagram. For example, for diagram b) one would perform the following operations on a function \(\Psi(x_1,x_2,x_3)\) \[\left[ \Psi(x_1,x_2,x_3)+\Psi(x_2,x_1,x_3)\right]- \left[\Psi(x_3,x_2,x_1)+\Psi(x_3,x_1,x_2)\right].\] Hence, one symmetrizes with respect to an interchange of the first two variables, because they appear in the same row of the Young diagram, and anti-symmetrizes with respect to the first and third variable, as they appear on the same column. We note that we are referring to the position of the variable in the list, not the index, and that symmetrization always comes before anti-symmetrization. At the end of these operations one calculates the norm of the resulting wave function. If it passes a certain threshold, by default set to \(10^{-5}\), one keeps the obtained function as a proper fermionic spatial part. If the norm is below the threshold, one continues with the next allowed diagram until either a norm larger than the threshold is found or all diagrams are used up. If a solution \(\Psi_j\) does not yield a norm above the threshold for any diagram it is removed since it corresponds to a wave function with only bosonic or other non-fermionic character. Generally, as the number of forbidden diagrams increases with \(N\), the number of wave functions that need to be removed also increases quickly with \(N\), in particular in the lowest part of the spectrum. The case of two electrons is specific, as all solutions of eq. (\ref{eq:SENd}) correspond to allowed fermionic wave functions: the symmetric ones to the singlet states and the anti-symmetric ones to the triplet states.