Abstract. 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-Berezinskiȋ (KTB) transition and the PDF becomes progressively less skewed with increasing temperature until it closely approximates a Gaussian function above T ≈ 4π. Near to that temperature we find some evidence of a phase transition, although our observations appear to exclude a thermodynamic singularity.