The dynamics of the neutral atoms in Hall thrusters affects several plasma processes, from ionization to electrons' mobility. In the context of Hall thruster's particle-in-cell (PIC) modeling, the neutrals are often treated kinetically, similar to the plasma species, and their interactions with themselves and the ions are resolved using the direct-simulation Monte–Carlo (DSMC) algorithm. However, the DSMC approach is computationally resource demanding. Therefore, modeling the neutrals as a 1D fluid has been also pursued in simulations that do not involve the radial coordinate and, hence, do not resolve the neutrals' radial expansion. In this article, we present an extensive study on the sensitivity of the PIC simulations of Hall thruster discharge to the model used for the neutral dynamics. We carried out 1D axial PIC simulations with various fluid and kinetic models of the neutrals as well as self-consistent quasi-2D axial-azimuthal simulations with different neutrals’ fluid descriptions. Our results show that the predictions of the simulations in either 1D or 2D configurations are highly sensitive to the neutrals' model, and that different treatments of the neutrals change the spatiotemporal evolution of the discharge. Moreover, we observed that considering the ion-neutral collisions causes a significant variation in the neutral temperature, thus requiring that the neutrals' energy equation to be included as well in their fluid system of equations. Finally, the self-consistent axial-azimuthal simulations highlighted that a neutrals’ model based on the continuity conservation equation only is not an appropriate choice and leads to physically unexpected high-frequency global discharge oscillations.