  We utilize Bardeen's tunneling theory to calculate intra- and interorbital
hopping integrals between phosphorus donors in silicon using known orbital wave
functions. While the two-donor problem can be solved directly, the knowledge of
hoppings for various pairs of orbitals is essential for constructing
multi-orbital Hubbard models for chains and arrays of donors. To assure
applicability to long-range potentials, we rederive Bardeen's formula for the
matrix element without assuming non-overlapping potentials. Moreover, we find a
correction to the original expression allowing us to use it at short distances.
We also show that accurate calculation of the lowest donor-pair eigenstates is
possible based on these tunnel couplings, and we characterize the obtained
states. The results are in satisfactory quantitative agreement with those
obtained with the standard H\"uckel tight-binding method. The calculation
relies solely on the wave functions in the barrier region and does not
explicitly involve donor or lattice potentials, which has practical advantages.
We find that neglecting the central correction potential in the standard method
may lead to qualitatively incorrect results, while its explicit inclusion
raises severe numerical problems, as it is contained in a tiny volume. In
contrast, using wave functions obtained with this correction in the proposed
method does not raise such issues. Nominally, the computational cost of the
method is to calculate a double integral along the plane that separates donors.
For donor separation in directions where valley interference leads to
oscillatory behavior, additional averaging over the position of the integration
plane is needed. Despite this, the presented approach offers a competitive
computational cost as compared to the standard one. This work may be regarded
as a benchmark of a promising method for calculating hopping integrals in
lattice models.
