PHYSICAL REVIEW B 112, 075158 (2025)
Atomistic tight-binding Hartree-Fock calculations of multielectron configurations
in P-doped silicon devices: Wavefunction reshaping
Maicol A. Ochoa ,
1,2,*
Keyi Liu ,
1,3
Piotr Ró
˙
za
´
nski ,
4
Michal Zieli
´
nski ,
4
and Garnett W. Bryant
1,3,†
1
National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA
2
Department of Chemistry and Biochemistry, University of Maryland, College Park, Maryland 20742, USA
3
Joint Quantum Institute, University of Maryland, College Park, Maryland 20742, USA
4
Institute of Physics, Faculty of Physics, Astronomy and Informatics, Nicolaus Copernicus University,
ul. Grudziadzka 5, 87-100 Toru´n, Poland
(Received 4 March 2025; revised 18 July 2025; accepted 12 August 2025; published 28 August 2025)
Donor-based quantum devices in silicon are attractive platforms for universal quantum computing and analog
quantum simulations. The nearly atomic precision in dopant placement promises great control over the quantum
properties of these devices. We present atomistic calculations and a detailed analysis of many-electron states
in a single phosphorus atom and selected phosphorus dimers in silicon. Our self-consistent method involves
atomistic calculations of the electron energies utilizing representative tight-binding Hamiltonians, computations
of Coulomb and exchange integrals without any reference to an atomic orbital set, and solutions to the
associated Hartree-Fock equations. First, we assess the quality of our tight-binding Hartree-Fock protocol against
configuration-interaction calculations for two electrons in a single phosphorus atom, finding that our formalism
provides an accurate estimation of the electron-electron repulsion energy requiring smaller computational
boxes and self-consistent single-electron wavefunctions. Then, we compute charging and binding energies in
phosphorus dimers observing their variation as a function of impurity-impurity separation. Our calculations
predict an antiferromagnetic ground state for the two-electron system and a weakly bound three-electron state in
the range of separations considered. We rationalize these results in terms of the single-electron energies, charging
energies, and the wavefunction reshaping.
DOI: 10.1103/14lg-249q
I. INTRODUCTION
Silicon nanostructures that incorporate dopants forming
arrays with atomic precision are promising platforms for de-
veloping quantum simulators, quantum materials, and other
quantum technologies. Soon after the proposal by Kane [1]
of a silicon-based quantum computer encoding information in
the nuclear spin of
31
P donor atoms, and the implementation
proposed by Loss and DiVincenzo [2] based on the spin states
of coupled single-electron quantum dots, initial attempts to
fabricate impurity-based devices with the required atomic
precision appeared in the literature [3–7]. This, eventually,
led to the fabrication of the single-atom transistor [8] and
other devices comprising a few phosphorus atoms [9–13].
Recently, reports on the analog quantum simulations of the
many-body Su-Schrieffer-Heeger (SSH) model [14] and the
3 × 3 extended Fermi-Hubbard model [15] in dopant-based
silicon devices appeared.
Electronic structure calculations of bound electrons in
impurity-based devices in semiconductors require efficient
methods that can account for the device’s atomistic details,
including the impurity array confinement, the host material
electronic properties, potential energy variations along the
*
Contact author: maicol@nist.gov
†
Contact author: garnett.bryant@nist.gov
device resulting from external gates, and, importantly, the
electron-electron interactions. This characterization is critical
to understand the performance, operation, and potential of
these devices in applications such as quantum simulations of
the Fermi-Hubbard model. A typical dopant array may cover a
few nanometers, resulting in electronic structure calculations
that must incorporate the multimillion host-material atoms
surrounding the array. Large computational boxes extend-
ing beyond the array’s domain are also needed to account
for the extended structure of the host material and avoid
artificial distortions due to the computational boundary. At
present, such types of calculations are forbidden for the most
accurate electronic structure methods and are limited to a
few thousand atoms [16]. Thus it is common practice to
use semiempirical tight-bind methods (TB) to describe these
systems. For the case of phosphorus-atom arrays in silicon,
tight-binding calculations provide a complete picture of the
electron density of a single-electron in dopant arrays [17,18]
at the same time that it reproduces the relevant band gaps
and effective masses of the host material. One can solve the
coupled Poisson-Schrödinger equations [19,20] or implement
configuration-interaction (CI) calculations [21] starting from
the TB single-particle basis to obtain the energy and struc-
ture of multielectron states in the dopant array. As is the
case with most common algorithms for many-body quantum
methods, the CI computational complexity grows exponen-
tially with electron number and single-particle basis size,
2469-9950/2025/112(7)/075158(9) 075158-1 ©2025 American Physical Society
