Heterogeneous treatment effects are the center of gravity in many modern causal inference applications. In this paper, we investigate the estimation and inference of heterogeneous treatment effects with precision in a general nonparametric setting. To this end, we enhance the classical $k$-nearest neighbor method with a simple algorithm, extend it to a distributional setting, and suggest the two-scale distributional nearest neighbors (DNN) estimator with reduced finite-sample bias. Our recipe is first to subsample the data and average the 1-nearest neighbor estimators from each subsample. With appropriately chosen subsampling scale, the resulting DNN estimator is proved to be asymptotically unbiased and normal under mild regularity conditions. We then proceed with combining DNN estimators with different subsampling scales to further reduce bias. Our theoretical results on the advantages of the new two-scale DNN framework are well supported by several Monte Carlo simulations. The newly suggested method is also applied to a real-life data set to study the heterogeneity of treatment effects of smoking on children's birth weights across mothers' ages.