  The equations of state of the neutron star at the hadron-quark crossover
region are interpolated with the Gaussian process regression (GPR) method,
which can reduce the randomness of present interpolation schemes. The
relativistic mean-field (RMF) model and Nambu-Jona-Lasinio (NJL) model are
employed to describe the hadronic phase and quark phase, respectively. In the
RMF model, the coupling term between $\omega$ and $\rho$ mesons is considered
to control the density-dependent behaviors of symmetry energy, i.e. the slope
of symmetry energy, $L$. Furthermore, the vector interaction between quarks is
included in the NJL model to obtain the additional repulsive contributions.
Their coupling strengths and the crossover windows are discussed in the present
framework under the constraints on the neutron star from gravitational wave
detections, massive neutron star measurements, mass-radius simultaneous
observation of NICER collaboration, and the neutron skin thickness of
$^{208}$Pb from PREX-II. It is found that the slope of symmetry energy, $L$
should be around $50-90$ MeV and the crossover window is $(0.3,~0.6)~\rm
fm^{-3}$ with these observables. Furthermore, the uncertainties of neutron star
masses and radii in the hadron-quark crossover regions are also predicted by
the GPR method.
