The mean value of $y \log{y}$ over the ordinates of the CM points Let $-D < 0$ be a negative fundamental discriminant and let $y$ range over the values $y = y_Q = \frac{\sqrt{|D|}}{2a}$, as the values $(a,b,c)$ run through the reduced binary quadratic forms $Q = aX^2+bXY+cY^2$ of discriminant $b^2-4ac = -D$ (thus $\mathrm{gcd}(a,b,c) = 1$, $-a < b \leq a \leq c$). Those are the ordinates of the CM points $z_Q = \frac{b + \sqrt{-D}}{2a}$ in the standard fundamental domain for $\mathcal{H} / \mathrm{PSL}(2,\mathbb{Z})$. 
By Duke's theorem, the $z_Q$ are equidistributed in the measure $\frac{3}{\pi} \frac{dx \, dy}{y^2}$. The highest lying point has ordinate $k := \sqrt{|D|}/2$, and if $\varphi: [1,\infty) \to \mathbb{R}$ is any compactly supported or fast enough decaying function, it follows that the mean value of $\varphi(\frac{\pi}{3}y_Q)$ (over all reduced forms $Q$ if the given discriminant $-D$) is equal  to $\int_1^k \varphi(t) \frac{dt}{t^2} + o_{\varphi}(1)$, as $D \to \infty$. The asymptotic $\mathrm{Avg}_Q \varphi(\frac{\pi}{3}y_Q) \sim \int_1^k 
\varphi(t) \frac{dt}{t^2}$ turns out to be also satisfied for $\varphi(t) = t$ (assuming $L(s,\chi_D)$ satisfies the Riemann hypothesis),  for a very special reason: the Kronecker-Chowla-Selberg limit formula expresses 
$$
\mathrm{Avg}_Q \big(\frac{\pi}{3} y_Q - \log{y_Q} \big) = \log{k} + \frac{L'}{L}(1,\chi_D) + O(1),
$$
with an absolutely bounded $O(1)$ term.
On the other hand it is plain that $\mathrm{Avg}_Q \varphi(\frac{\pi}{3}y_Q) \sim \int_1^k 
\varphi(t) \frac{dt}{t^2}$ fails to hold with $\varphi(t) = t^{1+\epsilon}$, for any fixed $\epsilon > 0$: already the principal form (with $a = 1$) contributes an $\sim k^{\epsilon} \frac{k}{h}$ to this average ($h$ denoting the class number), a contribution that can be as big as $k^{\epsilon} \log{\log{k}}$.
I am interested in what happens for $\varphi(t) = t\log{t}$.
Question. Assume GRH, and possibly some other standard analytic hypotheses. What can be said about the mean value of $\frac{\pi}{3} y_Q \log{y_Q}$, asymptotically as $D \to \infty$?
Evidently this mean value lies somewhere between $\log{k}$ and $(\log{k})^2$. Should the limit $\lim_{D \to \infty} \mathrm{Avg}_Q(y_{Q}\log{y_Q}) \big/ (\log{k})^2$ even exist?
[Of course, the same question could be asked for the positive discriminants (indefinite binary quadratic forms and closed geodesics on the modular surface). ] 
 A: 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.
A: 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$.
