In this article, we analyze the signal-to-interference ratio (SIR) coverage probability of a multi-cell downlink system with random traffic arrivals. Based on a theoretical model, in which the base stations (BSs) and devices are deployed as independent Poisson point processes (PPPs), we first present the sufficient and necessary conditions of the ε-stable region for the queueing system. Then, by taking into account the impact of spatially queueing interactions among BSs, we focus on the steady status of the queues and introduce a two-queue-length approximated model for the system. Specifically, in the studied new model, the BSs are separated into two sets in the steady state: the long-queue BS set and the short-queue BS set. The locations of the former are modeled by a Neyman-Scott process and those of the latter are modeled as a residual hole process. By applying the first-order statistic approximation, we further approximate the hole process by a homogeneous PPP with the same density. To model the deployment of the BSs precisely, the related parameters in the approximated model are fitted. Using tools from stochastic geometry and queueing theory, we derive the SIR coverage probability in the steady state. An iterative algorithm is proposed to calculate the active probability of the BSs. To reduce the computational complexity, Beta distribution is applied to approximate the probability density function of the service rate in each iteration of the algorithm. Finally, the effect of the fitting parameters and the accuracy of our analysis are presented via Monte Carlo simulations. INDEX TERMS Temporal traffic dynamics, spatially interacting queues, Poisson point process, Neyman-Scott process, stochastic geometry.