  We investigate the nonrotating neutron stars in $f(T)$ gravity with
$f(T)=T+\alpha{T}^2$, where $T$ is the torsion scalar in the teleparallel
formalism of gravity. In particular, we utilize the SLy and BSk family of
equations of state for perfect fluid to describe the neutron stellar matter and
search for the effects of the $f(T)$ modification on the models of neutron
stars. For positive $\alpha$, the modification results in a smaller stellar
mass in comparison to general relativity, while the neutron stars will contain
larger amount of matter for negative $\alpha$. Moreover, there seems to be an
upper limit for the central density of the neutron stars with $\alpha>0$,
beyond which the effective $f(T)$ fluid would have a steplike phase transition
in density and pressure profiles, collapsing the numerical system. We obtain
the mass-radius relations of the realistic models of neutron stars and subject
them to the joint constraints from the observed massive pulsars PSR J0030+0451,
PSR J0740+6620, and PSR J2215+5135, and gravitational wave events GW170817 and
GW190814. For the neutron star model in $f(T)$ gravity to be able to
accommodate all the mentioned data, the model parameter $\alpha$ needs to be
smaller than $-4.295$, $-6.476$, $-4.4$, and $-2.12$ (in the unit of
${G}^2M_\odot^2/c^4$) for SLy, BSk19, BSk20, and BSk21 equations of state,
respectively. If one considers the unknown compact object in the event GW190814
not to be a neutron star and hence excludes this dataset, the constraints can
be loosened to $\alpha<-0.594$, $-3.5$, $0.4$ and $1.9$ (in the unit of
${G}^2M_\odot^2/c^4$), respectively.
