The question below originates from another post on MathOverflow and a (successful) attempt to prove the Chudnovsky algorithm with modular forms arising from quaternion algebra.
Motivation: As is stated in the former post, a left ideal $I$ of a maximal order $\mathcal{O}$ in a quaternion algebra ramified at $p$ ($p$ is a prime) and $\infty$ can be used to construct modular form $\theta_I$ of weight 2 with Fricke eigenvalue $-1:$
$$\theta_I(\tau)=\sum_{x\in I}e^{2\pi i\tau\frac{N(x)}{N(I)}}$$ One can use the quotient of different linear combinations of theta functions $\theta_I$ to construct explicit equations of modular curves $X_0^{+}(p)$. It is of great interest to investigate the exact positions of the poles of certain quotients (or the zeros of certain linear combination of theta functions from quaternion algebra).
Experiment: Let $p=37$ (a case which is known to E. Hecke). Suppose the maximal order $\mathcal{O}$ in $A(37)$ is the order given by the Proposition 5.2 in A. Pizer's paper, then one can find the class number of left-$\mathcal{O}$ ideals in $A(37)$ is $3$, and there are two different theta functions associated to these ideals:
$$\theta_{I_i}(\tau)=\sum_{x\in\mathbb{Z}^4}q^{x^{T}M_ix},q=e^{\pi i\tau},i=1,2$$
while the $M_1$ and $M_2$ are two Gram matrices attached to two even lattices of dimension $4$:
$$M_1=\left(\begin{matrix}2 & 0 & 1 & 1\\ 0 & 4 & 1 &2\\ 1 & 1 & 10 & 1\\ 1 & 2 & 1 & 20 \end{matrix}\right),M_2=\left(\begin{matrix}4 & 1 & 0 & 1\\ 1 & 6 & 3 &1\\ 0 & 3 & 8 & 2\\ 1 & 1 & 2 & 10 \end{matrix}\right).$$
Then $\phi=(\theta_{I_1}-\theta_{I_2})/2$ is a cusp form on $\Gamma_0(p)$ with Fricke eigenvalue $-1$($\phi$ is such cusp form with zeros at cusps of the possible highest order).
We note that $$\Phi(\tau)=\phi(\tau)\prod_{k=0}^{p-1}\phi\left(\frac{\tau+k}{p}\right)$$ is a modular form of weight $2(p+1)$ on $\Gamma(1)$. As $\Gamma_0(37)$ has two elliptic point of order $2$ and two elliptic points of order $3$, then $$\varphi=\frac{\Phi}{E_4^4E_6^2\Delta^4}$$ is a modular function on $\Gamma(1)$ (everywhere holomorphic except a double pole at $\infty$), $E_4,E_6,\Delta$ are Eisenstein series and modular discriminant, respectively. It is possible to determine numerically the first few coefficients of the $q$-expansion of $\varphi$, and one can determine that $$\varphi(\tau)=(j(\tau)-66^3)^2,$$ which leads to the conclusion that $\phi(\tau)$ vanishes at $(\pm12+2i)/37$.
Experiment Results: One can extend these calculations to some other primes. I tried to construct modular forms of weight $2$ on $\Gamma_0(p)$ with Fricke eigenvalue $-1$ and zeros of highest possible order at the cusp of $\Gamma_0(p)$ for a few primes, and all the zeros of such cusp forms are quadratic irrationals. The results are listed in the table below.
Notations: We denote $\#0$ to be the highest possible order of zeros at the cusp of $\Gamma_0(p)$ for the linear combinations of theta functions. Let $-D,D>0$, $(a_i,b_i,c_i)$ be all the possible reduced binary quadratic forms with discriminant $-D$ and $\tau_{(i,-D)}=(-b_i+\sqrt{-D})/(2a_i)$. We denote a zero of those modular forms by $(\tau_{(i,-D)},k),k\in\mathbb{Z}$, since those modular forms vanish at $(\tau_{(i,-D)}+k)/p$.
\begin{array}{|c|r|c|c|c|c|c|c|} \hline p & \#0 & 1 & 2 & 3 & 4 & 5 & 6\\ \hline 37&1& (\tau_{(1,-16)},12) & & & & & \\ 43&2& (\tau_{(1,-12)},13)& & & & & \\ 67& 3& (\tau_{(1,-12)},8) & (\tau_{(1,-7)},12)& & & & \\ 73& 4& (\tau_{(1,-27)},26)& & & & & \\ 97& 6& (\tau_{(1,-12)},26)& & & & & \\ 163& 7& (\tau_{(1,-7)},26)& (\tau_{(1,-11)},29)& (\tau_{(1,-27)},13)& (\tau_{(1,-44)},57)& (\tau_{(2,-44)},35)& (\tau_{(3,-44)},90)\\ 193& 9& (\tau_{(1,-8)},34)& (\tau_{(1,-16)},31) & (\tau_{(1,-24)},90)& (\tau_{(2,-24)},45)& (\tau_{(1,-75)},37)& (\tau_{(2,-75)},77)\\ \hline \end{array}
Question: What is the reason that all the zeros of these modular forms are quadratic irrationals?