  We present a detailed investigation of the probability density function (PDF)
of order parameter fluctuations in the finite two-dimensional XY (2dXY) model.
In the low temperature critical phase of this model, the PDF approaches a
universal non-Gaussian limit distribution in the limit T-->0. Our analysis
resolves the question of temperature dependence of the PDF in this regime, for
which conflicting results have been reported. We show analytically that a weak
temperature dependence results from the inclusion of multiple loop graphs in a
previously-derived graphical expansion. This is confirmed by numerical
simulations on two controlled approximations to the 2dXY model: the Harmonic
and ``Harmonic XY'' models. The Harmonic model has no
Kosterlitz-Thouless-Berezinskii (KTB) transition and the PDF becomes
progressively less skewed with increasing temperature until it closely
approximates a Gaussian function above T ~ 4\pi. Near to that temperature we
find some evidence of a phase transition, although our observations appear to
exclude a thermodynamic singularity.
