Singularity-free stationary electro-vacuum solution Let's say we have a spherically symmetric fluid:
$$
T^{\alpha \beta} = \begin{bmatrix}
  \rho & 0 & 0 & 0 \\
     0 & p & 0 & 0 \\
     0 & 0 & p & 0 \\
     0 & 0 & 0 & p    
\end{bmatrix}
$$
where $\rho$ and $p$ are functions of the radial coordinate r.
If $p=-\rho/3$ this would be traceless, and possibly compatible with an electro-vacuum solution (Are there other necessary conditions for compatibility besides the stress-energy tensor being traceless?). While it would probably be unstable, is there a singularity free electrovacuum solution of GR which is stationary?
It seems that this has enough symmetry that someone skilled with differential-geometry and GR could actually probably just solve this themselves.  I'm interested in not just the answer, but also in seeing the process of how to solve this in GR.  So please include any work in your solution, thank you!
 A: You may be interested in the Bertotti-Robinson solution to the field equations. It is the only conformally flat solution to the electrovacuum equations with a non-null (in the sense of the electromagnetic invariants) Faraday tensor. (See also Chapter 7 of Griffiths and Podolsky, Exact space-times in Einstein's general relativity.) (Note that in the form presented in the ScienceWorld link it is not obviously static. See below.)
This solution, however, has a uniform electromagnetic field, and is not asymptotically flat. In the asymptotically flat case, using the Jebsen-Birkhoff-Eiseland rigidity theorem, a spherically symmetric solution (as you postulated in the first line of your question) must admit an additional symmetry, and the field equations reduce to ODEs. Asymptotic flatness allows you to prescribe asymptotic data "at infinity" and integrate back. In this case you have a stronger result that for any (possibly nonlinear) model of electromagnetism that reduces to Maxwell theory in the weak field limit, a spherically symmetric solution must be static and must also be singular along the symmetry axis. For a review of this fact and other related results, see the recent paper by Shadi Tavildar-Zadeh.  
Similar to the Bertotti-Robinson solution is also the Melvin Universe, which is a static non-singular, cylindrical symmetric space-time with an electric field directed along the axis of the cylinder. (You can find more about this solution in Griffiths and Podolsky's book mentioned above.) This solution isn't really physical: it was shown by Thorne, and by Melvin and Wallingford that the geometry is very strange in the Melvin Universe: time-like geodesics with only radial components to the velocity cannot escape to radial infinity ever; the geodesics oscillate and repeatedly cross the central symmetry axis. In otherwords, the gravitational effects induced by the electric field in this solution is very strong. (For a non-linear electrodynamics version of the Melvin Universe, there's a paper of Gibbons and Herdeiro on the subject.) 

Since you also asked how one coms up with these solutions: a good method is guess and check. Guess an ansatz (for example, a symmetry class), reduce the equation, and try to solve. For example, this would be how I derive the Bertotti-Robinson solution:
We will make the horribly unphysical ansatz that the solution is spherically symmetric and that the space-time is a direct product space-time. So in particular the geometry decomposes into $\mathbb{S}^2 \times (\mathbb{R}^2,g)$. Jebsen-Birkhoff-Eiseland still applies, so we can take one coordinate of $\mathbb{R}^2$ to be the Killing coordinate. 
The metric then takes the canonical form
$$ -A(r) dt^2 + B(r) dr^2 + ds^2_{\mathbb{S}^2}$$
As a direct product metric, the Ricci curvatures are just the sums of the component Ricci curvatures. So using the standard formulae for curvatures of two dimensional surfaces, we quickly compute that in the spherical directions the Ricci curvature is that of the standard sphere, which is equal to the metric (2 dimensional manifold with constant curvature 1). In the $\mathbb{R}^2$ direction, the Gaussian curvature is 
$$ - \frac{1}{2\sqrt{A(r)B(r)}} \left( \frac{A'(r)}{\sqrt{A(r)B(r)}}\right)' $$
The trace-free condition from Einstein-Maxwell requires that the Gaussian curvature equal to -1, which you can either solve explicitly for by integrating in $r$, or by "knowledge" that the only 2 dimensional Lorentzian manifold with constant negative curvature is AdS2. 
So if you take $r$ to be the "inverse radial" coordinate (where $r = 0$ corresponds to spatial infinity), the usual presentation of AdS2 gives that $A(r) = B(r) = r^{-2}$. 
For the Maxwell field, you observe that by spherical symmetry $F$ only has $dr\wedge dt$ component (if you toss out magnetic monopoles; you can gain back the other term using a duality rotation anyway). So $F$ is given by a scalar function $f$. The curvature condition implies that $R_{rr} = T_{rr} = F_{rt}F_{rt}g^{tt} = -f^2 A^{-1} = -1g_{rr} = -1/r^2$ (observe that you can say the same thing using $R_{tt}$ or the spherical components). This means that $f^2 = r^{-4}$, or $f = r^{-2}$. 
In remains to check that such a solution verifies Maxwell's equation. $dF = f' dr\wedge dr\wedge dt = 0$ clearly. while the Hodge dual $*F = \sqrt{|g|}g^{tt}g^{rr} f dvol_{\mathbb{S}^2} = - dvol_{\mathbb{S}^2}$, so $d*F = 0$ also. So the Bertotti-Robinson solution is indeed a static solution to the field equations. That the solution is non-singular follows from the fact that it is a direct product of two constant curvature spaces, so the curvature is bounded, and that the electromagnetic invariant $g(F,F)$ is constant valued, so the Maxwell field does not blowup either. 
A: This doesn't qualify as electrovac, but you can superpose arbitrarily placed extremal black holes, and making a soup of infinitely many infinitesimal charge/mass BH's, you get a static solution which is nontrivial, and only sourced by these infinitesimal sources which are nonsingular in the infinitesimal limit. The metric and electrostatic potential are determined in isotropic coordinates by a single solution to Laplaces equation in flat space. I forget the name of these, but it's in Hawking and Ellis
A: The answer to this question is likely no, because the Robertson-Walker spacetime with a radiation equation of state satisfies Einstein's equation for that stress energy tensor, and, famously, it only has a static solution in the case of the Einstein Static Universe, where the cosmological constant balances out the gravitational attraction of the expanding universe.  The staticity of the Einstein Static universe is known to be unstable under the influence small perturbations.  
The only reason I say 'likely no' and not no is that it may be possible to save a solution like this by relaxing the requirement that space be homogenous to just spherical symmetry, but I'm inclined to say that it's unlikely that this will end up working out.  
