  Among all the available observational techniques for studying magnetic fields
in the dense cold phase of the interstellar medium, linear polarization of
spectral lines, referred to in the literature as the Goldreich-Kylafis effect
(Goldreich & Kylafis 1981; hereafter "GK effect"), remains one of the most
underutilized methods. In this study, we implement the GK effect into the
multilevel, non-local thermodynamic equilibrium radiative transfer code PyRaTE.
Different modes of polarized radiation are treated individually with separate
optical depths computed for each polarization direction. We benchmark our
implementation against analytical results and provide tests for various
limiting cases. In agreement with previous theoretical results, we find that in
the multilevel case the amount of fractional polarization decreases when
compared to the two-level approximation, but this result is subject to the
relative importance between radiative and collisional processes. Finally, we
post-process an axially symmetric, non-ideal magnetohydrodynamic
chemo-dynamical simulation of a collapsing prestellar core and provide
theoretical predictions regarding the shape (as a function of velocity) of the
polarization fraction of CO during the early stages in the evolution of
molecular clouds. The code is freely available to download.
