In this work, our main concern is to study the existence and multiplicity of solutions for the following sub-elliptic system with Hardy type potentials and multiple critical exponents on Carnot group

$$\{\begin{array}{rlr}& -{\mathrm{\Delta }}_{\mathbb{G}}u=\frac{{\psi }^{\alpha }|u{|}^{2^{\ast }(\alpha )-2}u}{d(z{)}^{\alpha }}+\frac{p_1}{2^{\ast }(\gamma )}\frac{{\psi }^{\gamma }|u{|}^{p_1-2}u|v{|}^{p_2}}{d(z,z_0{)}^{\gamma }}+\lambda h(z)\frac{{\psi }^{\sigma }|u{|}^{q-2}u}{d(z{)}^{\sigma }}& \text{in}\mathrm{\Omega },\\ & -{\mathrm{\Delta }}_{\mathbb{G}}v=\frac{{\psi }^{\beta }|v{|}^{2^{\ast }(\beta )-2}v}{d(z{)}^{\beta }}+\frac{p_2}{2^{\ast }(\gamma )}\frac{{\psi }^{\gamma }|u{|}^{p_1}|v{|}^{p_2-2}v}{d(z,z_0{)}^{\gamma }}+\lambda h(z)\frac{{\psi }^{\sigma }|v{|}^{q-2}v}{d(z{)}^{\sigma }}& \text{in}\mathrm{\Omega },\\ & u=v=0& \text{on}\mathrm{\partial }\mathrm{\Omega },\end{array}$$

where $-{\mathrm{\Delta }}_{\mathbb{G}}$ is a sub-Laplacian on Carnot group $\mathbb{G}$, $\alpha ,\beta ,\gamma ,\sigma \in [0,2)$, $d$ is the ${\mathrm{\Delta }}_{\mathbb{G}}$-natural gauge, $\psi =|{\mathrm{\nabla }}_{\mathbb{G}}d|$ and ${\mathrm{\nabla }}_{\mathbb{G}}$ is the horizontal gradient associated to ${\mathrm{\Delta }}_{\mathbb{G}}$. The positive parameters $\lambda$, $q$ satisfy $0<\lambda <\mathrm{\infty }$, $1<q<2$, and $p_1$, $p_2>1$ with $p_1+p_2=2^{\ast }(\gamma )$, here $2^{\ast }(\alpha ):=\frac{2(Q-\alpha )}{Q-2}$, $2^{\ast }(\beta ):=\frac{2(Q-\beta )}{Q-2}$ and $2^{\ast }(\gamma )=\frac{2(Q-\gamma )}{Q-2}$ are the critical Hardy-Sobolev exponents, $Q$ is the homogeneous dimension of the space $\mathbb{G}$. By means of variational methods and the mountain-pass theorem of Ambrosetti and Rabonowitz, we study the existence of multiple solutions to the sub-elliptic system.

We are devoted to the study of the following sub-Laplacian system with Hardy-type potentials and critical nonlinearities

$$\{\begin{array}{r}-{\mathrm{\Delta }}_{\mathbb{G}}u-{\mu }_1\frac{{\psi }^{2}u}{\text{d}(z{)}^2}={\lambda }_1\frac{{\psi }^{\alpha }|u{|}^{2^{\ast }(\alpha )-2}u}{\text{d}(z{)}^{\alpha }}+\beta p_{1}f(z)\frac{{\psi }^{\gamma }|u{|}^{p_1-2}u|v{|}^{p_2}}{\text{d}(z{)}^{\gamma }}\text{in}\mathbb{G},\\ -{\mathrm{\Delta }}_{\mathbb{G}}v-{\mu }_2\frac{{\psi }^{2}v}{\text{d}(z{)}^2}={\lambda }_2\frac{{\psi }^{\alpha }|v{|}^{2^{\ast }(\alpha )-2}v}{\text{d}(z{)}^{\alpha }}+\beta p_{2}f(z)\frac{{\psi }^{\gamma }|u{|}^{p_1}|v{|}^{p_2-2}v}{\text{d}(z{)}^{\gamma }}\text{in}\mathbb{G},\end{array}$$

where $-{\mathrm{\Delta }}_{\mathbb{G}}$ is the sub-Laplacian on Carnot group $\mathbb{G}$, ${\mu }_1$, ${\mu }_2\in [0,{\mu }_{\mathbb{G}})$, $\alpha ,\gamma \in (0,2)$, ${\lambda }_1$, ${\lambda }_2$, $\beta$, $p_1$, $p_2>0$ with $1<p_1+p_2<2$, $\text{d}(z)$ is the ${\mathrm{\Delta }}_{\mathbb{G}}$-gauge, $\psi =|{\mathrm{\nabla }}_{\mathbb{G}}\text{d}(z)|$, $2^{\ast }(\alpha ):=\frac{2(Q-\alpha )}{Q-2}$ is the critical Sobolev-Hardy exponents, and ${\mu }_{\mathbb{G}}=(\frac{Q-2}{2}{)}^2$ is the best Hardy constant on $\mathbb{G}$. By combining a variant of the symmetric mountain pass theorem with the genus theory, we prove the existence of infinitely many weak solutions whose energy tends to zero when $\beta$ or ${\lambda }_1$, ${\lambda }_2$ belong to a suitable range.