The energy and flux budget (EFB) closure theory for a passive scalar (non-buoyant and non-inertial particles or gaseous admixtures) is developed for stably stratified turbulence. The physical background of the EFB turbulence closures is based on the budget equations for the turbulent kinetic and potential energies and turbulent fluxes of momentum and buoyancy as well as the turbulent flux of particles. The EFB turbulence closure is designed for stratified geophysical flows from neutral to very stable stratification, and it implies that turbulence is maintained by the velocity shear at any stratification. In a steady-state, expressions for the turbulent flux of the passive scalar and the anisotropic non-symmetric turbulent diffusion tensor are derived, and universal flux Richardson number dependencies of the components of this tensor are obtained. The diagonal component in the vertical direction of the turbulent diffusion tensor is suppressed by strong stratification, while the diagonal components in the horizontal directions are not suppressed, but they are dominant in comparison with the other components of the turbulent diffusion tensor. This implies that any initially created strongly inhomogeneous particle cloud is evolved into a thin pancake in a horizontal plane with very slow increase in its thickness in the vertical direction. The turbulent Schmidt number (the ratio of the eddy viscosity and the vertical turbulent diffusivity of the passive scalar) linearly increases with the gradient Richardson number. The physics of such a behavior is related to the buoyancy force that causes a correlation between fluctuations of the potential temperature and the particle number density. This correlation that is proportional to the product of the vertical turbulent particle flux and the vertical gradient of the mean potential temperature reduces the vertical turbulent particle flux. Considering the applications of these results to the atmospheric boundary-layer turbulence, the theoretical relationships are derived, which allows us to determine the turbulent diffusion tensor as a function of the vertical coordinate measured in the units of the local Obukhov length scale. The obtained relations are potentially useful in modeling applications of particle dispersion in the atmospheric boundary-layer turbulence and free atmosphere turbulence.