The SPH method has become widely applicable due to accuracy improvements achieved through high-accuracy approximation models. Nevertheless, boundary treatment errors near the free surface still lead to accuracy degradation due to particle deficiency and the resulting instability of moment matrices. To address this issue, a regularization scheme has been proposed to stabilize computation. In this study, we employ a novel, physically consistent regularization scheme proposed by the authors. Furthermore, free surface boundary conditions can be imposed as constraints within the least squares fitting, enabling stable computation of high-accuracy differential operators even near the free surface. In addition to these numerical considerations, accurate representation of free surface behavior at the microscale requires appropriate modeling of surface tension. Various surface tension models have been proposed for this purpose. However, several existing approaches are sensitive to particle deficiency and density errors. As a result, their accuracy often relies on empirical corrections, which limits the robustness of curvature evaluation. Although high-accuracy surface tension models based on least squares fitting that evaluate curvature through local reconstruction of free surface geometry have been proposed, an SPH scheme that consistently integrates such models with numerical stability near the free surface and appropriate boundary condition treatments has not yet been fully established. In this study, we adopt least squares fitting for curvature evaluation. Moreover, by incorporating a physically consistent regularization scheme and the enhanced boundary treatment into the LS-SPH method, we propose a highly accurate and numerically stable framework for analyzing droplet dynamics.