Abstract The nebular phase of lanthanide-rich ejecta of a neutron star merger (NSM) is studied by using a one-zone model, in which the atomic properties are represented by a single species, neodymium (Nd). Under the assumption that β-decay of r-process nuclei is the heat and ionization source, we solve the ionization and thermal balance of the ejecta under non-local thermodynamic equilibrium. The atomic data including energy levels, radiative transition rates, collision strengths, and recombination rate coefficients, are obtained by using atomic structure codes, GRASP2K and HULLAC. We find that both permitted and forbidden lines roughly equally contribute to the cooling rate of Nd II and Nd III at the nebular temperatures. We show that the kinetic temperature and ionization degree increase with time in the early stage of the nebular phase while these quantities become approximately independent of time after the thermalization break of the heating rate because the processes relevant to the ionization and thermalization balance are attributed to two-body collision between electrons and ions at later times. As a result, in spite of the rapid decline of the luminosity, the shape of the emergent spectrum does not change significantly with time after the break. We show that the emission-line nebular spectrum of the pure Nd ejecta consists of a broad structure from 0.5 μm to 20 μm with two distinct peaks around 1 μm and 10 μm.