Heterogeneous cellular networks (HCNs) play a significant role in 5G networks, but have mainly been modeled and analyzed on a two-dimensional (2D) plane. In this paper, we model a three-dimensional (3D) two-tier HCN consisting of macro-cells and small-cells via a stochastic geometry approach, where the base stations (BSs) in each tier are assigned an appropriate antenna downtilt to enhance the downlink signal and suppress the inter-cell interference. Based on the 3D HCN model, we derive the downlink spatially-averaged coverage probability and area spectral efficiency (ASE) as functions of the BS antenna downtilts, BS heights, BS densities and cell-association bias. To facilitate fast numerical evaluation, we propose accurate approximations for the integral parts in the coverage probability expression. We then obtain the optimal BS antenna downtilts for both tiers that maximize the downlink spatially-averaged coverage probability by using partial derivative and the bisection method. Analytical and simulation results show that for given BS heights, densities of small-cell BSs (SBSs) and cell-association biases of SBSs, our optimized BS antenna downtilts of both tiers can significantly enhance the downlink spatially-averaged coverage probability and ASE in comparison with the fixed BS antenna downtilts network. Useful insights into the 3D deployment of HCNs considering BS antenna downtilts are obtained based on the analytical and simulation results.INDEX TERMS 3D antenna downtilt, base station height, coverage probability, cell-assciation bias, heterogeneous cellular networks