Abstract The thermal structure of the chromosphere is regulated through a complex interaction of various heating processes, radiative cooling, and the ionization degree of the plasma. Here, we study the impact on the thermal properties of the chromosphere when including the combined action of nonequilibrium ionization (NEI) of hydrogen and helium and ion–neutral interaction effects. We have performed a 2.5D radiative magnetohydrodynamic simulation using the Bifrost code. This model includes ion–neutral interaction effects by solving the generalized Ohm’ s law (GOL) as well as NEI for hydrogen and helium. The GOL equation includes ambipolar diffusion and the Hall term. We compare this simulation with another simulation that computes the ionization in local thermodynamic equilibrium (LTE) including ion–neutral interaction effects. Our numerical models reveal substantial thermal differences in magneto-acoustic shocks, the wake behind the shocks, spicules, low-lying magnetic loops, and the transition region. In particular, we find that heating through ambipolar diffusion in shock wakes is substantially less efficient, while in the shock fronts themselves it is more efficient, under NEI conditions than when assuming LTE.