Assuming $T_{0}$ to be an m-accretive operator in the complex Hilbert space ${\mathcal{H}}$, we use a resolvent method due to Kato to appropriately define the additive perturbation $T=T_{0}+W$ and prove stability of square root domains, that is, $$\begin{eqnarray}\text{dom}((T_{0}+W)^{1/2})=\text{dom}(T_{0}^{1/2}).\end{eqnarray}$$ Moreover, assuming in addition that $\text{dom}(T_{0}^{1/2})=\text{dom}((T_{0}^{\ast })^{1/2})$, we prove stability of square root domains in the form $$\begin{eqnarray}\text{dom}((T_{0}+W)^{1/2})=\text{dom}(T_{0}^{1/2})=\text{dom}((T_{0}^{\ast })^{1/2})=\text{dom}(((T_{0}+W)^{\ast })^{1/2}),\end{eqnarray}$$ which is most suitable for partial differential equation applications. We apply this approach to elliptic second-order partial differential operators of the form $$\begin{eqnarray}-\text{div}(a{\rm\nabla}\,\cdot )+(\vec{B}_{1}\cdot {\rm\nabla}\,\cdot )+\text{div}(\vec{B}_{2}\,\cdot )+V\end{eqnarray}$$ in $L^{2}({\rm\Omega})$ on certain open sets ${\rm\Omega}\subseteq \mathbb{R}^{n}$, $n\in \mathbb{N}$, with Dirichlet, Neumann, and mixed boundary conditions on $\partial {\rm\Omega}$, under general hypotheses on the (typically, non-smooth, unbounded) coefficients and on $\partial {\rm\Omega}$.