Calculate the integral of the Slater determinant

Question

This is a Slater determinant: $$ s=\left|\begin{array}{ll} \psi_{1 s}\left(r_1\right) \alpha & \psi_{1 s}\left(r_1\right) \beta \\ \psi_{1 s}\left(r_2\right) \alpha & \psi_{1 s}\left(r_2\right) \beta \end{array}\right|=\psi_{1 s}\left(r_1\right) \alpha \cdot \psi_{1 s}\left(r_2\right) \beta-\psi_{1 s}\left(r_1\right) \beta \cdot \psi_{1 s}\left(r_2\right) \alpha $$

where:

• $ r_1=(x_1,y_1,z_1), \,r_2=(x_2,y_2,z_2) $
• $\psi_{1s}(r_i)=e^{-\frac{1}{2} \sqrt{x_i^2+y_i^2+z_i^2}} / (2 \sqrt{2} \pi)$
• $\alpha=\left[\begin{array}{l}1 \\ 0\end{array}\right] \quad \beta=\left[\begin{array}{l}0 \\ 1\end{array}\right]$

I need to calculate:

$$ \left\langle s \left| \frac{1}{r_{12}} \right| s \right\rangle =\iiint \iiint s \cdot \frac{1}{\sqrt{\left(x_1-x_2\right)^2+\left(y_1-y_2\right)^2+\left(z_1-z_2\right)^2}} \cdot s \, \mathrm{d} x_1 \mathrm{d} y_1 \mathrm{d} z_1 \mathrm{d} x_2 \mathrm{d} y_2 \mathrm{d} z_2 $$

In order to do this, I used the following:

he[n_, l_, m_] := Module[
  (*Base transformation*)
  {xyz = CoordinateTransform["Cartesian" -> "Spherical", {x, y, z}]},
  (*Refer to external resource functions*)
  ResourceFunction["HydrogenWavefunction"][{n, l, m}, 2, xyz]]

(* Wave function of 1s*)
psi1s = he[1, 0, 0]

(*Define the wave function of two electrons*)
psi1 = psi1s /. {x -> x1, y -> y1, z -> z1};
psi2 = psi1s /. {x -> x2, y -> y2, z -> z2};

(*Define the spin of two electrons*)
spinUp = {1, 0};
spinDown = {0, 1};

(*Create Slater determinant*)
slater = {
  {psi1*spinUp, psi1*spinDown},
  {psi2*spinUp, psi2*spinDown}
  }

slater // Det

I encountered a problem when creating Slater determinants.

Det::matsq: Argument {...} at position 1 is not a non-empty square matrix.

I don't know how to fix the problem and proceed with calculating the integral.

======== update ========

Here is my reasoning process

Through my reasoning, the following conclusion can be drawn:

$$ \begin{aligned} & \left\langle s\left|\frac{1}{\left|\vec{r}_{12}\right|}\right| s\right\rangle \\ = & \int_{\phi_2=0}^{\phi_2=2 \pi} \int_{\phi_1=0}^{\phi_1=2 \pi} \int_{\theta_2 =0}^{\theta_2=\pi} \int_{\theta_1=0}^{\theta_1=\pi} \int_{\gamma_2=0}^{r_2=\infty} \int_{r_1=0}^{r_1=\infty}\left(\psi_{1 s}\left(r_1\right)\right)^2 \left(\psi_{1 s}\left(r_2\right)\right)^2 \cdot \frac{1}{\left|\vec{r}_{12}\right|} d r_1 d r_2 d \theta_1 d \theta_2 d \phi_1 d \phi_2 \end{aligned} $$

where:

• $$ \begin{array}{l} \psi_{1s}\left(r_1\right)=e^{-r_1 / 2} /(2 \sqrt{2} \pi) \\ \psi_{1s}\left(r_2\right)=e^{-r_2/2} /(2 \sqrt{2} \pi) \end{array} $$

• $$ \frac{1}{\left|\overrightarrow{r_{12}}\right|}=\frac{1}{\sqrt{r_1^2+r_2^2-2 r_1 \cdot r_2\left(\sin \theta_1 \cdot \sin \theta_2 \cdot \cos \left(\phi_1-\phi_2\right)+\cos \theta_1 \cos \theta_2\right)}} $$

It is not difficult to compute these integrals using Wolfram Language, as shown below

he[n_, l_, m_] :=
 ResourceFunction["HydrogenWavefunction"][{n, l, m}, 
  2, {r, \[Theta], \[Phi]}]

psi1s = he[1, 0, 0]

\[Psi]1 = 
 psi1s /. {r -> r1, \[Theta] -> \[Theta]1, \[Phi] -> \[Phi]1}
\[Psi]2 = psi1s /. {r -> r2, \[Theta] -> \[Theta]2, \[Phi] -> \[Phi]2}

eqj = \[Psi]1*\[Psi]2*1/
   Sqrt[(r1^2 + r2^2 - 
      2*r1*r2*(Sin[\[Theta]1]*Sin[\[Theta]2]*Cos[\[Phi]1 - \[Phi]2] + 
         Cos[\[Theta]1]*Cos[\[Theta]2]))]*\[Psi]1*\[Psi]2

j = NIntegrate[eqj,
  {r1, 0, Infinity},
  {\[Theta]1, 0, Pi},
  {\[Phi]1, 0, 2*Pi},
  {r2, 0, Infinity},
  {\[Theta]2, 0, Pi},
  {\[Phi]2, 0, 2*Pi}]

0.945383

It is not difficult to compute these integrals using Wolfram Language, but I want to use Wolfram Language to complete the reasoning process on this draft. What should I do?