We develop a novel method of measuring the lensing distortion profiles of clusters with stacking the scaled amplitudes of background galaxy ellipticities as a function of the scaled centric radius according to the NFW prediction of each cluster, base
d on the assumption that the different clusters in a sample follow the universal NFW profile. First we demonstrate the feasibility of this method using both the analytical NFW model and simulated halos in high-resolution $N$-body simulations. We then apply, as a proof of concept, this method to the Subaru weak lensing data and the XMM/Chandra X-ray observables for a sample of 50 massive clusters in the redshift range $0.15le zle 0.3$, where their halo masses range over an order of magnitude. To estimate the NFW parameters of each cluster, we use the halo mass proxy relation of X-ray observables, based on either the hydrostatic equilibrium or the gas mass, and then infer the halo concentration from the model $c(M)$ relation. We evaluate a performance of the NFW scaling analysis by measuring the scatters of 50 cluster lensing profiles relative to the NFW predictions over a range of radii, $0.14le R/[h^{-1}{rm Mpc}]le 2.8$. We found a 4 - 6$sigma$ level evidence of the universal NFW profile in 50 clusters, for both the X-ray halo mass proxy relations, although the gas mass appears to be a better proxy of the underlying true mass. By comparing the measurements with the simulations of cluster lensing taking into account the statistical errors of intrinsic galaxy shapes in the Subaru data, we argue that additional halo mass errors or intrinsic scatters of $sigma_{ln M_{500c}}sim 0.2$ - $0.3$ could reconcile a difference between the measurements and the simulations.
