Calculus of Variations
Using latitude, $\phi$ and longitude, $\theta$, we have
$$
1=\cos^2(\phi)\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)^2+\left(\frac{\mathrm{d}\phi}{\mathrm{d}s}\right)^2\tag1
$$
Integrating $(1)$ with respect to $s$, taking the first variation, and integrating by parts, we get
$$
\scriptsize\delta s=\int\left[\left(\color{#C00}{2\sin(\phi)\cos(\phi)\,\frac{\mathrm{d}\theta}{\mathrm{d}s}\,\frac{\mathrm{d}\phi}{\mathrm{d}s}-\cos^2(\phi)\,\frac{\mathrm{d}^2\theta}{\mathrm{d}s^2}}\right)\delta\theta-\left(\color{#090}{\sin(\phi)\cos(\phi)\,\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)^2+\frac{\mathrm{d}^2\phi}{\mathrm{d}s^2}}\right)\delta\phi\right]\mathrm{d}s\tag2
$$
Orthogonality requires that a curve where we have $\delta s=0$ for any $\delta\theta$ and $\delta\phi$ must satisfy
$$
\color{#C00}{\frac{\mathrm{d}^2\theta}{\mathrm{d}s^2}=2\tan(\phi)\,\frac{\mathrm{d}\theta}{\mathrm{d}s}\,\frac{\mathrm{d}\phi}{\mathrm{d}s}}\tag3
$$
and
$$
\color{#090}{\frac{\mathrm{d}^2\phi}{\mathrm{d}s^2}=-\sin(\phi)\cos(\phi)\,\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)^2}\tag4
$$
Solving the Equations
$$
\begin{align}
\log\left(\frac{\mathrm{d}\theta}{\mathrm{d}s}\right)
&=2\log(\sec(\phi))+\log(\cos(\epsilon))\tag{5a}\\
\frac{\mathrm{d}\theta}{\mathrm{d}s}
&=\cos(\epsilon)\sec^2(\phi)\tag{5b}\\
\left(\frac{\mathrm{d}\phi}{\mathrm{d}s}\right)^2
&=1-\cos^2(\epsilon)\sec^2(\phi)\tag{5c}\\
s
&=\int\frac{\mathrm{d}\phi}{\sqrt{1-\cos^2(\epsilon)\sec^2(\phi)}}\tag{5d}\\
&=\int\frac{\mathrm{d}\sin(\phi)}{\sqrt{\cos^2(\phi)-\cos^2(\epsilon)}}\tag{5e}\\
&=\int\frac{\mathrm{d}\sin(\phi)}{\sqrt{\sin^2(\epsilon)-\sin^2(\phi)}}\tag{5f}\\
&=\sin^{-1}\left(\frac{\sin(\phi)}{\sin(\epsilon)}\right)+s_0\tag{5g}\\[9pt]
\end{align}
$$
Explanation:
$\text{(5a)}$: divide $(3)$ by $\frac{\mathrm{d}\theta}{\mathrm{d}s}$ and integrate
$\text{(5b)}$: remove the log
$\text{(5c)}$: apply $(1)$ to $\text{(5b)}$
$\text{(5d)}$: separate the differential equation and integrate
$\text{(5e)}$: multiply the integrand by $\frac{\cos(\phi)}{\cos(\phi)}$
$\text{(5f)}$: $\cos^2(x)=1-\sin^2(x)$
$\text{(5g)}$: evaluate the integral
Solving $\text{(5g)}$ for $\sin(\phi)$ gives
$$
\bbox[5px,border:2px solid #C0A000]{\sin(\phi)=\sin(s-s_0)\sin(\epsilon)}\tag6
$$
Furthermore,
$$
\begin{align}
\theta
&=\int\frac{\cos(\epsilon)\,\mathrm{d}(s-s_0)}{1-\sin^2(s-s_0)\sin^2(\epsilon)}\tag{7a}\\
&=\int\frac{\cos(\epsilon)\,\mathrm{d}\tan(s-s_0)}{\sec^2(s-s_0)-\tan^2(s-s_0)\sin^2(\epsilon)}\tag{7b}\\
&=\int\frac{\cos(\epsilon)\,\mathrm{d}\tan(s-s_0)}{1+\tan^2(s-s_0)\cos^2(\epsilon)}\tag{7c}\\[6pt]
&=\tan^{-1}(\tan(s-s_0)\cos(\epsilon))+\theta_0\tag{7d}
\end{align}
$$
Explanation:
$\text{(7a)}$: apply $\text{(5h)}$ to $\text{(5b)}$, separate the differential equation, and integrate
$\text{(7b)}$: multiply the integrand by $\frac{\sec^2(s-s_0)}{\sec^2(s-s_0)}$
$\text{(7c)}$: $\sec^2(x)=1+\tan^2(x)$
$\text{(7d)}$: evaluate the integral
Solving $\text{(7d)}$ for $\tan(\theta-\theta_0)$ gives
$$
\bbox[5px,border:2px solid #C0A000]{\tan(\theta-\theta_0)=\tan(s-s_0)\cos(\epsilon)}\tag8
$$
Compare with a Great Circle
Given the parameterization in $\mathbb{R}^3$ of a great circle
$$
(x,y,z)=(\cos(s-s_0),\sin(s-s_0)\cos(\epsilon),\sin(s-s_0)\sin(\epsilon))\tag9
$$
we get
$$
\sin(\phi)=z=\sin(s-s_0)\sin(\epsilon)\tag{10}
$$
and
$$
\tan(\theta-\theta_0)=\frac yx=\tan(s-s_0)\cos(\epsilon)\tag{11}
$$
Equation $(10)$ matches equation $(6)$ and equation $(11)$ matches equation $(8)$. Thus, the shortest path between two points on a sphere is an arc of a great circle.
Plotting the Solved Equations
Since $\phi\in\left[-\frac\pi2,\frac\pi2\right]$, we can get $\phi$ by applying $\sin^{-1}$ to $(6)$:
$$
\bbox[5px,border:2px solid #C0A000]{\phi=\sin^{-1}(\sin(\epsilon)\sin(s-s_0))}\tag{12}
$$
However, to get $\theta-\theta_0$, we cannot simply apply $\tan^{-1}$ to $(8)$. Instead, we note that
$$
\begin{align}
\tan((\theta-\theta_0)-(s-s_0))
&=\frac{\tan(\theta-\theta_0)-\tan(s-s_0)}{1+\tan(\theta-\theta_0)\tan(s-s_0)}\tag{13a}\\
&=\frac{\tan(s-s_0)(\cos(\epsilon)-1)}{1+\cos(\epsilon)\tan^2(s-s_0)}\tag{13b}
\end{align}
$$
Explanation:
$\text{(13a)}$: formula for the tangent of a difference
$\text{(13b)}$: apply $(8)$
Thus, we get
$$
\bbox[5px,border:2px solid #C0A000]{\theta-\theta_0=(s-s_0)-\tan^{-1}\left(\frac{(1-\cos(\epsilon))\tan(s-s_0)}{1+\cos(\epsilon)\tan^2(s-s_0)}\right)}\tag{14}
$$
Great circles with $\epsilon=\left\{\color{#DF0000}{0^{\large\circ}},\color{#FF8F00}{15^{\large\circ}},\color{#EFDF00}{30^{\large\circ}},\color{#00BF00}{45^{\large\circ}},\color{#0000FF}{60^{\large\circ}},\color{#8000FF}{75^{\large\circ}},\color{#CC00FF}{90^{\large\circ}}\right\}$ plotted using $(12)$ and $(14)$.
I think it is a little more straightforward if you work in Cartesian coordinates.
Take an earth of radius one, for simplicity.
Then we can take the longitudes as $\pm l$ without loss of generality. Then
the halfway point between the two points, projected to the surface of the earth, will have the highest latitude (equivalently the highest $z$ component).
The two points are $(\cos \phi \cos l, \cos \phi \sin l, \sin \phi)$ and
$(\cos \phi \cos l, -\cos \phi \sin l, \sin \phi)$, and the mid point is
$(\cos \phi \cos l, 0, \sin \phi)$. To project to the surface, we divide by
the norm, to get a $z$ component (on our unit earth) of
$\sin \delta = {\sin \phi \over \sqrt{ (\cos \phi \cos l)^2 + \sin^2 \phi}} = { \tan \phi \over \sqrt{\cos^2l+\tan^2 \phi} }$.
To obtain the $\arctan$, we first need $\cos \delta$, which we get from
$\cos \delta = \sqrt{1 -\sin^2 \delta} = { \cos l \over \sqrt{ \cos^2l+\tan^2 \phi} } $, and hence
$\tan \delta = {\tan \phi \over \cos l}$, which is the desired result.
Best Answer
You are right. The angle between the two plane is ${90}°$. Here is two ways to see it.
First, the point $P$ is the one closest to the North pole. Since great circle are the straigh line of spherical geometry, the shortest distance between a point and a line is obtain with a perpendicular.
Second, the angle between the plane is given by the angle between their normal vector.
Let consider the sphere centered at the origin and the poles on the $z$ axis. Rotate the sphere so the point $P$ is on the $yz$-plane. Note: the point $Q$ is also on the $yz$-plane since $P$ and $Q$ are the ends of the same diameter.
The normal of the plane passing thru $P$, $Q$ and the two poles is on the $x$ axis. E.g. the normal could be the vector starting at $O$ pointing toward were the equtorial plane meet the oblique plane.
The normal vector of the oblique plane passing thru $P$ and $Q$ is in the $yz$-plane. The normal could be the vector starting at $O$ pointing toward the mid point between $P$ and $Q$ on the great circle passing by the poles.
The two vectors are perpendicular.
Hope it helps