Abstract. The vortex contribution to the probability density function of longitudinal magnetization fluctuations is examined in finite 2D XY systems close to the Kosterlitz-Thouless-Berezinskii transition temperature. Within the temperature range studied their relevance is limited to rare fluctuations, where they increase the probability of events exceeding four standard deviations below the mean magnetization.The characterization of fluctuation statistics is a central problem in the study of critical phenomena, as the break down of Landau theory on approaching the critical point, implies a non-Gaussian distribution for order parameter fluctuations [1]. From renormalization group theory it is customary to think of critical phenomena divided into universality classes, characterized by the symmetry group of the order parameter and the spatial dimension. One would therefore expect critical fluctuation statistics to be determined essentially by the set of critical exponents (say β and ν for a regular critical point) which describe the scaling behavior of derivatives of the singular part of the free energy. Evidence from, e.g. the Ising model [2,3] and Potts models [4] suggests that this is indeed generally the case.We have, however recently considered an exception to this established phenomenology, the 2D XY model [5,6]. Here the low temperature phase consists of a line of critical points, with one independent critical exponent η(T ) = T /2πJ (with J the exchange constant) extending down to zero temperature and separated from the high temperature paramagnetic phase by the Kosterlitz-Thouless-Berezinskii (KTB) phase transition [7,8,9]. It is well established that for all temperatures below the KTB transition temperature, T KTB , renormalisation group flows are to a quadratic effective Hamiltonian [10,11], with the result that the asymptotic behaviour, in the limit of large system size is perfectly captured by a harmonic model. The advantage of such a simple model is that the probability density function, P (m), for fluctuations of the order parameter, m, can be calculated analytically, without using either renormalisation group or the scaling hypothesis. Surprisingly, we find that the form of distribution is independent of η along the whole line of critical points. This universal scaling function, plotted as the solid line in Fig. 1, arises when σP (m) is plotted against µ = (m − m )/σ, where m and σ are the mean and standard ‡ pcwh@ens-lyon.fr § sellitto@ictp.trieste.it