The mean value of $y \log{y}$ over the ordinates of the CM points

Let $g : \Gamma \backslash \mathbb{H} \to \mathbb{C}$ be any bounded continuous function. Duke's theorem states that \[\frac{1}{h(D)} \sum_{A \in \mathrm{Cl}_K} g(z_A) = \frac{1}{\mathrm{vol}(\Gamma \backslash \mathbb{H})} \int_{\Gamma \backslash \mathbb{H}} g(z) \, d\mu(z) + o_g(1)\] as $D \to -\infty$ through negative fundamental discriminants, where $h(D) = \# \mathrm{Cl}_K$ the class number of $K = \mathbb{Q}(\sqrt{D})$, $\mathrm{Cl}_K$ the class group, and $d\mu(z) = \frac{dx \, dy}{y^2}$ the $\mathrm{SL}_2(\mathbb{R})$-invariant measure on the upper half-plane $\mathbb{H}$.

The proof uses the spectral decomposition \[g(z) = \frac{\langle g, 1\rangle}{\mathrm{vol}(\Gamma \backslash \mathbb{H})} + \sum_{f \in \mathcal{B}_0(\Gamma)} \langle g, f \rangle f(z) + \frac{1}{4\pi} \int_{-\infty}^{\infty} \left\langle g, E\left(\cdot, \frac{1}{2} + it\right)\right\rangle E\left(z, \frac{1}{2} + it\right) \, dt,\] where $\mathcal{B}_0(\Gamma)$ denotes an orthonormal basis of Hecke-Maass cusp forms and $E(z,s)$ denotes the real-analytic Eisenstein series.

We want to take $g(z) = y \log y$, but this is of course not square-integrable with respect to $d\mu(z)$. The trick is to let $g_s(z) = y^s \log y$ and $G_s(z) = g_s(z) - \frac{\partial}{\partial s} E(z,s)$ with $\Re(s) > 1/2$ and $s \neq 1$. Then $G_s(z)$ is square-integrable, with \[g_s(z) = \frac{\partial}{\partial s} E(z,s) + \frac{\langle G_s, 1\rangle}{\mathrm{vol}(\Gamma \backslash \mathbb{H})} + \sum_{f \in \mathcal{B}_0(\Gamma)} \langle G_s, f \rangle f(z) + \frac{1}{4\pi} \int_{-\infty}^{\infty} \left\langle G_s, E\left(\cdot, \frac{1}{2} + it\right)\right\rangle E\left(z, \frac{1}{2} + it\right) \, dt.\]

It is clear that $\langle G_s, f\rangle = \langle g_s, f\rangle$. Using Zagier's method of renormalisation of integrals, one can show that \[\langle G_s, 1\rangle = \int_{\sqrt{3}/2}^{1} y^{s - 2} \log y \left(1 - 2\sqrt{1 - y^2}\right) \, dy + \frac{1}{(s - 1)^2},\] and that \[\left\langle G_s, E\left(\cdot, \frac{1}{2} + it\right)\right\rangle = 0.\] (Both of these are somewhat nontrivial to show, however.) It follows by analytic continuation that \[g(z) = \lim_{s \to 1} \left(\frac{\partial}{\partial s} E(z,s) + \frac{1}{\mathrm{vol}(\Gamma \backslash \mathbb{H}) (s - 1)^2}\right) + C + \sum_{f \in \mathcal{B}_0(\Gamma)} \langle g, f \rangle f(z),\]

where \[C = \frac{1}{\mathrm{vol}(\Gamma \backslash \mathbb{H})} \int_{\sqrt{3}/2}^{1} y^{-1} \log y \left(1 - 2\sqrt{1 - y^2}\right) \, dy \approx -0.00181698.\]

