In this paper, systems of fractional Laplacian equations are investigated, which involve critical homogeneous nonlinearities and Hardy-type terms as follows $$\begin{aligned} {\left\{ \begin{array}{ll} &{}(- \Delta )^s u - \gamma \frac{u}{|x|^{2 s}} = \frac{\lambda \mu }{\mu + \eta } \frac{|u|^{\mu - 2} |v|^{\eta } u}{|x|^{\alpha }}+ \frac{Q(x)}{2^*_s (\alpha )}\frac{ H_u(u,v)}{|x|^{\alpha }}, ~ in~ {\mathbb {R}}^N,\\ \\ &{} (- \Delta )^s v - \gamma \frac{v}{|x|^{2 s}} = \frac{\lambda \eta }{\mu + \eta } \frac{|u|^{\mu } |v|^{\eta - 2} v}{|x|^{\alpha }} + \frac{Q(x)}{2^*_s (\alpha )}\frac{ H_v(u,v)}{|x|^{\alpha }}, ~ in ~{\mathbb {R}}^N, \end{array}\right. } \end{aligned}$$ where $$0< s < 1$$ , $$0< \alpha< 2s <N$$ , $$0< \gamma < \gamma _H$$ with $$\begin{aligned} \gamma _H =4^s \frac{\Gamma ^2(\frac{N+2s}{4})}{\Gamma ^2(\frac{N-2s}{4})}, \;\;\; \text {for} \; i = 1,\ldots ,l \end{aligned}$$ being the fractional best Hardy constant on $${\mathbb {R}}^N$$ , $$2^*_s (\alpha ) = \frac{2 (N -\alpha )}{N - 2 s}$$ is critical Hardy–Sobolev exponent, $$\mu , \eta > 1$$ with $$\mu + \eta = 2^*_s(\alpha )$$ , Q is G-symmetric functions (G is a closed subgroup of O(N), see Sect. 2 for details) satisfying some appropriate conditions which will be specified later, $$\lambda $$ is real parameter, $$H_u$$ , $$H_v$$ are the partial derivatives of the 2-variable $$C^1$$ -functions H(u, v) and $$(-\Delta )^s$$ is the fractional Laplacian operator which (up to normalization factors) may be defined as $$\begin{aligned} (- \Delta )^s u (x) = - \frac{1}{2}\int _ {{\mathbb {R}}^{N} } \frac{u (x + y) + u (x - y) - 2 u (x) }{|x-y|^{N + 2 s}}{\text {d}}y. \end{aligned}$$ By variational methods and local compactness of Palais–Smale sequences, the extremals of the corresponding best Hardy–Sobolev constant are found and the existence of solutions to the system is established.