We extend the range of validity of the ARTIS 3D radiative transfer code up to hundreds of days after explosion, when Type Ia supernovae are in their nebular phase. To achieve this, we add a non-local thermodynamic equilibrium (non-LTE) population and ionisation solver, a new multi-frequency radiation field model, and a new atomic dataset with forbidden transitions. We treat collisions with non-thermal leptons resulting from nuclear decays to account for their contribution to excitation, ionisation, and heating. We validate our method with a variety of tests including comparing our synthetic nebular spectra for the well-known one-dimensional W7 model with the results of other studies. As an illustrative application of the code, we present synthetic nebular spectra for the detonation of a sub-Chandrasekhar white dwarf in which the possible effects of gravitational settling of Ne22 prior to explosion have been explored. Specifically, we compare synthetic nebular spectra for a 1.06 M$_\odot$ white dwarf model obtained when 5.5 Gyr of very-efficient settling is assumed to a similar model without settling. We find that this degree of Ne22 settling has only a modest effect on the resulting nebular spectra due to increased Ni58 abundance. Due to the high ionisation in sub-Chandrasekhar models, the nebular [Ni II] emission remains negligible, while the [Ni III] line strengths are increased and the overall ionisation balance is slightly lowered in the model with Ne22 settling. In common with previous studies of sub-Chandrasekhar models at nebular epochs, these models overproduce [Fe III] emission relative to [Fe II] in comparison to observations of normal Type Ia supernovae.

Recent state-of-the-art calculations of A-values and electron impact excitation rates for Fe III are used in conjunction with the Cloudy modeling code to derive emission line intensity ratios for optical transitions among the fine-structure levels of the 3d$^6$ configuration. A comparison of these with high resolution, high signal-to-noise spectra of gaseous nebulae reveals that previous discrepancies found between theory and observation are not fully resolved by the latest atomic data. Blending is ruled out as a likely cause of the discrepancies, because temperature- and density-independent ratios (arising from lines with common upper levels) match well with those predicted by theory. For a typical nebular plasma with electron temperature $T_{\rm e} = 9000$ K and electron density $\rm N_{e}=10^4 \, cm^{-3}$, cascading of electrons from the levels $\rm ^3G_5$, $\rm ^3G_4$ and $\rm ^3G_3$ plays an important role in determining the populations of lower levels, such as $\rm ^3F_4$, which provide the density diagnostic emission lines of Fe III, such as $\rm ^5D_4$ - $\rm ^3F_4$ at 4658 \AA. Hence further work on the A-values for these transitions is recommended, ideally including measurements if possible. However, some Fe III ratios do provide reliable $N_{\rm e}$-diagnostics, such as 4986/4658. The Fe III cooling function calculated with Cloudy using the most recent atomic data is found to be significantly greater at $T_e$ $\simeq$ 30000 K than predicted with the existing Cloudy model. This is due to the presence of additional emission lines with the new data, particularly in the 1000--4000 \AA\ wavelength region.