  We study both analytically, using the renormalization group (RG) to two loop
order, and numerically, using an exact polynomial algorithm, the
disorder-induced glass phase of the two-dimensional XY model with quenched
random symmetry-breaking fields and without vortices. In the super-rough glassy
phase, i.e. below the critical temperature $T_c$, the disorder and thermally
averaged correlation function $B(r)$ of the phase field $\theta(x)$, $B(r) =
\bar{<[\theta(x) - \theta(x+ r) ]^2>}$ behaves, for $r \gg a$, as $B(r) \simeq
A(\tau) \ln^2 (r/a)$ where $r = |r|$ and $a$ is a microscopic length scale. We
derive the RG equations up to cubic order in $\tau = (T_c-T)/T_c$ and predict
the universal amplitude ${A}(\tau) = 2\tau^2-2\tau^3 + {\cal O}(\tau^4)$. The
universality of $A(\tau)$ results from nontrivial cancellations between
nonuniversal constants of RG equations. Using an exact polynomial algorithm on
an equivalent dimer version of the model we compute ${A}(\tau)$ numerically and
obtain a remarkable agreement with our analytical prediction, up to $\tau
\approx 0.5$.
