The hydrodynamic interaction between two particles suspended in shear flows is fundamental to the macroscopic characterization of suspension flows. Although such interaction in quiescent or linear shear flow is well understood, studies on that in a nonlinear shear field are rare. In this study, the hydrodynamic interaction between two neutrally-buoyant smooth spheres moving at negligible Reynolds numbers in an unbounded plane Poiseuille flow has been calculated by three-dimensional boundary element method (BEM) simulations. The BEM results have been compared with the analytical results obtained with the method-of-reflection (MoR) approximations. The BEM simulations have been found to provide satisfactory predictions if the number of elements on the spheres are more than 200, whereas the MoR approximations provide satisfactory predictions only when the minimum separation between the spheres is relatively large although this MoR method has the advantage to easily calculate the hydrodynamic interaction between two spheres freely moving at negligible Reynolds numbers in unbounded quadratic flow by solving ordinary differential equations. Furthermore, it is found that there is a preferential cross-streamline migration of the center-of-gravity of the sphere-pair in the plane of shear in plane Poiseuille flow which does not arise in simple shear flow. This migration is always directed towards low shear regions when the sphere having larger translational velocity approaches the other sphere, and reverses towards high shear regions when the faster sphere leads the other sphere in the plane of shear. There is also a cross-streamline migration of the center-of-gravity of the sphere-pair in the plane of vorticity, but this migration does not have a preferential direction. These migrations are symmetric about the point where the spheres are at the minimum separation, and are only significant when the hydrodynamic interaction of the spheres is strong. These results show that the migration of the center-of-gravity of the sphere-pair can be attributed to the nonlinearity of the shear field, which agrees with the MoR approximations. The hydrodynamic interaction between the two spheres has been quantified under various conditions by the BEM simulations for both identical and disparate spheres.