Suppose I have $(\theta, \phi)\in [0, \pi)\times [0, 2\pi)$ as spherical coordinates for the 2-sphere. That is, I parameterize a point on the sphere according to \begin{align} x_1 &= \sin\theta\cos\phi \\ x_2 &= \sin\theta\sin\phi \\ x_3 &= \cos\theta. \end{align} Suppose that I put box constraints on my spherical coordinates so that I actually have $\theta_l \leq \theta\leq \theta_u$ and $\phi_l\leq\phi\leq\phi_u$. I would like to find the maximal possible distance between two such points in this box. That is, what are $(\theta, \phi)$ and $(\theta',\phi')$ such that the distance between $(x_1,x_2,x_3)$ and $(x_1',x_2',x_3')$ is maximized, subject to the constraint $\theta_l\leq \theta,\theta'\leq\theta_u$ and $\phi_l\leq\phi,\phi'\leq\phi_u$.
When it comes to the notion of distance, it seems natural to measure distance by the great circle distance on the sphere, $d((x_1, x_2,x_3), (x_1',x_2',x_3')) = \cos^{-1}(x_1x_2' + x_2x_2' + x_3x_3')$. However, I just end up with a cumbersome expression in trigonometric functions when I take this approach. Therefore, I would also be willing to use a somewhat less natural measure of distance such as \begin{align} d((x_1, x_2,x_3), (x_1',x_2',x_3')) = \sqrt{(x_1-x_1')^2 + (x_2-x_2')^2 + (x_3-x_3')^2} \end{align} for the sake of simplicity.
My intuition says that the maximal distance should be found on the boundary of this box in spherical coordinates. I have not been able to prove this however, for either distance. Does anyone have any suggestions?
Your conjecture is partly almost correct. :)
Consider a narrow band around the equator (so "all phis, but theta near pi/2, or, in more familiar terms, latitudes from 5 degrees N to 5 degrees S"). Clearly the max distance here is 2, with that distance achieved by the ends of a diameter. Naturally, two opposite points on the equator are ideal choices. But neither point would be on the boundary. If we tilt our diameter a little, we can put one end at 5N, and the other end at 5S and still achieve distance 2, with both points on the boundary.
But now...let's make our band run from 5N to 4S. The antipode of any point on the upper circle is NOT on the lower circle, and vice versa. So for this (degenerate) case, your conjecture is wrong --- you get no PAIR of boundary points achieving maximal euclidean distance.
But using "all phis" seems like a cheap trick, so let's let longitude run from 179W to 179E (i.e., $\theta$ goes from about $0.1$ to about $3.0$. This gives us two new boundary arcs. And if we go to a point that's at 4N on one of these arcs, its antipode will be on the lower circle (4S) at longitude near $0$.
So in the degenerate case of "all longitudes", your claim is (maybe) false --- it's not clear to me what the "boundary of this box in spherical coordinates" is exactly! --- but in other cases...it sure looks true.
(I've adjusted the answer above by swapping the roles of $\theta$ and $\phi$, because I didn't realize that you were using what I call the "physics parameterization", the "math one" being where $\theta$ is longitude and $\phi$ is latitude.)
I'm now going to make another adjustment, and define $\alpha = \theta - \frac{\pi}{2}$, and I'll do everything in $(\phi, \alpha)$ coordinates. This makes $\alpha$ (which I'll call 'latitude') run from $-\pi/2$ to $\pi/2$, which is easier for my explanation.
So let me prove your theorem.
We're looking at rectangles $0 \le \phi_L \le \phi \le \phi_U < 2\pi$. By subtracting $\phi_L$ from everything (which amounts to rotating the sphere on its N/S axis!) we can reduce to the case $$ 0 \le \phi \le K < 2\pi. $$ What about $\alpha$? It's in some range $$ \pi/2 \le \alpha_L \le \alpha \le \alpha_U \le \pi/2, $$ and here there's no easy simplification.
The boundary of this square (on the sphere) is made up of partial latitude circles (lat latitudes $\alpha_L$ and $\alpha_U$, and partial longitude arcs (at longitudes $0$ and $K$). This picture shows an example (in purple), where $\alpha_L$ is something like $-1.5$, and $\alpha_U$ is something like $+0.9$. $K$ is about $2.8$, say, so the two latitude circles go almost all the way around the globe.
In the case where $\alpha_L$ or $\alpha_U$ is $\pm \pi/2$, the lower or upper partial circle becomes a single dot at the north or south pole. If both are at their limits, the boundary becomes a pair of longitude arcs from pole to pole. In this situation, we see that the two poles are on the boundary, and achieve the maximal distance (namely 2), and we're done.
In what follows, I'll assume that one of the two is not at its limit, and by symmetry, I'll make that $\alpha_U$...so from now on, $\alpha_U < \pi/2$, while $\alpha_L$ could be $-pi/2$.
Look again at that figure, with the purple boundary. I've drawn another set, the "antipode" of the purple curve, i.e., for each purple point, I've drawn the antipodal point in green. You can see the purple and green sets intersect at a point in the lower right, marked with a red dot; that means, by symmetry, that they also intersect at the upper left (also marked in red). These two red points again are distance two apart, so they are boundary points that achieve the maximum possible distance.
This is the only insight we need to handle the case where $K \ge \pi$. Let's prove that.
Suppose that $K \ge \pi$. Then the four boundary arcs are, in $(\phi,\alpha)$ coordinates,
\begin{align} A:& (0, \alpha) & \alpha_L \le \alpha \le \alpha_U < \pi/2\\ B:& (K, \alpha) & \alpha_L \le \alpha \le \alpha_U < \pi/2\\ C:& (\phi, \alpha_L) & 0 \le \phi \le K < 2\pi\\ D:& (\phi, \alpha_U) & 0 \le \phi \le K < 2\pi \end{align}
We split into cases.
Case 1: $\alpha_L$ and $\alpha U$ have opposite signs
Case 2: $\alpha_L$ and $\alpha U$ have the same sign (which, by symmetry, we may assume is $+$).
For Case 1, one of the two numbers $|\alpha_L|$ and $\alpha_U$ is as large or larger than the other. By symmetry again, we maysuppose that $|\alpha_L| \ge |\alpha_U|$, i.e., that as in the picture, the lower circle is smaller than the upper. Because we are in Case 1, we have that $\alpha_L < 0 < \alpha_U$ (again as in the picture: the "rectangle" reaches both above and below the equator.)
Then the point $(\phi, \alpha) = (\pi, \alpha_U)$ lies on arc $D$, while the antipodal point, $(0, -\alpha_U)$ lies on arc $A$. (These are the two red dots in the figure). The distance between them is 2, which is the largest possible, so in Case 1, the maximum Euclidean distance is achieved on the boundary.
Case 2 can be handled by the same method as the $K < \pi$ case, which I'll address now.
Second situation: $K < \pi$ or $0 < \alpha_L < \alpha_U < \pi/2$.
I'll characterize both of these by observing that the "rectangle" on the sphere, in both cases, contains no antipodal points, so the max distance between points in the rectangles is less than two. The key characteristic, however, is that the rectangle contains no antipodal points.
Lemma: if the closed subset $A$ of the sphere contains no antipodal points, and $P, Q\in A$, with $P \in interior(A)$, then there's a point $P' \in A$ with $d(P', Q) > d(P, Q)$.
Corollary: the greatest distance between pairs of points of $A$ occurs between a pair of boundary points.
NB: this greatest distance exists because it's the maximum of a continuous function on a compact set (namely $A \times A \subset S^2 \times S^2$).
So all that remains is to prove the lemma.
Let me start by observing that if we let $f(X) = d(X, S)$, where $S$ is the south pole, and $X \in S^2$, then this real-valued function has nonzero gradient at every point of $S^2$ except the north and the south pole, which you can readily prove by explicit computation. The gradient direction, at each point, is the direction pointing "north", which is hardly surprising: to get more distant from the south pole, you head north! Stated differently, for any $X$ except the north and south pole, if you travel along the great circle arc from the south pole through $X$, the function $f$ is strictly increasing.
Generalizing, for any $A \in S^2$, the "distance from $A$" function, defined on $S^2$, has nonzero gradient at every point $X$ except $\pm A$, for any $A \in S^2$, and as you traverse the great circle arc from $A$ through $X$, this distance-from-$A$ function is strictly increasing until you reach the antipode of $A$.
Lemma: if the closed subset $A$ of the sphere contains no antipodal points, and $P, Q\in A$, with $P \in interior(A)$, then there's a point $P' \in A$ with $d(P', Q) > d(P, Q)$.
Proof: Because $A$ has no antipodal pairs,$P$ is not antipodal to $Q$. The distance-from-$Q$ function is strictly increasing on the great circle arc from $Q$ through $P$. If we write $\mathbf v = P - [P\cdot Q] Q$, then $v$ is perpendicular to the ray from the origin to $Q$, and nonzero, so $\mathbf u = \frac{1}{\|v\|}v$ is a unit vector orthogonal to that ray. Letting $S = O + \mathbf u$, where $O$ is the origin, and $c = \arccos(P \cdot Q)$, we can describe the great circle arc from $Q$ to $P$ by $$ \gamma(t) = \cos(ct) Q + \sin(ct) S. $$ This gives $\gamma(0) = Q$ and $\gamma(1) = P$.
Now $P$ is in the interior of $A$, so there's some open disk $D$ such that $P \in D \subset A$. If we let $$ E = \{t \in \Bbb R \mid \gamma(t) \in D\}, $$ then $E$, being the preimage of an open set under a continuous function, must be open. So there's some open interval $I = (1-a, 1 + b)$ ($a, b > 0$) around $t = 1$ contained in $E$. Let's look at $\gamma(1 + \frac{b}2)$. Now $1 + \frac{b}2$ is in the interval $I$, hence in $E$, hence $P' = \gamma( 1+ \frac{b}2)$ is in $D$, hence in $A$.
But we know that "distance from $Q$" is increasing along the great-circle arc, i.e., that $f(t) = d(Q, \gamma(t))$ is increasing as a function of $t$. Hence $f(P')$ is greater than $f(P)$, and we're done.
I apologize for this proof not being a good deal prettier, but it's the best I can do in limited time.