The structural properties of neutron stars, namely, maximum mass is a subject that theoretical astrophysicists have desired to study. Theoretically, investigation of structure for neutron stars necessities obtaining the perfect hydrostatic equilibrium equations. Up to now, different models and metrics have been used to calculate these equations. It is expected higher order energy-momentum gravity plays crucial role in the investigation of neutron stars structure due to high density cores of neutron stars. Recently, the energy-momentum squared gravity (EMSG) is exposed .
In this paper, we will compute mass-radius relation for neutron stars using EMSG and employing realistic equations of state (EOS). We have obtained the EOS of neutron star matter using the lowest order constrained variational (LOCV) formalism as a microscopic approach [2, 3]. The LOCV method is a powerful tool for determination of the properties of the nucleonic matter at zero and finite temperatures [4-5].
 O. Akarsu, J. D. Barrow, S. Çıkıntoğlu, K. Yavuz Ekşi and N. Katirci,
Phys. Rev. D 97, 124017 (2018).
 Z. Asadi and M. Bigdeli, J. Phys. G: Nucl. Part. Phys. 45, 065101 (2018).
 Z. Asadi and M. Bigdeli, (2019) submitted to J. Phys. G: Nucl. Part. Phys (2019)
 M. Bigdeli, Phys. Rev . C 82, 054312 (2010).
 G. H. Bordbar and M. Modarres, Phys. Rev. C 57 (1998) 714.