Let $\vec{x}_1$ and $\vec{x}_2$ be the two points. Let $r = |\vec{x}_1 - \vec{x}_2|$ be the distance between them. By elementary geometry, if you draw two circle of radius $r$ using these two points as center, the area of their intersection is given by $(\frac{2\pi}{3} - \frac{\sqrt{3}}{2})r^2$. Notice the picking of two points are independent, we have:
$$E\left[ \vec{x}_1 \cdot \vec{x}_2 \right] = E\left[\vec{x}_1\right] \cdot E\left[\vec{x}_2\right] = \vec{0} \cdot \vec{0} = 0$$
This implies
$$E\left[|\vec{x}_1 - \vec{x}_2|^2\right] = E\left[|\vec{x}_1|^2 + |\vec{x}_2|^2\right] = 2\frac{\int_0^R r^3 dr}{\int_0^R rdr} = R^2$$
As a result, the expected area of the intersection is $(\frac{2\pi}{3} - \frac{\sqrt{3}}{2})R^2$.
Update for those who are curious
Let $\mathscr{C}$ be the set of events such that the intersection contains the origin, then:
$$\begin{align}
\operatorname{Prob}\left[\,\mathscr{C} \right] &= \frac{2\pi + 3\sqrt{3}}{6\pi}\\
E\left[\,|\vec{x}_1 - \vec{x}_2|^2 : \mathscr{C}\right] &= \frac{20\pi + 21\sqrt{3}}{6(2\pi + 3\sqrt{3})}
\end{align}$$
and the expected area of intersection conditional to containing the center is given by:
$$\frac{(4\pi - 3\sqrt{3})(20\pi + 21\sqrt{3})}{36(2\pi + 3\sqrt{3})}$$
To evaluate $E\left[ \varphi(\vec{x}_1,\vec{x}_2) ) : \mathscr{C} \right]$ for any function $\varphi( \vec{x}_1, \vec{x}_2 )$ which is symmetric and rotational invariant w.r.t its argument, you need to compute an integral of the from:
$$\int_{\frac{\pi}{3}}^{\pi} \frac{d\theta}{\pi} \left[2\int_{0}^{R} \frac{2udu}{R^2} \left( \int_{\alpha(\theta)u}^{u} \frac{2vdv}{R^2} \phi( \vec{x}_1, \vec{x}_2 ) \right) \right] $$
where $u \ge v$ are $|\vec{x}_1|$ and $|\vec{x}_2|$ sorted in descending order.
$\theta$ is the angle between $\vec{x}_1$ and $\vec{x}_2$. The mysterious $\alpha(\theta)$ is $\max(2\cos(\theta),0)$ for $\theta \in [\frac{\pi}{3},\pi]$.
The integral is a big mess and I need a computer algebra system to crank that out.
I won't provide more details on this part not relevant to the main answer.
Rewrite your first hyperboloid equation as
$$
\sqrt{(x-x_1)^2+(y-y_1)^2+(z-z_1)^2} - d_{12}=
\sqrt{(x-x_2)^2+(y-y_2)^2+(z-z_2)^2}.
$$
Square both sides and rearrange to get
$$
\tag{1}
2x(x_2-x_1)+2y(y_2-y_1)+2z(z_2-z_1)
+d_1^2 - d_2^2 + d_{12}^2=
2d_{12}\sqrt{(x-x_1)^2+(y-y_1)^2+(z-z_1)^2},
$$
where $d_1^2=x_1^2+y_1^2+z_1^2$ and $d_2^2=x_2^2+y_2^2+z_2^2$.
Do the same for the other hyperboloid:
$$
\tag{2}
2x(x_3-x_1)+2y(y_3-y_1)+2z(z_3-z_1)
+d_1^2 - d_3^2 + d_{13}^2=
2d_{13}\sqrt{(x-x_1)^2+(y-y_1)^2+(z-z_1)^2},
$$
where $d_3^2=x_3^2+y_3^2+z_3^2$.
If you now multiply $(1)$ by $d_{13}$, $(2)$ by $d_{12}$, and subtract $(2)$ from $(1)$, you end up with a linear equation:
$$
\tag{3}
2x[d_{13}(x_2-x_1)-d_{12}(x_3-x_1)]+
2y[d_{13}(y_2-y_1)-d_{12}(y_3-y_1)]+
2z[d_{13}(z_2-z_1)-d_{12}(z_3-z_1)]
+d_{13}(d_1^2 - d_2^2 + d_{12}^2)
-d_{12}(d_1^2 - d_3^2 + d_{13}^2)= 0.
$$
This is the equation of a plane: the locus you need is then the intersection between this plane and either hyperboloid, that is a conic section (hyperbola or ellipse)
As you can see the equations are getting more and more complicated, so I will just sketch a possible strategy to parameterize that conic.
Square equation $(1)$ (or $(2)$, it doesn't matter) and simplify the result.
Insert into that equation the expression for $z$ obtained from $(3)$ (if equation $(3)$ doesn't contain $z$ just solve for $x$ or $y$ instead).
The equation you thus get is that of a conic section in plane $xy$, projection onto that plane of the original conic. You can find its center and axes with the usual formulas and thus find parametric equations $x(t)$ and $y(t)$ for its coordinates.
Plug $x(t)$ and $y(t)$ into the expression for $z$ obtained above, to get also z(t).
Best Answer
I'll calculate the result for two points randomly uniformly picked the unit disk; it scales with $R^2$.
Consider $p_1$ at distance $r_1$ from the origin. Where can $p_2$ lie such that the origin is closer to either point than the other point? For the origin to be closer to $p_1$ than $p_2$, we must have $p_2$ outside the circle around $p_1$ with radius $r_1$. For the origin to be closer to $p_2$ than $p_1$, we must have $p_2$ on the origin's side of the line equidistant from the origin and $p_2$.
Thus, given $p_1$ at $(0,r_1)$, the admissible region for $p_2$ is the difference of two circular segments, the segment of the circle around the origin with radius $1$ up to $y=r_1/2$ minus the segment of the circle around $p_1$ with radius $r_1$ up to $y=r_1/2$. The area whose expected value we seek is proportional to the squared distance $d^2$ between $p_1$ and $p_2$. Thus, to calculate the desired expected value, we need the area of a circular segment and the integral over a circular segment of the squared distance to a given point on its axis. In detail:
$$ \begin{align} A(r,a)&:=\int_{-r}^a\mathrm dy\int_{-\sqrt{r^2-y^2}}^{\sqrt{r^2-y^2}}\mathrm dx\;, \\ B(r,a,b)&:=\int_{-r}^a\mathrm dy\int_{-\sqrt{r^2-y^2}}^{\sqrt{r^2-y^2}}\mathrm dx\left(x^2+(y-b)^2\right)\;, \\ \langle d^2\rangle&=\frac{\displaystyle\int_0^1\mathrm dr_1r_1\left(B\left(1,\frac{r_1}2,r_1\right)-B\left(r_1,-\frac{r_1}2,0\right)\right)}{\displaystyle\int_0^1\mathrm dr_1r_1\left(A\left(1,\frac{r_1}2\right)-A\left(r_1,-\frac{r_1}2\right)\right)}\;. \end{align} $$
Some help from our electronic friends yields the integrals:
$$ \begin{align} A(r,a)=&a\sqrt{r^2-a^2}+r^2\arctan\frac a{\sqrt{r^2-a^2}}+\frac\pi2r^2\;, \\ B(r,a,b)=&\frac16\left( \left(6b^2a+8b\left(a^2-r^2\right)+r^2a+2a^3\right)\sqrt{r^2-a^2}\right. \\ &\left.+3r^2\left(r^2+2b^2\right)\left(\frac\pi2+\arctan\frac a{\sqrt{r^2-a^2}}\right)\right)\;. \end{align} $$
In the negative terms with $r=r_1$, the arctangent is
$$\arctan\frac{-r_1/2}{\sqrt{r_1^2-(r_1/2)^2}}=-\arctan\frac1{\sqrt3}=-\frac\pi6\;.$$
Then further integral crunching yields
$$ \begin{align} \langle d^2\rangle&=\frac{\displaystyle\frac1{288}(88\pi+81\sqrt3)-\frac1{288}\left(8\pi-3\sqrt3\right)}{\displaystyle\frac1{16}(4\pi+3\sqrt3)-\frac14\left(\frac\pi3-\frac{\sqrt3}4\right)} \\ &= \frac{\displaystyle\frac5{18}\pi+\frac7{24}\sqrt3}{\displaystyle\frac\pi6+\frac{\sqrt3}4} \\ &= \frac16\frac{20\pi+21\sqrt3}{2\pi+3\sqrt3} \\ &\approx1.44034\;, \end{align} $$
which is confirmed by numerical simulations. The average squared distance of two points uniformly randomly picked in the unit disk without conditions is $1$, so the condition that the lune determined by the points contain the origin increases the average squared distance considerably.
The area to be averaged is $2A(d,-d/2)=d^2(4\pi-3\sqrt3)/6$, so the expected area is
$$ \frac{4\pi-3\sqrt3}6\langle d^2\rangle=\frac1{36}\frac{(4\pi-3\sqrt3)(20\pi+21\sqrt3)}{2\pi+3\sqrt3}\approx1.76927\;. $$
The problem with three points seems considerably more involved.
[Edit:]
Some more details on the calculation, to respond to the comments.
The basic idea is that the conditional expectation value of $d^2$ is the integral of $d^2$ over the volume satisfying the condition, divided by the volume satisfying the condition. Because we're dividing two integrals, we don't need to worry about constant factors. If we were instead calculating the probability for the condition to be fulfilled, we'd have to multiply the denominator by $2\pi$ for the missing integration over $\phi_1$ and divide it by $\pi^2$ for normalizing the density, for an overall factor of $2/\pi$, so the probability for the condition to be fulfilled is
$$\frac2\pi\left(\frac\pi6+\frac{\sqrt3}4\right)=\frac13+\frac{\sqrt3}{2\pi}\approx0.609\;.$$
But since we're dividing two integrals, this factor $2/\pi$ cancels. The integration is performed in polar coordinates for $p_1$ and in Cartesian coordinates for $p_2$; the origin for the coordinates is the centre of the unit circle in case of the polar coordinates for $p_1$ and the centre of the respective circle with radius $r$ in case of the Cartesian coordinates for $p_2$, and those are two different circles in the two terms of the difference.
The difference is the difference of the integrals over two different circular segments (as discussed in the comments). Both segments have the line equidistant between $p_1$ and the centre of the unit circle as their linear boundary. The positive term is the integral over the lower segment of the unit circle $C(O,1)$ defined by that line, and the negative term is the integral over the lower segment of the circle $C(p_1,r_1)$ defined by that line, where $O$ is the origin and $C(p,r)$ denotes the circle with origin $p$ and radius $r$.
Since the segments both have the same form and just differ in their centre and radius, the calculation is organized into calculating integrals $A$ and $B$ with general radius $r$ and boundary at $a$ and then substituting specific values in the two terms of the difference. Since Cartesian coordinates relative to the centre of the respective circle are used, the $y$ coordinate $b$ of $p_1$ relative to the respective centre also needs to be taken into account in calculating the squared distance $d^2$ from $p_1$.
So in the integrals $A$ and $B$ over a circular segment ($A$ being the integral of $1$, i.e. the volume, and $B$ being the integral of $d^2$), the parameter $r$ is the radius of the circle, the parameter $a$ is the upper $y$ coordinate (relative to the circle's centre) at which the circular segment ends, and the parameter $b$ is the $y$ coordinate of $p_1$ (relative to the circle's centre).
For the positive term, which is for the circular segment of the unit circle $C(O,1)$, the radius $r$ is $1$, the $y$ coordinate $a$ of the bisector relative to the centre at $O$ is $r_1/2$, and the $y$ coordinate $b$ of $p_1$ relative to the centre at $O$ is $r_1$ (since $p_1$ is at $(0,r_1)$). For the negative term, which is for the circular segment of the circle $C(p_1,r_1)$, the radius $r$ is $r_1$, the $y$ coordinate $a$ of the bisector relative to the centre at $p_1$ is $-r_1/2$, and the $y$ coordinate $b$ of $p_1$ relative to the centre at $p_1$ is $0$.
Since we're integrating over the lower segments of the respective circles, the $y$ coordinate runs from its minimal value $-r$ to its value $a$ at the linear boundary, and the $x$ coordinate takes its entire range of values from $-\sqrt{r^2-y^2}$ to $\sqrt{r^2-y^2}$.
This is all for the integration over $p_2$ given $p_1$ fixed at $(0,r_1)$. Then we also have to integrate over $p_1$. The integration over the polar angle $\phi_1$ of $p_1$ just yields a constant factor $2\pi$ that cancels (see above), since it's only the relative orientation of $p_1$ and $p_2$ that matters, not the absolute orientation of $p_1$. So we can just integrate the result we got for $p_1$ at $(0,r_1$) over $r_1$, taking into account the Jacobian $r_1$ of the polar coordinates.