Counting intersection numbers of closed geodesics on Shimura curves

Let $\Gamma\subseteq\text{PSL}(2, \mathbb{R})$ correspond to the group of units of norm $1$ in an Eichler order $\mathrm{O}$ of an indefinite quaternion algebra over $\mathbb{Q}$. Closed geodesics on $\Gamma\backslash\mathbb{H}$ correspond to optimal embeddings of real quadratic orders into $\mathrm{O}$. The weighted intersection numbers of pairs of these closed geodesics conjecturally relates to the work of Darmon-Vonk on a real quadratic analogue to the difference of singular moduli. In this paper, we study the total intersection number over all embeddings of a given pair of discriminants. We precisely describe the arithmetic of each intersection, and produce a formula for the total intersection. This formula is a real quadratic analogue of the work of Gross and Zagier on factorizing the difference of singular moduli. The results are fairly general, allowing for a large class of non-maximal Eichler orders, and non-fundamental/non-coprime discriminants. The paper ends with some explicit examples illustrating the results of the paper.