The $β$-model is a powerful tool for modeling large and sparse networks driven by degree heterogeneity, where many network models become infeasible due to computational challenge and network sparsity. However, existing estimation algorithms for $β$-model do not scale up. Also, theoretical understandings remain limited to dense networks. This paper brings several significant improvements over existing results to address the urgent needs of practice. We propose a new $\ell_2$-penalized MLE algorithm that can comfortably handle sparse networks of millions of nodes with much-improved memory parsimony. We establish the first rate-optimal error bounds and high-dimensional asymptotic normality results for $β$-models, under much weaker network sparsity assumptions than best existing results. Application of our method to large COVID-19 network data sets discovered meaningful results.