The Debye series has been a key tool for the understanding of light scattering features, and it is also a convenient method for understanding and improving the design of optical instruments aimed at optical particle sizing. Gouesbet has derived the Debye series formulation for generalized Lorenz-Mie theory (GLMT). However, the scattering object is a homogeneous sphere, and no numerical result is provided. The Debye series formula for plane-wave scattering by multilayered spheres has been derived before. We have devoted our work to the Debye series of Gaussian beam scattering by multilayered spheres. The integral localized approximation is employed in the calculation of beam-shape coefficients (BSCs) and allows the study of the scattering characteristics of particles illuminated by the strongly focused beams. The formula and code are verified by the comparison with the results produced by GLMT and also by the comparison with the result for the case of plane-wave incidence. The formula is also employed in the simulation of the first rainbow by illuminating the particle with one or several narrow beams.