We consider M/Ph/n + M queueing systems in steady state. We prove that the Wasserstein distance between the stationary distribution of the normalized system size process and that of a piecewise Ornstein Uhlenbeck (OU) process is bounded by C/root T., where the constant C is independent of the arrival rate A and the number of servers n as long as they are in the HalfinWhitt parameter regime. For each integer m > 0, we also establish a similar bound for the difference of the mth steady-state moments. For the proofs, we develop a modular framework that is based on Stein's method. The framework has three components: Poisson equation, generator coupling, and state space collapse. The framework, with further refinement, is likely applicable to steady-state diffusion approximations for other stochastic systems.