The burst approximation is a widely used technique to simplify stochastic gene expression models. However, the dynamics and analytical properties of the protein number distribution in gene expression models under the burst approximation are barely studied. In this study, we propose and systematically analyze surrogate models with multiple gene states and arbitrary burst size distributions. An analytical time-dependent solution to the chemical master equation is derived and then exploited in two directions. Theoretically, several fine properties of the protein number distribution are established using functional analysis. For geometrically distributed burst sizes, the distribution is dominated by a scaled negative binomial distribution, and is light-tailed in certain parameter regimes. Computationally, we develop efficient algorithms in three settings, enabling fast calculation of the protein number distribution. Furthermore, the approximation error relative to full gene expression models is estimated in terms of low-order moments of the distribution, thereby clarifying the validity of the burst approximation.