  This paper deals with existence of solutions to the following fractional
$p$-Laplacian system of equations
  \begin{equation*}
  %\tag{$\mathcal P$}\label{MAT1}
  \begin{cases}
  (-\Delta_p)^s u =|u|^{p^*_s-2}u+
\frac{\gamma\alpha}{p_s^*}|u|^{\alpha-2}u|v|^{\beta}\;\;\text{in}\;\Omega,
  (-\Delta_p)^s v =|v|^{p^*_s-2}v+
\frac{\gamma\beta}{p_s^*}|v|^{\beta-2}v|u|^{\alpha}\;\;\text{in}\;\Omega,
  %
  % u,\;v\in\wsp,
  \end{cases}
  \end{equation*}
  where $s\in(0,1)$, $p\in(1,\infty)$ with $N>sp$, $\alpha,\,\beta>1$ such that
$\alpha+\beta = p^*_s:=\frac{Np}{N-sp}$ and $\Omega=\mathbb{R}^N$ or smooth
bounded domains in $\mathbb{R}^N$. For $\Omega=\mathbb{R}^N$ and $\gamma=1$, we
show that any ground state solution of the above system has the form $(\lambda
U, \tau\lambda V)$ for certain $\tau>0$ and $U,\;V$ are two positive ground
state solutions of $(-\Delta_p)^s u =|u|^{p^*_s-2}u$ in $\mathbb{R}^N$. For all
$\gamma>0$, we establish existence of a positive radial solution to the above
system in balls. For $\Omega=\mathbb{R}^N$, we also establish existence of
positive radial solutions to the above system in various ranges of $\gamma$.
