  Finite-size scaling is investigated in detail around the critical point in
the heavy-quark region of nonzero temperature QCD. Numerical simulations are
performed with large spatial volumes up to the aspect ratio $N_s/N_t=12$ at a
fixed lattice spacing with $N_t=4$. We show that the Binder cumulant and the
distribution function of the Polyakov loop follow the finite-size scaling in
the $Z(2)$ universality class for large spatial volumes with $N_s/N_t \ge 9$,
while, for $N_s/N_t \le 8$, the Binder cumulant becomes inconsistent with the
$Z(2)$ scaling. To realize the large-volume simulations in the heavy-quark
region, we adopt the hopping parameter expansion for the quark determinant: We
generate gauge configurations using the leading order action including the
Polyakov loop term for $N_t=4$, and incorporate the next-to-leading order
effects in the measurements by the multipoint reweighting method. We find that
the use of the leading-order configurations is crucially effective in
suppressing the overlapping problem in the reweighting and thus reducing the
statistical errors.
