  We introduce Monte Carlo-based non-LTE line radiative transfer calculations
in the 3D dust radiative transfer code SKIRT, which was originally set up as a
dust radiative transfer code. By doing so, we develop a generic and powerful 3D
radiative transfer code that can self-consistently generate spectra with
molecular and atomic lines against the underlying continuum. We test the
accuracy of the non-LTE line radiative transfer module in the extended SKIRT
code using standard benchmarks. We find excellent agreement between the SKIRT
results, the published benchmark results, and results obtained using the
ray-tracing non-LTE line radiative transfer code MAGRITTE, which validates our
implementation. We apply the extended SKIRT code on a 3D hydrodynamic
simulation of a dusty AGN torus model and generate multi-wavelength images with
CO rotational-line spectra against the underlying dust continuum. We find that
the low-J CO emission traces the geometrically thick molecular torus, whereas
the higher-J CO lines originate from the gas with high kinetic temperature
located in the innermost regions of the torus. Comparing the calculations with
and without dust radiative transfer, we find that higher-J CO lines are
slightly attenuated by the surrounding cold dust when seen edge-on. This shows
that atomic and molecular lines can experience attenuation, an effect that is
particularly important for transitions at mid- and near-infrared wavelengths.
Therefore, our self-consistent dust and non-LTE line radiative transfer
calculations can help interpret the observational data from Herschel, ALMA, and
JWST.
