The observation of gravitational waves from compact objects has now become an active part of observational astronomy. For a sound interpretation, one needs to compare such observations against detailed Numerical Relativity simulations, which are essential tools to explore the dynamics and physics of compact binary mergers. To date, essentially all simulation codes that solve the full set of Einstein's equations are performed in the framework of Eulerian hydrodynamics. The exception is our recently developed Numerical Relativity code SPHINCS_BSSN which solves the commonly used BSSN formulation of Einstein equations on a structured mesh and the matter equations via Lagrangian particles. We show here, for the first time, SPHINCS_BSSN neutron star merger simulations with piecewise polytropic approximations to four nuclear matter equations of state. We introduce some further methodological refinements (a new way of steering dissipation, an improved particlemesh mapping) and we explore the impact of the exponent that enters in the calculation of the thermal pressure contribution. We find that it leaves a noticeable imprint on the gravitational wave amplitude (calculated via both quadrupole approximation and the Ψ 4 -formalism) and has a noticeable impact on the amount of dynamic ejecta. Consistent with earlier findings, we only find a few times 10 −3 M as dynamic ejecta in the studied equal mass binary systems, with softer equations of state (which are more prone to shock formation) ejecting larger amounts of matter. We also see a credible high-velocity (∼ 0.5..0.7c) ejecta component of ∼ 10 −4 M in all our cases. Such a high-velocity component has been suggested to produce an early, blue precursor to the main kilonova emission and it could also potentially cause a kilonova afterglow.