  We study the weakly coupled critical elliptic system \begin{equation*}
\begin{cases} -\Delta u=\mu_{1}|u|^{2^{*}-2}u+\lambda\alpha
|u|^{\alpha-2}|v|^{\beta}u & \text{in }\Omega,\\ -\Delta
v=\mu_{2}|v|^{2^{*}-2}v+\lambda\beta |u|^{\alpha}|v|^{\beta-2}v & \text{in
}\Omega,\\ u=v=0 & \text{on }\partial\Omega, \end{cases} \end{equation*} where
$\Omega$ is a bounded smooth domain in $\mathbb{R}^{N}$, $N\geq 3$,
$2^{*}:=\frac{2N}{N-2}$ is the critical Sobolev exponent, $\mu_{1},\mu_{2}>0$,
$\alpha, \beta>1$, $\alpha+\beta =2^{*}$ and $\lambda\in\mathbb{R}$.
  We establish the existence of a prescribed number of fully nontrivial
solutions to this system under suitable symmetry assumptions on $\Omega$, which
allow domains with finite symmetries, and we show that the positive least
energy symmetric solution exhibits phase separation as $\lambda\to -\infty$.
  We also obtain existence of infinitely many solutions to this system in
$\Omega=\mathbb{R}^N$.
