I'm assuming a uniform continuous distribution, assuming the probability distributions of the two points are independent, and assuming you mean "close to orthogonal" rather than strictly orthogonal.
It is true that "all locations on the sphere are equally probable candidates for both points." Assuming independent distributions,
all locations are equally probable for the second point even given a known position for the first point.
The point is that the part of the sphere close to the "equator" that you randomly selected when you placed the first point
is much larger than the part that is close to that point or to the antipodal point.
If you go to a high enough dimension, most of the sphere will be "close" to the "equator," for whatever your definition of "close" is.
The remarks above, of course, amount merely to a restatement of the claim in slightly different terms.
Better would be an actual proof.
Suppose we have a hypersphere of $n$ dimensions.
To simplify the calculations, let the radius of the sphere be $1.$
The total volume inside this sphere is therefore
$$
V_n(1) = \frac{\pi^{n/2}}{\Gamma\left(\frac n2 + 1\right)}.
$$
We assume a uniform distribution over the volume inside the sphere. That is, when we pick a random point according to this distribution, the probability that the point lands inside a particular region of hyperspace inside the sphere is proportional to the volume of that region.
Now arbitrarily pick a number $q$ with $0 < q \leq 1$ and consider the probability that the distance $r$ between the center of the hypersphere and the random point will be less than $q.$
This probability is the volume inside the hypersphere of radius $q$ divided by the total volume inside the hypersphere, namely,
$$
P(r \leq q) = \frac{V_n(q)}{V_n(1)}
= \frac{\frac{\pi^{n/2}}{\Gamma\left(\frac n2 + 1\right)} q^n}
{\frac{\pi^{n/2}}{\Gamma\left(\frac n2 + 1\right)}}
= q^n.
$$
Now suppose $n$ is very large, for example, $n = 10000,$ and let's see what is the probability that the random point lands within a distance $q = 0.999$ from the center of the hypersphere. This is
$$
P(r \leq 0.999) = 0.999^{10000} \approx 0.000045173.
$$
So the random point will land in this region less than $0.005\%$ of the time.
The other $99.995\%$ of the time, the point lands in the thin hyperspherical shell within just $0.001$ units from the surface.
This doesn't happen because the point was somehow attracted to the surface of the hypersphere or because we skewed the distribution to something non-uniform to make it close to the surface. This is a distribution that is strictly uniform by volume, with no other influences. But it almost always lands very close to the surface simply because that's where almost all the volume is.
High dimensions are weird!
Now let's randomly choose two points and measure the angle between them.
In order to measure the angle, after the points have already been chosen
(so that we can no longer unduly influence the random placement of the points), we choose one of the points as the pole of a set of spherical coordinates.
We do this only so that we can use the "latitude" in that coordinate system to measure the angle between the two points.
The chance that the angle will be between $\phi_1$ and $\phi_2$
(where these are both angles between $0$ and $180$ degrees)
is then proportional to how much of the sphere lies between the lines of latitude at angles $\phi_1$ and $\phi_2$ from the pole.
A general formula for this probability in high dimensions is rather ugly, I think. It is related to the question
what is the surface area of a cap on a hypersphere?
The probability density, however, is dealt with in the answer to
Distribution of an angle between a random and fixed unit-length $n$-vectors.
If the angle $\phi$ is measure in radians, the density over the range
$0 \leq \phi\leq \pi$ is
$$
f(\phi)=\frac{\sin^{n-2}\phi}{\int_0^\pi\sin^{n-2}\theta\,\mathrm d\theta}.
$$
Now let's take our $10000$-dimensional hypersphere again and let's see what the probability is that the angle is within $1/50$ radian (a little more than one degree) from a right angle. This is
$$
P(\frac\pi2 - \frac1{50} \leq \phi \leq \frac\pi2 + \frac1{50})
= \frac{\int_{\pi/2 - 1/{50}}^{\pi/2 + 1/{50}}\sin^{9998}\phi\,\mathrm d\phi}{\int_0^\pi\sin^{9998}\theta\,\mathrm d\theta}
\approx 0.95449.
$$
So there is a better than $95\%$ chance that the angle will be between
$88.85$ and $91.15$ degrees.
I interpreted the problem as asking for a function that, for general $R$ and $h$, gives the distance between the two points for any reflection angle (i.e., not only for the case in which the beams are tangent), and allows to calculate the maximal reflection angle and the maximal distance. This function can then be used to calculate the maximal angle and the maximal distance in the specific scenario provided in the OP, with $R=6371$ Km and $h=500$ Km.
Let us consider the Earth circumference as represented by the circle $x^2+y^2=R^2$, with the center in the origin and where $R=6371$. We can place our object in $A(0,R+h)$ on the $y$-axis, representing a point that is $h$ Km up the Earth surface.
Now let us draw two lines passing through $A$, symmetric with respect to the $y$-axis and intersecting the circumference. For each line, let us consider the intersection point that is nearer to the $y$-axis. Let us call the two new points $B$ (in the first quadrant) and $C$ (in the second quadrant). These represent the two points on Earth surface.
Due to the symmetry of the construction, we can continue by analyzing only one of these two points, e.g. $C$. The equation of the line containing $AC$ can be written as $y=sx+R+h$, where $s$ is its positive slope. To determine where this line crosses the circumference, we can set
$$sx+R+h=\sqrt{R^2-x^2}$$
whose solutions are
$$x=\frac{-sR-sh \pm \sqrt{\Delta}}{s^2 + 1}$$
where $\Delta=s^2 R^2 - 2h R - h^2$.
As stated above, we are interested in the less negative solution for $x$, as it is that nearer to the $y$-axis. So we get that the $x$-coordinate of $C$ is
$$X_C=\frac{-sR-sh + \sqrt{\Delta}}{s^2 + 1}$$
and the $y$-coordinate is
$$Y_C=\frac{s(-sR-sh + \sqrt{\Delta})}{s^2 + 1}+R+h$$
As a result, the equation $y=tx$ of the line $OC$ has slope
$$t=\frac{Y_C}{X_C}x=s-\frac{(s^2+1)(R+h)}{sh+sR-\sqrt{\Delta}}$$
Now setting $\angle{BAC}=\alpha$ and $\angle{BOC}=\beta$, we have $s=\cot(\alpha/2)$ and $t=-\cot(\beta/2)$.
Thus, we get
$$ \cot(\beta/2) =\frac{(\cot^2(\alpha/2)+1)(R+h)}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}}
- \cot(\alpha/2)$$
and then
$$ \beta =2\cot^{-1}\left[\frac{(\cot^2(\alpha/2)+1)(R+h)}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}} - \cot(\alpha/2)\right] $$
So the length of the arc $D$ corresponding to $ \beta$, which is the distance along the spherical surface asked in the OP, is
$$ D =2R\cot^{-1}\left[\frac{(\cot^2(\alpha/2)+1)(R+h)}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}} - \cot(\alpha/2)\right] $$
where $\Delta=\cot^2(\alpha/2)R^2 - 2h R - h^2$.
The last equation can be simplified as
$$ D =2R\cot^{-1}\left[\frac{(R+h)+\sqrt{\Delta}}{\cot(\alpha/2)(R+h)-\sqrt{\Delta}} \right] $$
For example, for $\alpha=\pi/2$ and $h= (\sqrt{2}-1)R$, as expected we have $\Delta=0$ (this is the situation where $\alpha$ is a right angle and the light beams are tangent to the surface). In this case, $\beta$ is also a right angle and $D=\pi/2\,R$. Accordingly the formula above gives this result, as shown by WA here.
For any value of $R$ and $h$, the maximal angle $\alpha_{max}$ and the maximal distance $D_{max}$ (i.e., those obtained with the beams tangent to the surface) can be determined by considering the case in which $\Delta=0$. This case occurs when $\cot^2(\alpha/2)R^2 - 2h R - h^2=0$. Solving for $\alpha$ in the range $0 \leq \alpha \leq \pi$ we get
$$\alpha_{max} = 2 \cot^{-1}\left(\frac{\sqrt{h(h+2R)}}{R}\right)$$
Interestingly, when $\Delta=0$, the formula for the distance is considerably simplified, and by few calculations reduces to
$$ D_{max} =2R\cot^{-1}\left[\tan(\alpha/2)\right] $$
As shown here, in the specific scenario described in the OP, substituting $R=6371$ and $h=500$, we get
$$\alpha= 2 \cot^{-1}\left(\frac{10 \sqrt{66210}}{6371}\right) \approx 2.3739 \,\,\text{radians}$$
which corresponds to about $136$ degrees. Here is the plot of the distance $D$ (in Km) as a function of $\alpha$ (in radians) for $R=6371$ and $h=500$, as obtained by WA. The plot confirms the maximal real value of $\alpha$, concordant with the predicted value of $2.3739$. The blue and red lines indicate the real and imaginary part, respectively.
Lastly, from the simplified formula for the maximal distance, taking $R=6371$ and $\alpha=2.3739$, we get
$$D_{max}\approx 4891 \, \text{Km}$$
Best Answer
These spherical trigonometry problems can be hard on the brain. But in this case the answer is surprisingly simple. You ask:
We can turn this problem inside out, and ask:
So let us take two points on the equator, with coordinates $A=(0,0)$ (see Null Island) and $(0,X)$, where $X\le 180$. Any point $P$ on the sphere defines a partition into two hemispheres, by considering that point as a pole. We may assume that the longitude $\theta$ of $P$ is in the range $[90^\circ,270^\circ]$, because if not, we can replace $P$ by its antipodal point. And for such a point $P=(\theta,\phi)$ to separate the points $A$ and $B$ into distinct hemispheres requires $90^\circ\le\theta\le X+90^\circ$. The probability of this is clearly $X/180$.
Hence the answer to your question is $1-X/180$.