  The zero-temperature, classical $XY$-model on an $L \times L$ square-lattice
is studied by exploring the distribution $\Phi_L(y)$ of its centered and
normalized magnetization $y$ in the large $L$ limit. An integral representation
of the cumulant generating function, known from earlier works, is used for the
numerical evaluation of $\Phi_L(y)$, and the limit distribution $\Phi_{L
\rightarrow \infty} (y) = \Phi_0(y)$ is obtained with high precision. The two
leading finite-size corrections $\Phi_L (y) -\Phi_0 (y) \approx a_1(L)\,
\Phi_1(y) + a_2(L)\,\Phi_2(y)$ are also extracted both from numerics and from
analytic calculations. We find that the amplitude $a_1(L)$ scales as
$\ln(L/L_0) /L^2$ and the shape correction function $\Phi_1 (y)$ can be
expressed through the low-order derivatives of the limit distribution, $\Phi_1
(y) = [\,y\, \Phi_0 (y) + \Phi'_0 (y)\,]'$. The second finite-size correction
has an amplitude $a_2(L)\propto 1/L^2$ and one finds that $a_2\,\Phi_2(y) \ll
a_1 \,\Phi_1(y)$ already for small system size ($L> 10$). We illustrate the
feasibility of observing the calculated finite-size corrections by performing
simulations of the $XY$-model at low temperatures, including $T = 0$.
