We develop a general theory for describing phase coexistence between nonequilibrium steady states in Brownian systems, based on power functional theory [M. Schmidt and J. M. Brader, J. Chem. Phys. 138, 214101 (2013)]. We apply the framework to the special case of fluid-fluid phase separation of active soft sphere swimmers. The central object of the theory, the dissipated free power, is calculated via computer simulations and compared to a simple analytical approximation. The theory describes well the simulation data and predicts motility-induced phase separation due to avoidance of dissipative clusters. DOI: 10.1103/PhysRevLett.117.208003 Phase transitions in soft matter occur both in equilibrium and in nonequilibrium situations. Examples of the latter type include the glass transition [1], various types of shearbanding instabilities observed in colloidal suspensions [2,3], shear-induced demixing in semidilute polymeric solutions [4], and motility-induced phase separation in assemblies of active particles [5,6]. In contrast to phase transitions in equilibrium, which obey the statistical mechanics of Boltzmann and Gibbs, very little is known about general properties of transitions between out-of-equilibrium states. A corresponding universal framework for describing nonequilibrium soft matter is lacking at present.Theoretical progress has recently been made for the case of many-body systems governed by overdamped Brownian dynamics, encompassing a broad spectrum of physical systems [7]. It has been demonstrated that the dynamics of such systems can be described by a unique time-dependent power functional R t ½ρ; J, where the arguments are the space-and time-dependent one-body density distribution, ρðr; tÞ, and the one-body current distribution, Jðr; tÞ, in the case of a simple substance [8,9]. Both these fields are microscopically sharp and act as trial variables in a variational theory. The power functional theory is regarded to be "important, [as it] provides (i) a rigorous framework for formulating dynamical treatments within the [density functional theory] formalism and (ii) a systematic means of deriving new approximations" [10].The physical time evolution is that which minimizes R t ½ρ; J at time t with respect to Jðr; tÞ, while keeping ρðr; tÞ fixed. Hence,at the minimum of the functional. Here the variation is performed at fixed time t with respect to the positiondependent current. The density distribution is then obtained from integrating the continuity equation, ∂ρðr; tÞ=∂t ¼ −∇ · Jðr; tÞ, in time. The power functional possesses units of energy per time and can be split according towhere P t ½ρ; J accounts for the irreversible energy loss due to dissipation, _ F½ρ is the total time derivative of the intrinsic (Helmholtz) free energy density functional [7,11], and X t ½ρ; J is the external power, given bywhere _ V ext ðr; tÞ is the partial time derivative of the external potential V ext ðr; tÞ, and F ext ðr; tÞ is the external one-body force field, which in general consists of a sum of a conservative c...