How do we make symmetry assumptions rigorous? I have, for instance, a problem with a spherically symmetric charge distribution. I deduce here, in order to solve the problem easily, that the corresponding electric field must be symmetric. How is this type of deduction rigorously justified?
 A: Let me first answer this question in a particular class of electrostatics problems.
In the case of a localized charge distribution in electrostatics (one in which the charge density vanishes outside of some ball around the origin), the general solution for a potential that vanishes at $\infty$ is
$$
  \mathbf \Phi(\mathbf x) = k\int_{\mathbb R^3} d^3x' \frac{\rho(\mathbf x')}{|\mathbf x - \mathbf x'|} 
$$
In this case, we often want to show that if $\rho$ has certain symmetries, then $\Phi$ shares those symmetries.  Let's first derive a little result:

Suppose that we make an invertible transformation $T$ on the charge distribution (say a translation for a rotation for example), then the new (transformed) charge distribution $\rho_T$ will be
$$
  \rho_T(\mathbf x) = \rho(T^{-1}\mathbf x).
$$
What will the resulting potential $\Phi$ be?  Well, let's compute
$$
  \Phi_T(\mathbf x) 
  = k\int_{\mathbb R^3} d^3 x' \frac{\rho_T(\mathbf x')}{|\mathbf x - \mathbf x'|}
  = k\int_{\mathbb R^3} d^3 x' \frac{\rho(T^{-1}(\mathbf x'))}{|\mathbf x - \mathbf x'|}
$$
we can perform the resulting integral via a change of variables
$$
  \mathbf u = T^{-1}(\mathbf x') \implies \mathbf x' = T(\mathbf u)
$$
and the formula for changing the variable of integration for volume integrals gives
$$
  d^3 x' = J_T(\mathbf u) d^3 u
$$
where $J_T$ is the jacobian of the transformation, so that the transformed potential becomes
$$
  \Phi_T(\mathbf x) = k\int_{\mathbb R^3} d^3u \,J_T(\mathbf u)\frac{\rho(\mathbf u)}{|\mathbf x - T(\mathbf u)|}
$$
Now back to the show.

To see how this formula for the transformed potential is used to answer your question about symmetries, let's consider a translationally symmetric charge density;
$$
  \rho(\mathbf x - \mathbf x_0) = \rho(\mathbf x), \qquad \text{for all $\mathbf x_0$}
$$
In this case, the transformation $T$ is $T(\mathbf x) = \mathbf x + \mathbf x_0$.  The Jacobian is just 1, and our the formula we derived above for the transformed potential gives
$$
  \Phi_T(\mathbf x) = k\int_{\mathbb R^3} d^3u \frac{\rho(\mathbf u)}{|\mathbf x - (\mathbf u-\mathbf x_0)|} = k\int_{\mathbb R^3} d^3u \frac{\rho(\mathbf u)}{|(\mathbf x - \mathbf x_0) - \mathbf u|} = \Phi(\mathbf x - \mathbf x_0)
$$
On the other hand, the translational invariance of the charge density tells us that
$$
  \Phi_T(\mathbf x) 
  = k\int_{\mathbb R^3} d^3 x' \frac{\rho_T(\mathbf x')}{|\mathbf x - \mathbf x'|}
  = k\int_{\mathbb R^3} d^3 x' \frac{\rho(\mathbf x')}{|\mathbf x - \mathbf x'|}
  = \Phi(\mathbf x)
$$
so combining these results gives
$$
  \Phi(\mathbf x - \mathbf x_0) = \Phi(\mathbf x)
$$
Namely, the potential is also translationally symmetric.  A similar procedure can be used for other symmetries.  Try rotational invariance for example on your own!
Hope that helps!
Physics Rocks.
Addendum. I think you can show similar things for generic Neumann or Dirichlet boundary value problems in electrostatics in which, for example, you don't just want to solve Poisson's equation for a localized charge distribution with vanishing potential at infinity, but I haven't worked out the details.
