The aim of this paper is to show the existence of $\mathcal{R}$-bounded solution operator families for two-phase Stokes resolvent equations in $\dot\Omega =\Omega _+\cup\Omega _-$, where $\Omega _\pm$ are uniform $W_r^{2-1/r}$ domains of $N$-dimensional Euclidean space ${\mathbf{R}^N}$ ($N\geq 2$, $N < r < \infty$). More precisely, given a uniform $W_r^{2-1/r}$ domain $\Omega $ with two boundaries $ \Gamma _\pm$ satisfying $ \Gamma _+\cap \Gamma _-=\emptyset$, we suppose that some hypersurface $ \Gamma $ divides $\Omega $ into two sub-domains, that is, there exist domains $\Omega _\pm\subset\Omega $ such that $ \Omega _+\cap\Omega _-=\emptyset$ and $\Omega \setminus \Gamma =\Omega _+\cup\Omega _-, $ where $ \Gamma \cap \Gamma _+=\emptyset$, $ \Gamma \cap \Gamma _-=\emptyset$, and the boundaries of $\Omega _\pm$ consist of two parts $ \Gamma $ and $ \Gamma _\pm$, respectively. The domains $\Omega _\pm$ are filled with viscous, incompressible, and immiscible fluids with density $\rho_\pm$ and viscosity $\mu_\pm$, respectively. Here, $\rho_\pm$ are positive constants, while $\mu_\pm=\mu_\pm(x)$ are functions of $x\in{\mathbf{R}^N}$. On the boundaries $ \Gamma $, $ \Gamma _+$, and $ \Gamma _-$, we consider an interface condition, a free boundary condition, and the Dirichlet boundary condition, respectively. We also show, by using the $\mathcal{R}$-bounded solution operator families, some maximal $L_p{\theta}xt{-}L_q$ regularity as well as generation of analytic semigroup for a time-dependent problem associated with the two-phase Stokes resolvent equations. This kind of problems arises in the mathematical study of the motion of two viscous, incompressible, and immiscible fluids with free surfaces. The essential assumption of this paper is the unique solvability of a weak elliptic transmission problem for $\mathbf{f}\in L_q(\Omega )^N$, that is, it is assumed that the unique existence of solutions ${\theta}\in\mathcal{W}_q^1(\Omega )$ to the variational problem: $ (\rho^{-1}\nabla{\theta},\nabla{\varphi})_{\dot\Omega }=(\mathbf{f},\nabla{\varphi})_{\Omega } $ for any ${\varphi}\in\mathcal{W}_{q'}^1(\Omega )$ with $1 < q < \infty$ and $q'=q/(q-1)$, where $\rho$ is defined by $\rho=\rho_+$ ($x\in\Omega _+$), $\rho=\rho_-$ ($x\in\Omega _-$) and $\mathcal{W}_q^1(\Omega )$ is a suitable Banach space endowed with norm $\|\cdot\|_{\mathcal{W}_q^1(\Omega )}:=\|\nabla\cdot\|_{L_q(\Omega )}$. Our assumption covers e.g. the following domains as $\Omega $: ${\mathbf{R}^N}$, $\mathbf{R}_\pm^N$, perturbed $\mathbf{R}_\pm^N$, layers, perturbed layers, and bounded domains, where $\mathbf{R}_+^N$ and $\mathbf{R}_-^N$ are the open upper and lower half spaces, respectively.