Abstract: There is an urgent need to reduce the computational costs and improve the convergence rate for the three-dimensional (3D) high-fidelity pin-by-pin full core simulation. Therefore, the efficient and robust acceleration method of 3D large-scale pin-by-pin neutron transport models is a primary objective of high-resolution reactor calculations. In this paper, comeSn_JFNK, an efficient unified parallel solver was developed for 3D steady-state and transient pin-by-pin neutron transport models. The comeSn_JFNK solver integrated the parallel discrete ordinate (S_N) neutron transport code comeSn into the parallel computational framework comeJFNK of Jacobian-Free Newton Krylov (JFNK). The comeSn code and comeJFNK framework were developed by the Virtual Reactor Coupling Analysis Team (VRCAT) at Huazhong University of Science and Technology (HUST). Therefore, comeSn_JFNK took advantages of the fast and robust convergence of the parallel JFNK framework and the high accuracy and efficiency of the S_N method based on the KBA algorithm. To further improve computational efficiency of comeSn_JFNK, in the JFNK solution, neutron scalar fluxes instead of neutron angular fluxes were chosen as the global solution variables. This can reduce the number of variables and minimize the computing scale. To rapidly construct the unified JFNK residuals, the parallel KBA transport sweep methods and physics-based preconditioning techniques were utilized to improve the computational efficiency. To rapidly construct the unified JFNK residuals, the parallel KBA transport sweep methods and physics-based preconditioning techniques were utilized. This way of constructing residuals, for both steady-state and transient models, can further improve the computational efficiency. Finally, the detailed analysis of the computational accuracy and acceleration characteristics of comeSn_JFNK was presented by solving the pin-by-pin steady-state KAIST-3A and homogenized-pin transient C5G7-TD2 benchmark cases. The relative errors of radial average power density, effective multiplication factor and relative powers as a function of time were shown in this paper. There are almost the same numerical solutions between the parallel comeSn code and the parallel solver called comSn_JFNK. Numerical results also show that the parallel solver called comeSn_JFNK can achieve an acceleration of over 10 times for KAIST-3A benchmark problems and approximately 30 times for C5G7-TD2 cases compared to the original parallel comeSn code using the traditional source iteration/power iteration methods. These solutions indicate that the JFNK methods offer significant acceleration for solving both steady-state and transient S_N neutron transport models compared to traditional source iteration/power iteration methods. In summary, the JFNK method provides reliable computational accuracy and high computational efficiency, which demonstrates the potential and advantages to accelerate neutron transport solutions. It also establishes a foundation for efficient simultaneous solution of more complicated pin-by-pin neutron transport and thermal-hydraulic coupling problems in the nuclear reactor cores using the unified JFNK framework in the coupling multiphysics environment (COME) developed by VRCAT at HUST.