In the study of optical properties of large atomic system, a weak laser
driving is often assumed to simplify the system dynamics by linearly coupled
equations. Here we investigate the light scattering properties of atomic
ensembles beyond weak-field excitation through cumulant expansion method. By
progressively incorporating higher-order correlations into the steady-state
equations, an enhanced accuracy can be achieved in comparison to the exact
solutions from solving a full density matrix. Our analysis reveals that, in the
regime of weak dipole-dipole interaction (DDI), the first-order expansion
yields satisfactory predictions for optical depth, while denser atomic
configurations necessitate consideration of higher-order correlations. As the
intensity of incident light increases, atom saturation effects become
noticeable, giving rise to significant changes of light transparency, energy
shift, and decay rate. This saturation phenomenon extends to subradiant atom
arrays even under weak driving conditions, leading to substantial deviations
from the linear model. Our findings demonstrate the potential of mean-field
models as good extensions to linear models as it balances both accuracy and
computational complexity, which can be an effective tool for probing optical
properties in large atom systems. However, the crucial role of higher-order
cumulants in large atom systems under finite laser field excitations remains
unclear since it is challenging theoretically owing to the
exponentially-increasing Hilbert space in such light-matter interacting
systems.Comment: 4 figure