Now we let $z = z_A$ and sum over all $A \in \mathrm{Cl}_K$. Then the sum over $f \in \mathcal{B}_0(\Gamma)$ is $O((-D)^{1/2 - \delta})$ for some $\delta > 0$ via the period formula for Hecke-Maass cusp forms at Heegner points and subconvexity bounds for $L(1/2,f \otimes \chi_D)$. Next, we have that \[\sum_{A \in \mathrm{Cl}_K} E(z_A,s) = \frac{w_K}{2} \frac{\Lambda_K(s)}{\Lambda(2s)},\] where $w_K = \# \mathcal{O}_{K,\mathrm{tors}}^{\times}$ (which is equal to $2$ for $D \leq -4$), $\Lambda_K(s) = (2\pi)^{-s} \Gamma(s) \zeta(s) L(s,\chi_D)$, and $\Lambda(2s) = \pi^{-s} \Gamma(s) \zeta(2s)$. We may rewrite this as \[\frac{w_K \sqrt{-D}}{4} \frac{\exp\left((s - 1) \log \frac{\sqrt{-D}}{2}\right) \zeta(s) L(s,\chi_D)}{\zeta(2s)}.\] This has a simple pole at $s = 1$ with residue \[\frac{3}{\pi} \frac{w_K \sqrt{-D} L(1,\chi_D)}{2\pi} = \frac{h(D)}{\mathrm{vol}(\Gamma \backslash \mathbb{H})}\] by the class number formula. It follows that if $C' = C'(D)$ denotes the coefficient of $s - 1$ in the Laurent expansion of this, then \[\sum_{A \in \mathrm{Cl}_K} g(z_A) = C' + C h(D) + O((-D)^{1/2 -\delta}).\] The main term from $C'$ is \[\frac{3}{\pi} \frac{w_K \sqrt{-D} L(1,\chi_D)}{2\pi} \frac{1}{2} \left(\log \frac{\sqrt{-D}}{2}\right)^2 = \frac{h(D)}{\mathrm{vol}(\Gamma \backslash \mathbb{H})} \frac{1}{2} \left(\log \frac{\sqrt{-D}}{2}\right)^2.\] To prove that the remaining terms in $C'$ are smaller requires bounds of the form \[\frac{L'}{L}(1,\chi_D) = o(\log (-D)), \quad \frac{L''}{L}(1,\chi_D) = o((\log (-D))^2).\] Such bounds certainly follow from the generalised Riemann hypothesis, though I don't believe they are known unconditionally.

---

I should mention that this method works more generally for $g(y) = y (\log y)^m$ with $m$ a nonnegative integer, in which case you need to take $m$-th derivatives with respect to $s$, or even more generally for $g(y) = y^s (\log y)^m$ with $s \neq 1$, in which case you no longer need to take limits, but $h(D)$ will no longer naturally appear in the main term, so the answer may fluctuate wildly depending on the size of $L^{(k)}(s,\chi_D)/L(1,\chi_D)$ for $0 \leq k \leq m + 1$.

Can we do this by reducing it to a sum over ideals?

By Duke's theorem, we can easily deal with the low $y$ terms. Thus we can adjust the boundary condition to $-a < b \leq a < \sqrt{- D /3 }$.

Instead of CM points, we are now working simply with pairs of an oriented lattice of discriminant $-D$ and an primitive lattice point of norm $a < \sqrt{- D /3 }$ with the bijection that sets the fixed lattice point to $(1,0)$ and sets $(0,1)$ to be the unique other generator with positive orientation where $b$ is in the correct range.

In other words, these are $\mathcal O_{\mathbb Q(\sqrt{-D})}$-modules with a primitive element generating a submodule of index $a < \sqrt{- D /3 }$ or, dually, ideals of $\mathcal O_{\mathbb Q(\sqrt{-D})}$ of norm $a < \sqrt{-D/3}$ that are not divisible by any natural number.

So you are looking at something like the sum of $$\frac{ \sqrt{D}}{2|I|} \log \left( \frac{ \sqrt{D}}{2|I|} \right) = \frac{\sqrt{D}}{2|I|} \log \left( \frac{\sqrt{D}}{2} \right) - \frac{ \sqrt{D}}{2|I|} \log |I| $$ over all indivisible ideals $I$ of $\mathcal O_{\mathbb Q(\sqrt{-D})}$ satisfying $|I| < \sqrt{ -D /3}$.

It should be possible to evaluate these terms by contour-shifting the generating functions $ \frac{ \zeta_{ \mathbb Q(\sqrt{-D})} (s) }{ \zeta(2s)}$ and $\frac{d}{ds} \frac{ \zeta_{ \mathbb Q(\sqrt{-D})} (s) }{ \zeta(2s)}$ respectively, no? The leading terms should come from the pole of the zeta function at $1$, which when divided by the class number should give the desired $(\log k)^2/2$ factor.