Let $Z_1 \in \mathbb{R}^d$ and $Z_2\in\mathbb{R}^d$ be two nonzero vectors. Let $X$ be a random vector distributed uniformly on the hypersphere $\mathbb{S}^{d-1}$. In the case $d=2$, the marginal distribution of the angle $\gamma_1$ between $X$ and $Z_1$ is uniform on $[0,2\pi)$, as is the marginal distribution of the angle $\gamma_2$ between $X$ and $Z_2$.
What can be said about the joint distribution of $(\gamma_1, \gamma_2)$? In particular, for the case of general $d$, can we write the probability density function?
Let us begin with the definition of the Dirichlet distribution of $(U_1,\ldots,U_n)$ with parameters $(a_0,a_1,\ldots,a_n)$ all positive. It has density $$Cx_1^{a_1-1}\ldots x_n^{a_n-1}(1-x_1-\cdots-x_n)^{a_0-1}.$$ It is concentrated on $x_i>0$ $i=0,\ldots,n$ where we have written $x_0=1-x_1-\cdots-x_n$ for short. The constant $C$ is defined by $$\frac{1}{C}=B(a_0,a_1,\ldots,a_n)=\frac{\Gamma(a_0)\ldots\ \Gamma(a_n)}{\Gamma(a_0+\cdots+a_n)}.$$ The moments are easy to compute since $$\mathbb{E}(U_0^{b_0}\ldots U_n^{b_n})=\frac{B(a_0+b_0,a_1+b_1,\ldots,a_n+b_n)}{B(a_0,a_1,\ldots,a_n)}.$$ Here $U_0=1-U_1-\cdots-U_n.$ A typical example of Dirichlet distribution is given independent gamma random variables $Y_0,\ldots, Y_n$ with the same scale parmeter and shape parameters $(a_0,a_1,\ldots,a_n)$ by the formula $$U_i=\frac{Y_i}{Y_0+\cdots+Y_n}.$$ For instance if $Z=(Z_0,\ldots,Z_n)$ are independent $N(0,1)$ then $U_i=\frac{Z_i^2}{\|Z\|^2}$ defines a Dirichlet distribution with parameters $(a_0,a_1,\ldots,a_n)=(1/2,\ldots,1/2).$
Student@: the link of that stuff with your problem is the fact that $P=\frac{Z}{\|Z\|}$ is uniformly distributed on the unit sphere $S_n$ of $\mathbb{R}^{n+1}.$ Let me change your notations: I will call $A$ and $B$ your two points $Z_1$ and $Z_2$ and I understand that you are interested in the joint distribution of $(\gamma_1,\gamma_2)$ where $\cos \gamma_1=\langle A,P\rangle, \cos \gamma_2=\langle B,P\rangle .$ with $0<\gamma_i<\pi.$ If $e_1,\ldots,e_n,e_0$ is the canonical orthonormal basis of $\mathbb{R}^{n+1}$ then without loss of generality we may assume that there exists a number $\theta$ such that $$A=\cos \theta\, e_1+\sin \theta\, e_2,\ B=\cos \theta\, e_1-\sin \theta\, e_2$$ leading to $$\langle A,P\rangle=\frac{\cos \theta\, Z_1+\sin \theta\, Z_2}{\|Z\|},\ \langle B,P\rangle=\frac{\cos \theta\, Z_1-\sin \theta\, Z_2}{\|Z\|}.$$ As you can see this is easier to consider the joint distribution of $$V_1=\frac{1}{2}(\langle A,P\rangle+\langle B,P\rangle)=\cos \theta \frac{Z_1}{\|Z\|}, V_2=\frac{1}{2}(\langle A,P\rangle-\langle B,P\rangle)=\sin \theta \frac{Z_2}{\|Z\|}$$ Now observe that $V_1,V_2$ are dependent but symmetric. Therefore the knowledge of their distribution is given by the knowledge of the distribution of their squares. Since $$\frac{V_1^2}{\cos^2 \theta}=\frac{Z_1^2}{\|Z\|^2}, \frac{V_2^2}{\sin^2 \theta}=\frac{Z_2^2}{\|Z\|^2}$$ is Dirichlet distributed with parameters $ (1/2,1/2, \frac{n-1}{2})$, your problem is solved.