1. Introduction
Recently, the following Hartree-type system has attracted a lot of attention:
\begin{equation}
\left\{
\begin{aligned}
&-\Delta u=\left(\frac{1}{|x|^{n-\alpha}}\ast v^p\right)v^{p-1},\quad u \gt 0\ \text{in} \ \mathbb{R}^{n},\\
&-\Delta v=\left(\frac{1}{|x|^{n-\beta}}\ast u^q\right)u^{q-1},\quad v \gt 0\ \text{in} \ \mathbb{R}^{n},
\end{aligned}
\right.
\end{equation}where 
$n\geq3,$ 
$0 \lt \alpha, \beta \lt n,$ and 
$\min\{p,q\} \gt 1.$ This system can be viewed as a generalization of the static Hartree equation:
\begin{equation}
-\Delta u=\left(\frac{1}{|x|^{n-\alpha}}\ast u^p\right)u^{p-1},\quad u \gt 0\ \text{in} \ \mathbb{R}^{n},
\end{equation}which was studied extensively. Such an equation has many applications in the Hartree–Fock theory of the non-linear Schrödinger equations and the quantum theory of large systems of non-relativistic bosonic atoms and molecules (cf. [Reference Fröhlich and Lenzmann11] and [Reference Lieb and Simon30] and many others). It is also helpful in understanding the blowing up or the global existence and scattering of the solutions of the dynamic Hartree equation (cf. [Reference Li, Miao and Zhang28]).
When α = 2, Equation (1.2) has no positive solution if 
$1\leq p \lt \frac{n+2}{n-2},$ and all positive solutions are classified to the unique form 
$u(x)=c(\frac{t}{t^2+|x-x^*|^2})^{\frac{n-2}{2}}$ in the critical case 
$p=\frac{n+2}{n-2}$ (cf. [Reference Lei27]). Furthermore, the author also pointed out that the equation has positive stable solutions if and only if 
$p\geq 1+\frac{4}{n-4-2\sqrt{n-1}}.$ Afterwards, the same results for Equation (1.2) were obtained in [Reference Le23], and the author covered the full range for 
$0 \lt \alpha \lt n$ and 
$-\infty \lt p \lt \frac{n+\alpha}{n-2}.$ In addition, by the method of moving planes in integral forms, Du–Yang [Reference Du and Yang10] and Guo–Hu–Peng–Shuai [Reference Guo, Hu, Peng and Shuai17] gave the symmetry and uniqueness of the positive solutions of Equation (1.2) with the Sobolev-type critical exponent 
$p=2_\alpha^*:=\frac{n+\alpha}{n-2}.$ The existence of the super-solutions of Equation (1.2) and several sufficient conditions were studied in [Reference Moroz and Van Schaftingen36]. When 
$\alpha=n-4,$ p = 2, Liu [Reference Liu31] classified all 
$L^{\frac{2n}{n-2}}(\mathbb{R}^{n})$ solutions for the equivalent integral system of Equation (1.2). Afterwards, the integrability for all 
$L^{\frac{2n}{2+\alpha}}(\mathbb{R}^{n})$ solutions of Equation (1.2) with p = 2 was obtained, and decay rates of those solutions at infinity were estimated (see [Reference Lei26]). In addition, Equation (1.2) in the fractional setting was also studied (cf. [Reference Dai, Fang and Qin8, Reference Dai, Huang, Qin, Wang and Fang9, Reference Le21, Reference Ma, Shang and Zhang32]). Other results can refer to [Reference Le20, Reference Le22, Reference Yang and Yu43] and the references therein. Recently, Ghergu et al. [Reference Ghergu, Liu, Miyamoto and Moroz12] shows a necessary and sufficient condition of existence of super-solutions of
\begin{equation}
u(x)=|x|^{\alpha-n}*[u^{p-1}(|x|^{\beta-n}*u^p)], \quad u \gt 0 \ \ {\rm on} \ \ \mathbb{R}^n.
\end{equation}As a corollary of this result, we can obtain that Equation (1.2) has positive distributional super-solutions if
\begin{equation}p \gt \operatorname{max}\left\{2,\frac{2n-2+\alpha}{2(n-\alpha)},\frac{n+\alpha-2}{n-\alpha},\frac{2n}{2n-\alpha-2}\right\}.\end{equation} There are few researches on Equation (1.1) unlike on Equation (1.2). Recently, Wang and Yang [Reference Wang and Yang42] proved that u and v must be radially symmetric if 
$(u, v)\in L^{\frac{2n}{n-2}}(\mathbb{R}^{n})\times L^{\frac{2n}{n-2}}(\mathbb{R}^{n})$ is a positive solution of Equation (1.1) with 
$\alpha=\beta\geq n-4$ and 
$p=q=\frac{n+\alpha}{n-2}.$ But they did not give the explicit form of the solutions and restricted the global integrability assumption. In 2021, Le [Reference Le24] showed that the system (1.1) has no positive solution if
\begin{equation}
1 \lt p\leq \frac{n+\alpha}{n-2},\quad
1 \lt q\leq \frac{n+\beta}{n-2}
\quad \textrm{and} \quad (p,q)\neq\left(\frac{n+\alpha}{n-2},\frac{n+\beta}{n-2}\right)
\end{equation}by the method of moving spheres in integral forms. He also classified all positive solutions in the critical case 
$(p,q)=(\frac{n+\alpha}{n-2},\frac{n+\beta}{n-2}).$
In this paper, we investigate the nonexistence of positive solutions of the Hartree–Poisson system:
\begin{equation}
\left\{
\begin{aligned}
&-\Delta u=\left(\frac{1}{|x|^{n-2}}\times v^p\right)v^{p-1},\quad u \gt 0\ \text{in} \ \mathbb{R}^{n},\\
&-\Delta v=\left(\frac{1}{|x|^{n-2}}\times u^q\right)u^{q-1},\quad v \gt 0\ \text{in} \ \mathbb{R}^{n},
\end{aligned}
\right.
\end{equation}where 
$n \geq 3$ and 
$\min\{p,q\} \gt 1$. Actually, it seems difficult to investigate directly the properties of Equation (1.6) in view of the convolution term. Noting the relation between the Newton potential and the convolution properties of Dirac function, we can see that Equation (1.6) can be studied by the following Hartree–Poisson system:
\begin{equation}
\left\{
\begin{aligned}
&-\Delta u=wv^{p-1},\quad &u \gt 0\ \text{in} \ \mathbb{R}^{n}, \\
& -\Delta w=v^p,\quad &w \gt 0\ \text{in} \ \mathbb{R}^{n},\\
&-\Delta v=zu^{q-1},\quad &v \gt 0\ \text{in} \ \mathbb{R}^{n},\\
&-\Delta z=u^q,\quad &z \gt 0\ \text{in} \ \mathbb{R}^{n},
\end{aligned}
\right.
\end{equation}and the integral system
\begin{equation}
\left\{
\begin{aligned}
&u(x)=C_1\int_{\mathbb{R}^n}\frac{v^{p-1}(y)w(y)}{|x-y|^{n-2}}\,{\rm d}y,\\
&w(x)=\int_{\mathbb{R}^n}\frac{v^p(y)}{|x-y|^{n-2}}\,{\rm d}y,\\
&v(x)=C_2 \int_{\mathbb{R}^n}\frac{u^{q-1}(y)z(y)}{|x-y|^{n-2}}\,{\rm d}y,\\
&z(x)=\int_{\mathbb{R}^n}\frac{u^q(y)}{|x-y|^{n-2}}\,{\rm d}y.
\end{aligned}
\right.
\end{equation}Here, C 1 and C 2 are positive constants.
System (1.7) is related to the Lane–Emden system:
\begin{equation}
\left\{
\begin{aligned}
&-\Delta u=v^p,\quad u, v \gt 0 \; \text{in} \ \mathbb{R}^{n}, \\
&-\Delta v=u^q,\quad p,q \gt 0,
\end{aligned}
\right.
\end{equation}which arises in chemical, biological and physical sciences. One of the most concerned issues with Equation (1.9) is the Lane–Emden conjecture, which is still open. That is, the system (1.9) admits no positive classical solutions in the Sobolev-type subcritical case:
\begin{equation}
\frac{1}{p+1}+\frac{1}{q+1} \gt \frac{n-2}{n}.
\end{equation} Fortunately, partial results have been obtained. For 
$n\leq2,$ the conjecture is a consequence of a relatively easier, and the known result (1.9) has no positive super-solution if pq < 1 or pq > 1 and 
$\max\left\{\frac{2(p+1)}{pq-1}, \frac{2(q+1)}{pq-1}\right\}\geq n-2$ (which is called the Serrin-type condition). This result can be found in [Reference Mitidieri35] or [Reference Serrin and Zou39]. Therefore, an interesting open case is 
$n\geq3.$ In 1996, Mitidieri [Reference Mitidieri35] proved that the conjecture is true for radial solutions in all dimensions. Combining with this result, Chen–Li [Reference Chen and Li4] settled the conjecture under some integrability condition. When n = 3, it was proved in the full range (1.10) but under the additional assumption that u and v have polynomial growth at infinity (cf. [Reference Serrin and Zou39]). Afterwards, Polacik–Quittner–Souplet [Reference Polacik, Quittner and Souplet37] removed this assumption and proved the Lane–Emden conjecture when n = 3. And, when n = 4, the conjecture was solved by Souplet [Reference Souplet40]. When 
$n\geq5,$ the non-existence of positive classical solutions to Equation (1.9) is still unknown.
Inspired by these results, we will study the Lane–Emden conjecture for Equation (1.6). We say that Equation (1.6) is in the critical case when the pair (p, q) satisfies
\begin{equation}
\frac{1}{p}+\frac{1}{q}=\frac{2(n-2)}{n+2},
\end{equation}which is the Sobolev hyperbola. Throughout the paper, the case where the relation
\begin{equation*}
\frac{1}{p}+\frac{1}{q} \lt \frac{2(n-2)}{n+2},
\end{equation*}holds is referred to as the supercritical case, and the case
\begin{equation}
\frac{1}{p}+\frac{1}{q} \gt \frac{2(n-2)}{n+2},
\end{equation}holds is referred to as the subcritical case.
We always assume in this paper that 
$u,v \in L_{\rm loc}^1(\mathbb{R}^n)$ have the following slowly increasing properties, which implies the convolutions in Equation (1.6) make sense:
\begin{equation}
v^p, v^{p-1}u, u^q, u^{q-1}v \in L^1((1+|x|)^{2-n}dx, \mathbb{R}^n).
\end{equation} We say that (u, v) is a pair of positive distributional solution of Equation (1.6), if positive functions 
$u,v \in L_{\rm loc}^1(\mathbb{R}^n)$ satisfy Equation (1.13), and for any function 
$\phi \in C_0^\infty(\mathbb{R}^n)$, there holds
\begin{equation*}\begin{cases}
\displaystyle\int_{\mathbb{R}^n}u(x)[-\Delta \phi(x)]\,{\rm d}x=
\int_{\mathbb{R}^n}\phi(x) v^{p-1}(x)
\int_{\mathbb{R}^n}\frac{v^p(y)\,{\rm d}y}{|x-y|^{n-2}}\,{\rm d}x;\\
\displaystyle\int_{\mathbb{R}^n}v(x)[-\Delta \phi(x)]\,{\rm d}x=
\int_{\mathbb{R}^n}\phi(x) u^{q-1}(x)
\int_{\mathbb{R}^n}\frac{u^q(y)\,{\rm d}y}{|x-y|^{n-2}}\,{\rm d}x.
\end{cases}
\end{equation*} Furthermore, (u, v) is called a pair of positive classical solution of Equation (1.6), if positive functions 
$u,v \in C^2(\mathbb{R}^n)$ satisfy Equations (1.13) and (1.6) pointwise in 
$\mathbb{R}^n$.
Now, we give the main results.
First, the Serrin-type condition is a necessary condition of existence of positive solutions.
Theorem 1.1. Let 
$\min\{p,q\} \gt 1$, and
\begin{equation*}
\max\left\{\frac{8p}{(2p-1)(2q-1)-1}, \frac{8q}{(2p-1)(2q-1)-1}\right\}\geq n-2.
\end{equation*}Then, Equation (1.6) has no positive distributional solution.
Next, the Sobolev-type condition is a weaker one for existence of positive radial solutions.
Theorem 1.2. In the subcritical case (1.12), Equation (1.6) has no positive radial classical solutions.
Remark 1.1. (i) The Lane–Emden conjecture states that the Lane–Emden system has no positive classical solution under the subcritical condition. Mitidieri confirmed this conjecture with the radial classical solution in [Reference Mitidieri35]. Therefore, we only consider the non-existence of positive radial classical solutions here.
(ii) A natural question is when a distributed solution of Hartree–type equations is classical. To the best of our knowledge, there is no conclusion regarding Equation (1.6). For Equation (1.2), there are several results. First, Le [Reference Le25] pointed out that the solution 
$u \in L_{\rm loc}^{2n(p-1)/(\alpha+\beta)}(\mathbb{R}^n)$ of integral Equation (1.3) can be classified when p is the critical exponent, and hence u is classical. Similar to the proof of Theorem 1.1 in [Reference Dai, Huang, Qin, Wang and Fang9], we know that the distributional solutions of Equation (1.2) satisfy the integral equation when 
$\alpha \in (0,n/2)$. In addition, [Reference Dai, Huang, Qin, Wang and Fang9] shows that if 
$u \in H^{\alpha/2}(\mathbb{R}^n) \cap L^{2n/(n-\alpha)}(\mathbb{R}^n)$ is a distributional solution of 
$(-\Delta)^{\alpha/2}=u(|x|^{-2\alpha}*|u|^2)$, then 
$u \in C^{[\alpha]}(\mathbb{R}^n)$ as long as 
$\alpha \in [1,n/3)$. Thus, when 
$\alpha \geq 2$, those solutions are classical, and hence many elliptic methods (such as the Schauder estimation, the strong maximum principle, the method of moving planes, etc.) still work.
(iii) Although Equation (1.9) has no positive radial classical solution in the subcritical case, Equation (1.9) has other positive radial solutions, which are not continuous. For simplicity, we consider the case of 
$u \equiv v$. Now, both Equations (A.8) and (A.11) in [Reference Gidas and Spruck15] show that the Lane–Emden equation has radial solutions, which do not belong to 
$C^2(\mathbb{R}^n)$. We believe that analogous conclusions still hold true for Equation (1.6).
To proving Theorem 1.2, we apply the ideas in [Reference Mitidieri35]. To deal with convolution terms in Equation (1.6), we introduce two new unknown functions w and z, so Equation (1.6) is replaced by Equation (1.7), including four equations. Therefore, the process becomes more complicated when we combine those four equations into the Pohozaev equation.
In addition, the Sobolev-type condition is also a weaker one for existence of positive integrable solutions (i.e., 
$(u,v) \in
L^{r_0}(\mathbb{R}^n) \times L^{s_0}(\mathbb{R}^n)$), where
\begin{equation}
r_0=\frac{n[(2p-1))(2q-1)-1]}{8p},\qquad s_0=\frac{n[(2p-1))(2q-1)-1]}{8q}.
\end{equation}Theorem 1.3. In the subcritical case (1.12), Equation (1.6) has no positive classical solution in 
$L^{r_0}(\mathbb{R}^n) \times L^{s_0}(\mathbb{R}^n)$.
Remark 1.2. The Coulomb–Sobolev space 
$E^{2,s}(\mathbb{R}^n)$ is the vector space of functions 
$u \in \mathcal{D}^{1,2}(\mathbb{R}^n)$ such that (cf. Definition 2 in [Reference Mercuri, Moroz and Van Schaftingen33])
\begin{equation*}
\int_{\mathbb{R}^n}|\nabla u|^2\,{\rm d}x+
\int_{\mathbb{R}^n}||x|^{1-n}\times|u|^s|^2 \,{\rm d}x \lt \infty.
\end{equation*} If 
$(u,v) \in E^{2,s_1}(\mathbb{R}^n) \times E^{2,s_2}(\mathbb{R}^n)$, by Proposition 3.1 in [Reference Mercuri, Moroz and Van Schaftingen33], we have
\begin{equation*}
(u,v) \in L^{t_1}(\mathbb{R}^n) \times L^{t_2}(\mathbb{R}^n), \quad \textrm{for \ all} \
\ (t_1,t_2) \in I_1 \times I_2,
\end{equation*}where Ii (
$i=1,2$) are the closed intervals with endpoints 
$s_i+1$ and 
$2n/(n-2)$. When s 1 and s 2 satisfy 
$(r_0,s_0) \in I_1 \times I_2$, by Theorem 1.3, we know that Equation (1.6) has no positive classical solution in 
$E^{2,s_1}(\mathbb{R}^n) \times E^{2,s_2}(\mathbb{R}^n)$ under the subcritical condition (1.12).
To prove Theorem 1.3, we first convert Equation (1.6) to Equation (1.8). Next, using the method of the moving planes in integral form, we prove that all integrable solutions (i.e., 
$(u,v) \in L^{r_0}(\mathbb{R}^n)
\times L^{s_0}(\mathbb{R}^n)$ are radially symmetric. Thus, using Theorem 1.2, we see that Equation (1.6) has no classical solution in 
$L^{r_0}(\mathbb{R}^n) \times L^{s_0}(\mathbb{R}^n)$ under the subcritical condition (1.12).
Finally, we establish the equivalence between the Liouville theorem of Equation (1.6) and the estimate of boundary blowing-up rate for positive classical solutions of
\begin{equation}
\left\{
\begin{aligned}
&-\Delta u=\left(\frac{1}{|x|^{n-\alpha}}\ast v^p\right)v^{p-1} \quad \text{on} \ \Omega,\\
&-\Delta v=\left(\frac{1}{|x|^{n-\beta}}\ast u^q\right)u^{q-1} \quad \text{on} \ \Omega.
\end{aligned}
\right.
\end{equation}To obtain this result, we need a doubling lemma (cf. Lemma 2.2) which plays an important role in the study of the Lane–Emden conjecture. Polacik et al. [Reference Polacik, Quittner and Souplet37] proved that non-existence of bounded solutions of Equation (1.9) implies estimates of boundary blowing-up rate:
\begin{equation*}
u(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{2(p+1)}{pq-1}},\quad x\in\Omega,
\end{equation*}
\begin{equation*}
v(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{2(q+1)}{pq-1}},\quad x\in\Omega,
\end{equation*} where (u, v) solves the Lane–Emden system on the bounded domain Ω. Combining with the result of Serrin–Zou [Reference Serrin and Zou39], taking 
$\Omega=\partial B_R(x)$ and letting 
$R \rightarrow \infty$, they proved the Lane–Emden conjecture when n = 3. Therefore, we are interested in the boundary blowing-up rate of the system (1.15).
Theorem 1.4. Assume that Equation (1.1) with 
$\min\{p, q\} \gt 1$ has no bounded positive classical solution (u, v) in 
$\mathbb{R}^n$. Then, there exists 
$C=C(n,p,q) \gt 0$ such that any positive solution 
$(u,v)
\in C^2(\mathbb{R}^n) \times C^2(\mathbb{R}^n)$ of Equation (1.15) satisfies estimates of the boundary blow-up rates:
\begin{equation}
u(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{(2p-1)(\beta+2)+(\alpha+2)}{(2p-1)(2q-1)-1}},\quad x\in\Omega,
\end{equation}
\begin{equation}
v(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{(2q-1)(\alpha+2)+(\beta+2)}{(2p-1)(2q-1)-1}},\quad x\in\Omega.
\end{equation}On the contrary, if positive classical solutions of Equation (1.15) satisfy Equations (1.16) and (1.17), then Equation (1.1) has no positive classical solution.
Another analogous problem is
\begin{equation*}
\begin{cases}
0 \leq -\Delta u \leq (|x|^{-\alpha} * v)^\lambda,\\
0 \leq -\Delta v \leq (|x|^{-\beta} * u)^\sigma,
\end{cases}
\end{equation*}where 
$n \geq 3$, 
$\alpha,\beta \in (0,n)$ and 
$\lambda, \sigma \geq 0$. In 2015, Ghergu and Taliaferro studied the behaviour near the origin of positive solutions in 
$C^2(\mathbb{R}^n \setminus \{0\}) \cap L^1(\mathbb{R}^n)$ (cf. [Reference Ghergu and Taliaferro13]). This shows that asymptotic behaviour of the positive solution of Equation (1.1) is an interesting topic which will be investigated later.
2. Preliminaries
Recall a version of the Hardy–Littlewood–Sobolev inequality, which will be used in the method of moving planes in integral forms (introduced by Chen–Li–Ou [Reference Chen, Li and Ou5]).
Lemma 2.1. (Theorem 1 in Chapter 5 of [Reference Stein41])
Let 
$0 \lt \alpha \lt n$ and
\begin{equation*}f(x)=\int_{\mathbb{R}^n}\vert x-y\vert^{\alpha-n}g(y){}\mathrm dy.\end{equation*} Then, for any 
$s \gt \frac{n}{n-\alpha},$ we have
\begin{equation*}\|f\|_{L^s(\mathbb{R}^n)}\leq C(n,s,\alpha)\|g\|_{L^{\frac{ns}{n+\alpha s}}(\mathbb{R}^n)}.\end{equation*}Next we recall the doubling lemma by Polacik, Quittner and Souplet. Those ideas come from [Reference Hu19].
Lemma 2.2. (Lemma 5.1 in [Reference Polacik, Quittner and Souplet37], Doubling lemma)
Let (X, d) be a complete metric space, and let 
$\emptyset\neq D\subset\Sigma\subset X$ with Σ close. Set 
$\Gamma=\Sigma\setminus D.$ Finally, let 
$M: D\rightarrow(0,\infty)$ be bounded on compact subsets of D and fix a real k > 0. If 
$y\in D$ is such that
\begin{equation*}M(y) \mathrm{dist}(y, \Gamma) \gt 2k,\end{equation*}then there exists 
$x\in D$ such that
\begin{equation*}M(x) \mathrm{dist}(x, \Gamma) \gt 2k,\qquad M(x)\geq M(y),\end{equation*}and
\begin{equation*}M(z)\leq 2M(x),\quad {\rm for} \ {\rm all} \; z\in D\cap \overline{B_{k/M(x)}(x)}.\end{equation*}Now, we introduce a lemma which plays an important role in radial case, and it can be seen in [Reference Clément, Manasevich and Mitidieri6] or [Reference Mitidieri35].
Lemma 2.3. (Lemma 2.1 in [Reference Mitidieri35])
Assume n > 2 and 
$\Psi\in C^{2}(\mathbb{R}^{n}\setminus \{0\})$ is a positive radial function. Let
\begin{equation*}(-\Delta)^k\Psi\geq0,\quad {\rm in} \ \mathbb{R}^n,\; k=0,1,\end{equation*}then for every 
$r=|x| \in(0,\infty)$ we have
\begin{equation}
r\Psi^{\prime}(r)+(n-2)\Psi(r)\geq0.
\end{equation}Finally, we will use a result on the relation between Equation (1.6) and the system of partial differential equations (1.7) and the system of integral equations (1.8). The result can be found in [Reference Hayman and Kennedy18] (or [Reference Caristi, D’Ambrosio and Mitidieri2]).
Lemma 2.4. (Theorem 3.21 in [Reference Hayman and Kennedy18])
Let 
$n \geq 3$ and µ be a positive Radon measure on 
$\mathbb{R}^n$ and 
$l \in \mathbb{R}$. The following two statements are equivalent:
(a) u is a distributional solution of
$\Delta u+\mu=0$ on $\mathbb{R}^n$, and ess $\inf_{\mathbb{R}^n} u = l$.(b)
$u \in L^1_{\rm loc}(\mathbb{R}^n)$, and we have
\begin{equation*}
u(x) = l + c_*\int_{\mathbb{R}^n} \frac{{\rm d}\mu(y)}{|x-y|^{n-2}}, \quad \textrm{a.e.} \ x \in \mathbb{R}^n,
\end{equation*}
where 
$c_* := \Gamma((n-2)/2)(4\pi^{n/2}\Gamma(1))^{-1}$.
Further results can be found in reference [Reference D’ambrosio and Ghergu7]. D’ambrosio and Ghergu obtained their integral representation formulae for functions 
$u \in L^1_{\rm loc}
(\mathbb{R}^n)$, which satisfy 
$P(-\Delta)u =\mu$ in the sense of distributions, where P is a non-constant real non-homogeneous polynomial whose roots are non-positive. Those results can be applied to many non-homogeneous higher-order differential inequalities.
3. Liouville theorem under the Serrin condition
In this section, we give a necessary condition of existence of positive solutions, which provide an important ingredient in the proof of Theorem 5.1.
First, we say 
$(u,v,w,z)$ is a positive super-solution (sub-solution) of Equation (1.8) if 
$u,v,w,z \in L_{\rm loc}^1(\mathbb{R}^n)$ are positive so that the following inequalities make sense a.e. on 
$\mathbb{R}^n$
\begin{equation*}
\left\{
\begin{aligned}
&u(x) \geq (\leq) C_1\int_{\mathbb{R}^n}\frac{v^{p-1}(y)w(y)}{|x-y|^{n-2}}\,{\rm d}y,\\
&w(x) \geq (\leq) \int_{\mathbb{R}^n}\frac{v^p(y)}{|x-y|^{n-2}}\,{\rm d}y,\\
&v(x) \geq (\leq) C_2 \int_{\mathbb{R}^n}\frac{u^{q-1}(y)z(y)}{|x-y|^{n-2}}\,{\rm d}y,\\
&z(x) \geq (\leq) \int_{\mathbb{R}^n}\frac{u^q(y)}{|x-y|^{n-2}}\,{\rm d}y.
\end{aligned}
\right.
\end{equation*}Now, we use an idea in [Reference Caristi, D’Ambrosio and Mitidieri2] to prove the following theorem.
Theorem 3.1 Let 
$\min\{p,q\} \gt 1$ and
\begin{equation}
\max\left\{\frac{8p}{(2p-1)(2q-1)-1},
\frac{8q}{(2p-1)(2q-1)-1}\right\}\geq n-2.
\end{equation}Then, the integral system (1.8) has no positive super-solution.
Proof. If 
$(u(x), v(x))$ is a positive super-solution of Equation (1.8), we can deduce a contradiction. In fact, from the system (1.8), we have
\begin{equation}
u(x)\geq C_1\int_{B_R(0)}\frac{v^{p-1}(y)w(y)}{|x-y|^{n-2}}\,{\rm d}y\geq \frac{c}{(|x|+R)^{n-2}}\int_{B_R(0)}v^{p-1}(y)w(y)\,{\rm d}y.
\end{equation}
\begin{equation}
w(x)\geq \int_{B_R(0)}\frac{v^p(y)}{|x-y|^{n-2}}\,{\rm d}y\geq \frac{c}{(|x|+R)^{n-2}}\int_{B_R(0)}v^p(y)\,{\rm d}y.
\end{equation}
\begin{equation}
v(x)\geq C_2\int_{B_R(0)}\frac{u^{q-1}(y)z(y)}{|x-y|^{n-2}}\,{\rm d}y\geq \frac{c}{(|x|+R)^{n-2}}\int_{B_R(0)}u^{q-1}(y)z(y)\,{\rm d}y.
\end{equation}
\begin{equation}
z(x)\geq \int_{B_R(0)}\frac{u^q(y)}{|x-y|^{n-2}}\,{\rm d}y\geq \frac{c}{(|x|+R)^{n-2}}\int_{B_R(0)}u^q(y)\,{\rm d}y.
\end{equation} Taking q − 1 powers of Equation (3.2) and multiplying Equation (3.5) and then integrating on 
$B_R(0)$, we obtain
\begin{equation}
\begin{aligned}
&\quad \int_{B_R(0)}u^{q-1}(x)z(x)\,{\rm d}x\\
&\geq c \int_{B_R(0)}\frac{{\rm d}x}{(|x|+R)^{q(n-2)}}
\left(\int_{B_R(0)}v^{p-1}(y)w(y)\,{\rm d}y\right)^{q-1}\int_{B_R(0)}u^q(y)\,{\rm d}y\\
&\geq cR^{n-q(n-2)}\left(\int_{B_R(0)}v^{p-1}(y)w(y)\,{\rm d}y\right)^{q-1}\int_{B_R(0)}u^q(y)\,{\rm d}y.
\end{aligned}
\end{equation} Taking q powers of Equation (3.2) and integrating on 
$B_R(0)$, we get
\begin{equation*}\int_{B_R(0)}u^q(x)\,{\rm d}x\geq cR^{n-q(n-2)}\left(\int_{B_R(0)}v^{p-1}(y)w(y)\,{\rm d}y\right)^{q}.\end{equation*}Inserting this into Equation (3.6), we see that
\begin{equation}
\int_{B_R(0)}u^{q-1}(x)z(x)\,{\rm d}x\geq cR^{2n-2q(n-2)}\left(\int_{B_R(0)}v^{p-1}(y)w(y)\,{\rm d}y\right)^{2q-1}.
\end{equation}Similarly, we have
\begin{equation*}
\int_{B_R(0)}v^{p-1}(x)w(x)\,{\rm d}x\geq cR^{2n-2p(n-2)}\left(\int_{B_R(0)}u^{q-1}(y)z(y)\,{\rm d}y\right)^{2p-1}.
\end{equation*}Combining with Equation (3.7), there holds
\begin{equation}
\int_{B_R(0)}u^{q-1}(x)z(x)\,{\rm d}x\geq cR^{2q(n+2)-2p(2q-1)(n-2)}\left(\int_{B_R(0)}u^{q-1}(x)z(x)\,{\rm d}x\right)^{(2p-1)(2q-1)}.
\end{equation} In view of 
$(2p-1)(2q-1) \gt 1,$ the result above implies that
\begin{equation}
\int_{B_R(0)}u^{q-1}(x)z(x)\,{\rm d}x\leq cR^{-\frac{2q(n+2)-2p(2q-1)(n-2)}{(2p-1)(2q-1)-1}}.
\end{equation} When 
$\frac{8q}{(2p-1)(2q-1)-1} \gt n-2,$ there holds 
$2q(n+2)-2p(2q-1)(n-2) \gt 0.$ Letting 
$R\rightarrow\infty$ in Equation (3.9), we have 
$\int_{\mathbb{R}^n}u^{q-1}(x)z(x)\,{\rm d}x=0,$ which contradicts with 
$u, z \gt 0.$
When 
$\frac{8q}{(2p-1)(2q-1)-1}=n-2,$ there holds 
$2q(n+2)-2p(2q-1)(n-2)=0.$ Letting 
$R\rightarrow\infty$ in Equation (3.9), we have 
$u^{q-1}z\in L^1(\mathbb{R}^n).$ Similar to the derivation of Equation (3.8), we integrate on 
$A_R:=B_{2R}(0)\backslash B_R(0)$ instead of on 
$B_R(0)$. Thus,
\begin{equation*}
\int_{A_R}u^{q-1}(x)z(x)\,{\rm d}x\geq c\left(\int_{B_R(0)}u^{q-1}(x)z(x)\,{\rm d}x\right)^{(2p-1)(2q-1)}.
\end{equation*} Letting 
$R\rightarrow\infty$ and noting 
$u^{q-1}z\in L^1(\mathbb{R}^n),$ we see 
$\int_{\mathbb{R}^n}u^{q-1}(x)z(x)\,{\rm d}x=0.$ This is a contradiction.
In the same way, we can prove that Equation (1.8) has no positive super-solution if 
$\frac{8p}{(2p-1)(2q-1)-1}\geq n-2.$ Therefore, we complete the proof of Theorem 3.1.
Proof of Theorem 1.1
Let (u, v) be a pair of positive distributional solution of Equation (1.6) under the Serrin condition (3.1). Now, 
$\inf_{\mathbb{R}^n}u \geq 0$ and 
$\inf_{\mathbb{R}^n}v \geq 0$. According to Lemma 2.4, we have
\begin{equation*}
u(x) \geq c_*\int_{\mathbb{R}^n}\frac{w(y)v^{p-1}(y)\,{\rm d}y}{|x-y|^{n-2}}, \quad \textrm{a.e.} \ on \ \mathbb{R}^n,
\end{equation*}
\begin{equation*}
v(x) \geq c_*\int_{\mathbb{R}^n}\frac{z(y)u^{q-1}(y)\,{\rm d}y}{|x-y|^{n-2}} \quad \textrm{a.e.} \ on \ \mathbb{R}^n.
\end{equation*} Therefore, 
$(u,v,w,z)$ is a super-solution of Equation (1.8) with 
$C_1=C_2=c_*$. This contradicts with Theorem 3.1. Thus, Equation (1.6) has no positive solution under the Serrin condition (3.1).
Remark 3.1. By the same proofs of Theorem 3.1 and 1.1, we can also obtain that the Serrin-type condition of Equation (1.7) is Equation (3.1), and the Serrin-type condition of Equation (1.1) is
\begin{equation*}
\max\left\{\frac{\alpha+\beta(2p-1)+4p}{(2p-1)(2q-1)-1},
\frac{\beta+\alpha(2q-1)+4q}{(2p-1)(2q-1)-1}\right\} \geq n-2.
\end{equation*}4. Liouville theorems in subcritical case
In this section, we prove Theorem 1.2. Namely, we prove that Equation (1.6) has no positive radial classical solution in the Sobolev-type subcritical case. The ideas in [Reference Mitidieri35] are employed here. In fact, Mitidieri proved this non-existence by a contradiction argument. Assume Equation (1.9) has a pair of positive radial classical solution, one multiplies equations by the normal derivatives of solutions and integrates on a ball. Integrating by parts and combining them together, one can deduce a Pohozaev-type identity. In order to handle integrals on the boundary of the ball, one need to estimate decay rates of solutions at infinity where Lemma 2.3 plays an important role.
Now, we use the ideas in [Reference Mitidieri35] to deal with the non-existence of radial classical solutions of Equation (1.6). We first use those ideas to prove the non-existence of positive radial classical solutions of Equation (1.7) (rather than Equation (1.6)). The reason is that the convolution terms are not easy to handle when deducing directly the Pohozaev identity from Equation (1.6). Even if a Pohozaev identity of integral form can be derived from 
$w(x)=|x|^{2-n}\times v^p$ and 
$z(x)=|x|^{2-n}\times u^q$, new improper integrals will appear and their convergence is difficult to prove. Next, if (u, v) is a pair of classical solution of Equation (1.6), by the regularity theory of singular integrals (cf. § 4.2 in [Reference Gilbarg and Trudinger16]), from the Hölder continuity of u and v, we can derive the second-order differentiability of w and z, and hence 
$(u,v,w,z)$ is the classical solution of Equation (1.7). Thus, we can draw the desired conclusion.
Theorem 4.1. In the subcritical case (1.12), Equation (1.7) has no positive radial classical solutions.
Proof. If Equation (1.7) has positive radial solutions 
$(u, w, v, z)$, we can deduce a contradiction.
In fact, writing Equation (1.7) in radial coordinates, we obtain for 
$r \gt 0,$
\begin{equation*}-(ru^{\prime}(r)+(n-2)u(r))^{\prime}=rw(r)v^{p-1}(r),\qquad -(rw^{\prime}(r)+(n-2)w(r))^{\prime}=rv^{p}(r),\end{equation*}
\begin{equation*}-(rv^{\prime}(r)+(n-2)v(r))^{\prime}=rz(r)u^{q-1}(r),\qquad -(rz^{\prime}(r)+(n-2)z(r))^{\prime}=ru^{q}(r).\end{equation*} The first equation shows 
$(r^{n-2}u^{\prime})^{\prime} \lt 0$. Integrating from 0 to r yields 
$u^{\prime}(r) \lt 0$ for all r > 0. Similarly, 
$v^{\prime},w^{\prime},z^{\prime}$ are also negative for r > 0.
According to Equation (2.1), we have 
$(v(r)r^{n-2})^{\prime}, (w(r)r^{n-2})^{\prime}\geq0.$ Integrating the radial equations from s to t for 
$0 \lt s\leq t,$ we see that
\begin{equation*}
\begin{aligned}
su^{\prime}(s)+(n-2)u(s)&\geq w(s)s^{n-2}(v(s)s^{n-2})^{p-1}\int_s^t\xi^{1-p(n-2)}{\rm d}\xi\\
&\geq w(s)v^{p-1}(s)s^{p(n-2)+1}\frac{1}{1-p(n-2)}(t^{1-p(n-2)}-s^{1-p(n-2)}),
\end{aligned}
\end{equation*}and
\begin{equation*}
\begin{aligned}
sw^{\prime}(s)+(n-2)w(s)&\geq (v(s)s^{n-2})^{p}\int_s^t\xi^{1-p(n-2)}{\rm d}\xi\\
&\geq v^{p}(s)s^{p(n-2)+1}\frac{1}{1-p(n-2)}(t^{1-p(n-2)}-s^{1-p(n-2)}).
\end{aligned}
\end{equation*} Since 
$u^{\prime}, w^{\prime} \lt 0$ and 
$1-p(n-2) \lt 0$, we can see that for 
$r\geq r_0 \gt 0,$
\begin{equation*}u(r)\geq cr^2w(r)v^{p-1}(r),\qquad w(r)\geq cr^2v^p(r).\end{equation*}Similarly, we have
\begin{equation*}v(r)\geq cr^2z(r)u^{q-1}(r),\qquad z(r)\geq cr^2u^q(r).\end{equation*}Hence,
\begin{equation}
u(r)\leq cr^{-\frac{8p}{(2p-1)(2q-1)-1}},\qquad v(r)\leq cr^{-\frac{8q}{(2p-1)(2q-1)-1}},
\end{equation}and then
\begin{equation}
w(r)v^p(r),\; z(r)u^q(r)\leq cr^{-2}u(r)v(r)\leq cr^{-2\frac{(2p+1)(2q+1)-1}{(2p-1)(2q-1)-1}}.
\end{equation}According to the identity Equation (2.5) in [Reference Mitidieri34] (see also Equation (3.5) in [Reference Mitidieri35]), there holds
\begin{equation}
\begin{aligned}
&\quad-\int_0^Rv^{p-1}(r)v^{\prime}(r)w(r)r^n\,{\rm d}r-\int_0^Ru^{q-1}(r)u^{\prime}(r)z(r)r^n\,{\rm d}r\\
&=(n-2)\int_0^Ru^{\prime}(r)v^{\prime}(r)r^{n-1}\,{\rm d}r+R^nu^{\prime}(r)v^{\prime}(r).
\end{aligned}
\end{equation} Multiplying Equation (1.7)1 by v and Equation (1.7)3 by u, and integrating by parts on 
$(0,R)$, we get
\begin{equation}
-R^{n-1}u^{\prime}(R)v(R)+\int_0^Rr^{n-1}u^{\prime}(r)v^{\prime}(r)\,{\rm d}r=\int_0^Rr^{n-1}w(r)v^p(r)\,{\rm d}r,
\end{equation}
\begin{equation}
-R^{n-1}v^{\prime}(R)u(R)+\int_0^Rr^{n-1}u^{\prime}(r)v^{\prime}(r)\,{\rm d}r=\int_0^Rr^{n-1}z(r)u^q(r)\,{\rm d}r.
\end{equation}We claim that
\begin{equation}
\lim\limits_{R\rightarrow\infty}R^{n-1}u^{\prime}(R)v(R)=\lim\limits_{R\rightarrow\infty}R^{n-1}v^{\prime}(R)u(R)=0.
\end{equation}In fact, from Equations (2.1) and (4.1), it follows that
\begin{equation*}
\begin{aligned}
|R^{n-1}u^{\prime}(R)v(R)|,\; |R^{n-1}v^{\prime}(R)u(R)|&\leq (n-2)u(R)v(R)R^{n-2}\\
&\leq cR^{n-2-\frac{8(p+q)}{(2p-1)(2q-1)-1}}.
\end{aligned}
\end{equation*} In view of 
$n-2-\frac{8(p+q)}{(2p-1)(2q-1)-1} \lt 0$ (implied by Equation (1.12)), Equation (4.6) is true.
Using Equation (4.2), we also get
\begin{equation*}\int_0^\infty r^{n-1}w(r)v^p(r)\,{\rm d}r \lt \infty,\qquad \int_0^\infty r^{n-1}z(r)u^q(r)\,{\rm d}r \lt \infty.\end{equation*}Hence, from Equation (4.4)–(4.6), there holds
\begin{equation}
\int_0^\infty r^{n-1}u^{\prime}(r)v^{\prime}(r)\,{\rm d}r=\int_0^\infty r^{n-1}w(r)v^p(r)\,{\rm d}r=\int_0^\infty r^{n-1}z(r)u^q(r)\,{\rm d}r.
\end{equation} Integrating the left-hand side of Equation (4.3) by parts on 
$(0,R)$ yields
\begin{equation}
\begin{aligned}
-&\int_0^Rv^{p-1}(r)v^{\prime}(r)w(r)r^n\,{\rm d}r\\
&=-\frac{1}{p}w(R)v^p(R)R^n+\frac{n}{p}\int_0^Rw(r)v^p(r)r^{n-1\,}{\rm d}r+\frac{1}{p}\int_0^Rw^{\prime}(r)v^p(r)r^{n}\,{\rm d}r,
\end{aligned}
\end{equation}and
\begin{equation}
\begin{aligned}
-&\int_0^Ru^{q-1}(r)u^{\prime}(r)z(r)r^n\,{\rm d}r\\
&=-\frac{1}{q}z(R)u^q(R)R^n+\frac{n}{q}\int_0^Rz(r)u^q(r)r^{n-1}\,{\rm d}r+\frac{1}{q}\int_0^Rz^{\prime}(r)u^q(r)r^{n}\,{\rm d}r.
\end{aligned}
\end{equation} For any 
$\theta\in \mathbb{R},$ from Equations (4.4) and (4.5), it follows that
\begin{equation*}
\begin{aligned}
\int_0^Rr^{n-1}&u^{\prime}(r)v^{\prime}(r)\,{\rm d}r=\theta R^{n-1}u^{\prime}(R)v(R)+(1-\theta)R^{n-1}v^{\prime}(R)u(R)\\
&+\theta\int_0^Rr^{n-1}w(r)v^p(r)\,{\rm d}r+(1-\theta)\int_0^Rr^{n-1}z(r)u^q(r)\,{\rm d}r.
\end{aligned}
\end{equation*}Combining with Equations (4.8), (4.9) and (4.3), we obtain that
\begin{equation}
\begin{aligned}
&\quad \frac{1}{p}\int_0^Rw^{\prime}(r)v^p(r)r^{n}\,{\rm d}r+\frac{1}{q}\int_0^Rz^{\prime}(r)u^q(r)r^{n}\,{\rm d}r\\
=&\left[(n-2)\theta-\frac{n}{p}\right]\int_0^Rw(r)v^p(r)r^{n-1\,}{\rm d}r
+\left[ (n-2)(1-\theta)-\frac{n}{q}\right]\int_0^Rz(r)u^q(r)r^{n-1}\,{\rm d}r\\
&+\frac{1}{p}w(R)v^p(R)R^n+\frac{1}{q}z(R)u^q(R)R^n+R^nu^{\prime}(r)v^{\prime}(r)\\
&\quad+\theta(n-2) R^{n-1}u^{\prime}(R)v(R)+(1-\theta)(n-2)R^{n-1}v^{\prime}(R)u(R).
\end{aligned}
\end{equation} Let 
$0 \lt \theta \lt 1$. In view of 
$u^{\prime}(r), v^{\prime}(r) \lt 0,$ 
$ru^{\prime}(r)+(n-2)u(r)\geq0$ and 
$rv^{\prime}(r)+(n-2)v(r)\geq0$, there holds
\begin{equation}
\begin{aligned}
&\qquad R^nu^{\prime}(r)v^{\prime}(r)+\theta(n-2) R^{n-1}u^{\prime}(R)v(R)+(1-\theta)(n-2)R^{n-1}v^{\prime}(R)u(R)\\
&=\theta R^{n-1}u^{\prime}(R)[Rv^{\prime}(r)+(n-2)v(R)]+(1-\theta)R^{n-1}v^{\prime}(R)[Ru^{\prime}(r)+(n-2)u(R)]\leq0.
\end{aligned}
\end{equation}Therefore, Equation (4.10) reduces to
\begin{equation}
\begin{aligned}
& \frac{1}{p}\int_0^Rw^{\prime}(r)v^p(r)r^{n}\,{\rm d}r+\frac{1}{q}\int_0^Rz^{\prime}(r)u^q(r)r^{n}\,{\rm d}r\\
\leq&\left[(n-2)\theta-\frac{n}{p}\right]\int_0^Rw(r)v^p(r)r^{n-1\,}{\rm d}r
+\left[ (n-2)(1-\theta)-\frac{n}{q}\right]\int_0^Rz(r)u^q(r)r^{n-1\,}{\rm d}r\\
&\quad \quad +\frac{1}{p}w(R)v^p(R)R^n+\frac{1}{q}z(R)u^q(R)R^n.
\end{aligned}
\end{equation}Next, we claim that
\begin{equation}
\int_0^Rw^{\prime}(r)v^p(r)r^{n\,}{\rm d}r\geq \frac{2-n}{2}\int_0^Rv^p(r)w(r)r^{n-1\,}{\rm d}r.
\end{equation}Indeed, from Equation (1.7)2, it follows that
\begin{equation*}
-(r^{n-1}w^{\prime}(r))^{\prime}=v^p(r)r^{n-1}.
\end{equation*} Multiplying by rwʹ and integrating on 
$(0, R)$ yields
\begin{equation}
\begin{aligned}
\int_0^R&v^p(r)w^{\prime}(r)r^{n}\,{\rm d}r=-\int_0^R(r^{n-1}w^{\prime}(r))^{\prime}rw^{\prime}(r)\,{\rm d}r\\
&=-R^n(w^{\prime}(R))^2+\int_0^Rr^nw^{\prime}(r)w^{\prime\prime}(r)\,{\rm d}r+\int_0^Rr^{n-1}(w^{\prime}(r))^2\,{\rm d}r.
\end{aligned}
\end{equation}To handle the second term of the right-hand side, we notice that
\begin{equation*}
\begin{aligned}
\int_0^Rr^nw^{\prime}(r)w^{\prime\prime}(r)\,{\rm d}r&=\int_0^Rr^nw^{\prime}(r)\,{\rm d}w^{\prime}=R^n(w^{\prime}(R))^2-\int_0^Rw^{\prime}(r)(r^nw^{\prime}(r))^{\prime}\,{\rm d}r\\
&=R^n(w^{\prime}(R))^2-n\int_0^Rr^{n-1}(w^{\prime}(r))^2\,{\rm d}r-\int_0^Rr^nw^{\prime}(r)w^{\prime\prime}(r)\,{\rm d}r.
\end{aligned}
\end{equation*}Therefore,
\begin{equation}
\int_0^Rr^nw^{\prime}(r)w^{\prime\prime}(r)\,{\rm d}r=\frac{1}{2}R^n(w^{\prime}(R))^2-\frac{n}{2}\int_0^Rr^{n-1}(w^{\prime}(r))^2\,{\rm d}r.
\end{equation}To handle the third term of the right-hand side, we notice that
\begin{equation*}
\begin{aligned}
\int_0^Rv^p(r)w(r)r^{n-1\,}{\rm d}r&=-\int_0^R(r^{n-1}w^{\prime}(r))^{\prime}w(r)\,{\rm d}r\\
&=-R^{n-1}w^{\prime}(R)w(R)+\int_0^Rr^{n-1}(w^{\prime}(r))^2\,{\rm d}r.
\end{aligned}
\end{equation*}Combining this result with Equations (4.14) and (4.15), we get
\begin{equation*}
\begin{aligned}
&\int_0^Rv^p(r)w^{\prime}(r)r^{n\,}{\rm d}r\\
&=-\frac{1}{2}R^n(w^{\prime}(R))^2+\frac{2-n}{2}\int_0^Rv^p(r)w(r)r^{n-1\,}{\rm d}r+\frac{2-n}{2}R^{n-1}w^{\prime}(R)w(R)\\
&=-\frac{1}{2}R^{n-1}w^{\prime}(R)[Rw^{\prime}(R)+(n-2)w(R)]+\frac{2-n}{2}\int_0^Rv^p(r)w(r)r^{n-1}\,{\rm d}r\\
&\geq \frac{2-n}{2}\int_0^Rv^p(r)w(r)r^{n-1}\,{\rm d}r.
\end{aligned}
\end{equation*} Here, we use the fact of 
$w^{\prime}(R) \lt 0$ and 
$Rw^{\prime}(R)+(n-2)w(R)\geq0$ (implied by Equation (2.1)).
Similarly, we can also obtain
\begin{equation}
\int_0^Rz^{\prime}(r)u^q(r)r^{n}\,{\rm d}r\geq \frac{2-n}{2}\int_0^Ru^q(r)z(r)r^{n-1}\,{\rm d}r.
\end{equation} Letting 
$R\rightarrow\infty$ in Equation (4.12) and using Equations (4.2), (4.7), (4.13) and (4.16), we obtain
\begin{equation*}
0\leq \left\{n-2-\frac{n}{p}-\frac{n}{q}+\frac{n-2}{2}\cdot\frac{1}{p}
+\frac{n-2}{2}\cdot\frac{1}{q}\right\}\int_0^\infty v^p(r)w(r)r^{n-1\,}{\rm d}r.
\end{equation*}This contradicts with Equation (1.12), and hence Theorem 4.1 is proved.
Remark 4.1. On the contrary, when the subcritical case (1.12) is not true, i.e., 
$\frac{1}{p}+\frac{1}{q}\leq\frac{2(n-2)}{n+2},$ we can verify that
\begin{equation*}U(x)=a|x|^{-\frac{8p}{(2p-1)(2q-1)-1}},\qquad V(x)=b|x|^{-\frac{8q}{(2p-1)(2q-1)-1}},\end{equation*}
\begin{equation*}W(x)=c|x|^{-\frac{4(p+q)}{(2p-1)(2q-1)-1}},\qquad Z(x)=d|x|^{-\frac{4(p+q)}{(2p-1)(2q-1)-1}}\end{equation*} solve Equation (1.7) in 
$\mathbb{R}^n\setminus\{0\}$ for suitable 
$a,b,c,d \gt 0$. In addition, the classification result in [Reference Le24] shows that Equation (1.7) has an explicit radial solution on 
$\mathbb{R}^n$ when 
$p=q=\frac{n+2}{n-2}$.
Proof of Theorem 1.2
If Equation (1.6) has positive radial classical solution 
$u,v \in C^2(\mathbb{R}^n)$, we can use Lemma 2.4 to obtain that 
$(u,v,w,z)$ solves Equation (1.7) in distribution sense, where 
$w:=|x|^{2-n}*v^p$ and 
$z:=|x|^{2-n}*u^q$. Similar to the argument of regularity of the Newton potential in Section 4.2 of [Reference Gilbarg and Trudinger16], from 
$u,v \in C^2(\mathbb{R}^n),$ we can also deduce that 
$w,z
\in C^2(\mathbb{R}^n)$. Therefore, 
$(u,v,w,z)$ is the classical solution of Equation (1.7). This contradicts with Theorem 4.1.
5. Radial symmetry of integrable solutions
In this section, we employ the method of moving planes in integral forms introduced by Chen–Li–Ou [Reference Chen, Li and Ou5] to prove the radial symmetry of positive solutions of Equation (1.8). The methods of moving planes were founded by Alexanderoff in the early 1950s. Later, it was further developed by Serrin [Reference Serrin38], Gidas et al. [Reference Gidas, Ni and Nirenberg14], Caffarelli et al. [Reference Caffarelli, Gidas and Spruck1], Chen and Li [Reference Chen and Li3], Li and Zhu [Reference Li and Zhu29] and many others. Wang–Yang used this method to prove the radial symmetry of positive integrable solutions of Equation (1.1) (cf. [Reference Wang and Yang42]). Instead of the integrability condition in [Reference Wang and Yang42], we will consider another integrability condition (i.e., 
$(u,v)
\in L^{r_0}(\mathbb{R}^n)\times L^{s_0}(\mathbb{R}^n)$) to obtain the radial symmetry of positive solutions of Equation (1.8). In addition, we do not need the assumption that w(x) and z(x) are integrable.
Theorem 5.1. Let 
$(u, v)\in L^{r_0}(\mathbb{R}^n)\times L^{s_0}(\mathbb{R}^n)$ be a pair of positive solutions of Equation (1.8). Then (u, v) are radially symmetric and monotone decreasing about some point in 
$\mathbb{R}^n$.
Proof. First, we introduce some notation. For a given real number λ, let
\begin{equation*}\Sigma_\lambda=\{x=(x_1,x_2,\ldots,x_n)\mid x_1 \lt \lambda\},\end{equation*}and 
$x^\lambda=(2\lambda-x_1,x_2,\ldots,x_n)$ be the reflection point of x about the plane 
$x_1=\lambda$. Write
\begin{equation*}\Sigma_\lambda^u:=\{x\in \Sigma_\lambda\mid u(x) \gt u_\lambda(x):=u(x^\lambda)\},
\qquad \Sigma_\lambda^w:=\{x\in \Sigma_\lambda\mid w(x) \gt w_\lambda(x):=w(x^\lambda)\},\end{equation*}
\begin{equation*}\Sigma_\lambda^v:=\{x\in \Sigma_\lambda\mid v(x) \gt v_\lambda(x):=u(x^\lambda)\},
\qquad \Sigma_\lambda^z:=\{x\in \Sigma_\lambda\mid z(x) \gt z_\lambda(x):=z(x^\lambda)\}.\end{equation*} Assume (u, v) is a pair of positive solutions of Equation (1.8). Write 
$t=\frac{n[(2p-1)(2q-1)-1]}{4(p+q)}.$ According to Theorem 1.1, we know that
\begin{equation}
\max\left\{\frac{8p}{(2p-1)(2q-1)-1}, \frac{8q}{(2p-1)(2q-1)-1}\right\} \lt n-2,
\end{equation}which implies 
$t \gt \frac{n}{n-2}$ (due to 
$2\,\max\{p,q\} \gt p+q$). Therefore, by Lemma 2.1, we can deduce 
$w,z \in L^t(\mathbb{R}^n)$ from 
$(u,v) \in
L^{r_0}(\mathbb{R}^n) \times L^{s_0}(\mathbb{R}^n)$.
Step 1. We show that for λ sufficiently negative,
\begin{equation}
u_\lambda(x)\geq u(x),\qquad v_\lambda(x)\geq v(x) \quad {\rm for} \ {\rm all} \; x\in \Sigma_\lambda.
\end{equation} To show Equation (5.2), we will prove that 
$\Sigma_\lambda^u$ and 
$\Sigma_\lambda^v$ must have measure zero for λ sufficiently negative.
First, by the mean value theorem and the fact that for any 
$0 \lt a\leq b,$ 
$r \gt 0,$
\begin{equation*}a^r-b^r\geq \max\{r, 1\}b^{r-1}(a-b),\end{equation*}it follows that
\begin{equation*}
\begin{aligned}
v^{p-1}(y)&w(y)-v_\lambda^{p-1}(y)w_\lambda(y)
=[v^{p-1}(y)-v_\lambda^{p-1}(y)]w(y)+v_\lambda^{p-1}(y)[w(y)-w_\lambda(y)]\\
&\leq \max\{p-1, 1\}w(y)v^{p-2}(y)[v(y)-v_\lambda(y)]^++v_\lambda^{p-1}(y)[w(y)-w_\lambda(y)]^+
\end{aligned}
\end{equation*}and
\begin{equation*}v^{p}(y)-v_\lambda^{p}(y)\leq pv^{p-1}(y)[v(y)-v_\lambda(y)]^+.\end{equation*} Here, we denote 
$f^+=\max\{f, 0\}.$
Therefore, for 
$x\in \Sigma_\lambda^u$,
\begin{equation}
\begin{aligned}
0& \lt u(x)-u_\lambda(x)\\
&=C_1\int_{\Sigma_\lambda}\left(\frac{1}{|x-y|^{n-2}}
-\frac{1}{|x^\lambda-y|^{n-2}}\right)[v^{p-1}(y)w(y)-v_\lambda^{p-1}(y)w_\lambda(y)]\,{\rm d}y\\
&\leq C_1\,\max\{p-1, 1\}\int_{\Sigma_\lambda}\frac{1}{|x-y|^{n-2}}w(y)v^{p-2}(y)[v(y)-v_\lambda(y)]^+\,{\rm d}y\\
&\quad\qquad +C_1\int_{\Sigma_\lambda}\frac{1}{|x-y|^{n-2}}v_\lambda^{p-1}(y)[w(y)-w_\lambda(y)]^+\,{\rm d}y,
\end{aligned}
\end{equation}and for 
$x\in \Sigma_\lambda^w$,
\begin{equation}
\begin{aligned}
0 \lt w(x)-w_\lambda(x)&=\int_{\Sigma_\lambda}\left(\frac{1}{|x-y|^{n-2}}
-\frac{1}{|x^\lambda-y|^{n-2}}\right)[v^{p}(y)-v_\lambda^{p}(y)]\,{\rm d}y\\
&\leq p\int_{\Sigma_\lambda}\frac{1}{|x-y|^{n-2}}v^{p-1}(y)[v(y)-v_\lambda(y)]^+\,{\rm d}y.
\end{aligned}
\end{equation}Applying Lemma 2.1 and the Hölder inequalities, we obtain
\begin{equation}
\begin{aligned}
&\quad \parallel u-u_\lambda\parallel_{L^{r_0}(\Sigma_\lambda^u)}\\
&\leq c\parallel wv^{p-2}(v-v_\lambda)\parallel_{L^{\frac{nr_0}{n+2r_0}}(\Sigma_\lambda^v)}
+c\parallel v_\lambda^{p-1}(w-w_\lambda)\parallel_{L^{\frac{nr_0}{n+2r_0}}(\Sigma_\lambda^w)}\\
&\leq c\parallel w\parallel_{L^{t}(\Sigma_\lambda^v)}\parallel v\parallel^{p-2}_{L^{s_0}(\Sigma_\lambda^v)}
\parallel v-v_\lambda\parallel_{L^{s_0}(\Sigma_\lambda^v)}
+c\parallel v\parallel^{p-1}_{L^{s_0}(\mathbb{R}^n)}\parallel w-w_\lambda\parallel_{L^{t}(\Sigma_\lambda^w)}\\
&\leq c\parallel v\parallel^{2p-2}_{L^{s_0}(\Sigma_\lambda^v)}\parallel v-v_\lambda\parallel_{L^{s_0}(\Sigma_\lambda^v)}
+c\parallel v\parallel^{p-1}_{L^{s_0}(\mathbb{R}^n)}\parallel w-w_\lambda\parallel_{L^{t}(\Sigma_\lambda^w)},
\end{aligned}
\end{equation}and
\begin{equation}
\begin{aligned}
\parallel w-w_\lambda\parallel_{L^{t}(\Sigma_\lambda^w)}&\leq
c\parallel v^{p-1}(v-v_\lambda)\parallel_{L^{\frac{nt}{n+2t}}(\Sigma_\lambda^v)}
\leq c\parallel v\parallel^{p-1}_{L^{s_0}(\Sigma_\lambda^v)}\parallel v-v_\lambda\parallel_{L^{s_0}(\Sigma_\lambda^v)}.
\end{aligned}
\end{equation}It is easy to verify that
\begin{equation*}\frac{2}{n}+\frac{1}{t}=\frac{p}{s_0}=\frac{q}{r_0}.\end{equation*} Furthermore, by Equation (5.1), we see 
$r_0, s_0 \gt \frac{n}{n-2}$. Therefore, Lemma 2.1 can be used here.
Therefore, combining Equations (5.5) and (5.6) yields
\begin{equation}
\parallel u-u_\lambda\parallel_{L^{r_0}(\Sigma_\lambda^u)}
\leq c\left(\parallel v\parallel^{2p-2}_{L^{s_0}(\Sigma_\lambda^v)}
+\parallel v\parallel^{p-1}_{L^{s_0}(\mathbb{R}^n)}\parallel v\parallel^{p-1}_{L^{s_0}(\Sigma_\lambda^v)}\right)
\parallel v-v_\lambda\parallel_{L^{s_0}(\Sigma_\lambda^v)}.
\end{equation}Similarly, we have,
\begin{equation}
\parallel v-v_\lambda\parallel_{L^{s_0}(\Sigma_\lambda^v)}
\leq c\left(\parallel u\parallel^{2q-2}_{L^{r_0}(\Sigma_\lambda^u)}
+\parallel u\parallel^{q-1}_{L^{r_0}(\mathbb{R}^n)}\parallel u\parallel^{q-1}_{L^{r_0}(\Sigma_\lambda^u)}\right)
\parallel u-u_\lambda\parallel_{L^{r_0}(\Sigma_\lambda^u)}.
\end{equation} By the integrability condition 
$(u, v) \in L^{r_0}(\mathbb{R}^n)\times L^{s_0}(\mathbb{R}^n)$, for sufficiently negative λ, we arrive at
\begin{equation}
c\left(\parallel v\parallel^{2p-2}_{L^{s_0}(\Sigma_\lambda^v)}
+\parallel v\parallel^{p-1}_{L^{s_0}(\mathbb{R}^n)}\parallel v\parallel^{p-1}_{L^{s_0}(\Sigma_\lambda^v)}\right)
\leq \frac{1}{4},
\end{equation}and
\begin{equation}
c\left(\parallel u\parallel^{2q-2}_{L^{r_0}(\Sigma_\lambda^u)}
+\parallel u\parallel^{q-1}_{L^{r_0}(\mathbb{R}^n)}\parallel u\parallel^{q-1}_{L^{r_0}(\Sigma_\lambda^u)}\right)
\leq \frac{1}{4}.
\end{equation}It follows from Equations (5.7) and (5.8) that
\begin{equation*}\parallel u-u_\lambda\parallel_{L^{r_0}(\Sigma_\lambda^u)}=0,
\qquad \parallel v-v_\lambda\parallel_{L^{s_0}(\Sigma_\lambda^v)}=0;\end{equation*}hence, 
$\Sigma_\lambda^{u}$ and 
$\Sigma_\lambda^{v}$ must have measure zero. This completes Step 1.
Step 2. We move the plane 
$x_1=\lambda$ to the right as long as Equation (5.2) holds. Define
\begin{equation*}\lambda_0=\sup\{\mu\mid \rm{Equation}~(5.2) \;{\rm holds}\ {\rm for} \ {\rm any}\ \lambda\leq \mu\}.\end{equation*} Using a similar argument as in Step 1, one can see that 
$\lambda_0 \lt \infty$. Then we claim that
\begin{equation}
u(x)\equiv u_{\lambda_0}(x),\qquad v(x)\equiv v_{\lambda_0}(x)\quad \forall \; x\in \Sigma_{\lambda_0}.
\end{equation} Otherwise, we can move the plane further to the right. Indeed, if 
$v(x)\equiv v_{\lambda_0}(x)$ is not true, from the equalities in Equations (5.4) and (5.3), we deduce 
$u_{\lambda_0}(x) \gt u(x)$ in 
$\Sigma_{\lambda_0}$. Similarly, 
$v_{\lambda_0}(x) \gt v(x)$ in 
$\Sigma_{\lambda_0}.$ Write
\begin{equation*}
\widetilde{\Sigma_{\lambda_0}^u}=\{x\in \Sigma_{\lambda_0}\mid u(x)\geq u_{\lambda_0}(x)\},
\qquad \widetilde{\Sigma_{\lambda_0}^v}=\{x\in \Sigma_{\lambda_0}\mid v(x)\geq v_{\lambda_0}(x)\}.
\end{equation*} Then, obviously we have 
$\widetilde{\Sigma_{\lambda_0}^u}$ has measure zero, and 
$\lim_{\lambda\rightarrow\lambda_0^+}\Sigma_\lambda^u\subset\widetilde{\Sigma_{\lambda_0}^u}.$ The same is true for that of 
$v.$
By means of the integrability conditions 
$u\in L^{r_0}(\mathbb{R}^n)$ and 
$v\in L^{s_0}(\mathbb{R}^n)$, we can choose ɛ sufficiently small such that Equations (5.9) and (5.10) hold for all 
$\lambda\in[\lambda_0,\lambda_0+\varepsilon)$. Therefore, we have
\begin{equation*}
\parallel u-u_\lambda\parallel_{L^{r_0}(\Sigma_\lambda^u)}
\leq \frac{1}{2} \parallel u-u_\lambda\parallel_{L^{r_0}(\Sigma_\lambda^u)},
\end{equation*}which implies 
$\Sigma_\lambda^u$ must be measure zero. Similarly, 
$\Sigma_\lambda^v$ must also be measure zero. This contradicts with the definition of λ 0, and hence Equation (5.11) is proved.
Since the x 1 direction can be chosen arbitrarily, we deduce that u and v must be radially symmetric and decreasing about some point in 
$\mathbb{R}^n.$ This completes the proof of Theorem 5.1.
Remark 5.1. If (u, v) solves Equation (1.8) and u, v are radially symmetric, we claim that w must be radially symmetric. It can be easily seen from Equation (5.4) with 
$\lambda=\lambda_0$. Similarly, z is also radially symmetric.
From Theorem 4.1 and Theorem 5.1, we can prove Theorem 1.3.
Proof of Theorem 1.3
Let 
$(u,v) \in L^{r_0}(\mathbb{R}^n) \times L^{s_0}(\mathbb{R}^n)$, solve Equation (1.6) in classical sense (which implies that (u, v) also solves Equation (1.6) in distributional sense). The integrability of (u, v) leads to 
$\inf_{\mathbb{R}^n}u=\inf_{\mathbb{R}^n}v=0$. According to Lemma 2.4, we have
\begin{equation*}
u(x) = c_*\int_{\mathbb{R}^n}\frac{w(y)v^{p-1}(y)\,{\rm d}y}{|x-y|^{n-2}},
\quad
v(x) = c_*\int_{\mathbb{R}^n}\frac{z(y)u^{q-1}(y)\,{\rm d}y}{|x-y|^{n-2}},
\end{equation*}where 
$w=|x|^{2-n}*v^p$ and 
$z=|x|^{2-n}*u^q$. Therefore, 
$(u,v,w,z)$ is a solution of Equation (1.8) with 
$C_1=C_2=c_*$. According to Theorem 5.1, (u, v) is a pair of positive radial classical solution of Equation (1.6). Therefore, by Theorem 1.2, Equation (1.6) has no classical solution in 
$L^{r_0}(\mathbb{R}^n) \times L^{s_0}(\mathbb{R}^n)$ when Equation (1.12) holds.
Denote (1.8) with 
$C_1=C_2=c_\ast$ and with 2 replaced by 
$\alpha$ in 
$w$ and replaced by 
$\beta$ in 
$z$ by (1.8)'.
Corollary 5.2. Assume that Equation (1.8)’ has no positive radial solution. Then, Equation (1.1) has no positive radial solution in 
$C^2(\mathbb{R}^n) \times C^2(\mathbb{R}^n)$.
Proof. If 
$(u,v) \in C^2(\mathbb{R}^n) \times C^2(\mathbb{R}^n)$ is a pair of positive radial solution of Equation (1.1), we can deduce a contradiction.
In fact, if we write 
$w=|x|^{2-n}*v^p$ and 
$z=|x|^{2-n}*u^q$, w and z are positive. In addition, from Equation (1.1), it follows that
\begin{equation*}
w=-v^{1-p}\Delta u, \qquad z=-u^{1-q}\Delta v,
\end{equation*}which imply that 
$w,z$ are also radial.
According to Lemma 2.4, there holds
\begin{equation*}
u(x) =\inf_{\mathbb{R}^n}u+ c_*\int_{\mathbb{R}^n}\frac{w(y)v^{p-1}(y)\,{\rm d}y}{|x-y|^{n-2}},
\qquad
v(x) =\inf_{\mathbb{R}^n}v+ c_*\int_{\mathbb{R}^n}\frac{z(y)u^{q-1}(y)\,{\rm d}y}{|x-y|^{n-2}}.
\end{equation*} We claim 
$\inf_{\mathbb{R}^n}u=0$. Otherwise, we can find c > 0 such that 
$\inf_{\mathbb{R}^n}u \geq c$. Therefore, for any 
$x_0 \in \mathbb{R}^n$,
\begin{equation*}
z(x_0)=\int_{\mathbb{R}^n}\frac{u^q(y)\,{\rm d}y}{|x_0-y|^{n-\beta}}
\geq c^q \int_{\mathbb{R}^n \setminus B_{|x_0|}(0)}\frac{{\rm d}y}{|x_0-y|^{n-\beta}}
\geq c^q \int_{\mathbb{R}^n \setminus B_{|x_0|}(0)}\frac{{\rm d}y}{(2|y|)^{n-\beta}}
=\infty.
\end{equation*} It is impossible. Thus, 
$\inf_{\mathbb{R}^n}u=0$. Similarly, 
$\inf_{\mathbb{R}^n}v=0$. Therefore, 
$(u,v,w,z)$ solves Equation (1.8)’. This contradicts with the assumption of Corollary 5.2.
6. Boundary blowing-up rates
In this section, we use doubling lemma (Lemma 2.2) to prove Theorem 1.4.
Let 
$\Omega\neq\mathbb{R}^n$ be a domain of 
$\mathbb{R}^n$, and 
$u,v,w,z \in C^2(\mathbb{R}^n)$ be positive solutions of
\begin{equation}
\left\{
\begin{aligned}
&-\Delta u=wv^{p-1} & &\text{on} \ \Omega, \\
& -\Delta w=v^p & &\text{on} \ \mathbb{R}^n,\\
&-\Delta v=zu^{q-1} && \text{on} \ \Omega,\\
&-\Delta z=u^q && \text{on} \ \mathbb{R}^n.
\end{aligned}
\right.
\end{equation}Theorem 6.1. Let 
$\min\{p, q\} \gt 1.$ Assume that Equation (1.7) has no bounded positive classical solution in 
$\mathbb{R}^n.$ Then exists 
$C=C(n,p,q) \gt 0$ such that any positive classical solution 
$(u,v,w,z)$ of Equation (6.1) satisfies estimates of the boundary blow-up rates:
\begin{equation}
u(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{8p}{(2p-1)(2q-1)-1}},\quad x\in\Omega,
\end{equation}
\begin{equation}
v(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{8q}{(2p-1)(2q-1)-1}},\quad x\in\Omega,
\end{equation}and
\begin{equation}
w(x),\ z(x)\leq C[\mathrm{dist}(x, \partial\Omega)]^{-\frac{4(p+q)}{(2p-1)(2q-1)-1}},\quad x\in\Omega.
\end{equation}On the contrary, if positive classical solutions of Equation (6.1) satisfy Equations (6.2)–(6.4), then Equation (1.7) has no positive classical solution.
Proof. Write
\begin{equation*}
\sigma:=\frac{8p}{(2p-1)(2q-1)-1},\qquad \tau:=\frac{8q}{(2p-1)(2q-1)-1},
\qquad \gamma:=\frac{4(p+q)}{(2p-1)(2q-1)-1}.
\end{equation*}Then,
\begin{equation}
\sigma+2=\gamma+(p-1)\tau,\qquad \tau+2=\gamma+(q-1)\sigma, \qquad \gamma=p\tau-2=q\sigma-2.
\end{equation} Assume that one of the estimates (6.2), (6.3) and (6.4) fails. Then, there exists sequences 
$\Omega_k$, 
$(u_k, w_k, v_k, z_k)$, 
$y_k\in\Omega_k$, such that 
$(u_k, w_k, v_k, z_k)$ solves Equation (6.1) on 
$\Omega_k$ and
\begin{equation*}M_k:=u_k^{1/\sigma}+w_k^{1/\gamma}+v_k^{1/\tau}+z_k^{1/\gamma},\quad k=1,2,\ldots\end{equation*}satisfies
\begin{equation*}M_k(y_k) \gt 2k\mathrm{dist}^{-1}(y_k, \partial\Omega_k).\end{equation*} According to Lemma 2.2, there exists 
$x_k\in\Omega_k$ such that
\begin{equation*}M_k(x_k) \gt 2k\mathrm{dist}^{-1}(x_k, \partial\Omega_k),\end{equation*}
\begin{equation*}M_k(z)\leq 2M_k(x_k),\quad |z-x_k|\leq kM_k^{-1}(x_k).\end{equation*} Now we rescale 
$(u_k, w_k, v_k, z_k)$ by setting
\begin{equation*}
\widetilde{u}_k(y)=\lambda_k^\sigma u_k(x_k+\lambda_k y),
\qquad \widetilde{w}_k(y)=\lambda_k^\gamma w_k(x_k+\lambda_k y),\quad |y|\leq k,
\end{equation*}
\begin{equation*}
\widetilde{v}_k(y)=\lambda_k^\tau v_k(x_k+\lambda_k y),
\qquad \widetilde{z}_k(y)=\lambda_k^\gamma z_k(x_k+\lambda_k y),\quad |y|\leq k,
\end{equation*} with 
$\lambda_k=M_k^{-1}(x_k).$
In view of Equation (6.5), 
$\widetilde{u}_k,$ 
$\widetilde{w}_k,$ 
$\widetilde{v}_k$, 
$\widetilde{z}_k$ are also solutions of system (6.1) for 
$|y|\leq k.$
Moreover,
\begin{equation}
\left[\widetilde{u}_k^{1/\sigma}+\widetilde{w}_k^{1/\gamma}
+\widetilde{v}_k^{1/\tau}+\widetilde{z}_k^{1/\gamma}\right](0)=1,
\end{equation}
\begin{equation}
\left[\widetilde{u}_k^{1/\sigma}+\widetilde{w}_k^{1/\gamma}
+\widetilde{v}_k^{1/\tau}+\widetilde{z}_k^{1/\gamma}\right](y)\leq2,\quad |y|\leq k.
\end{equation} Applying the standard L 2 elliptic estimates and the embedding theorem, we know that the 
$C_{\rm loc}^{\delta_1}(\mathbb{R}^n)$-norm of 
$(\widetilde{u}_k, \widetilde{w}_k, \widetilde{v}_k, \widetilde{z}_k)$ is uniform bounded, where δ 1 is some number in 
$(0,1)$. Therefore, by the Schauder estimates, we can find some subsequence of 
$(\widetilde{u}_k, \widetilde{w}_k, \widetilde{v}_k, \widetilde{z}_k)$ converging to a solution 
$(\widetilde{u}, \widetilde{w}, \widetilde{v},
\widetilde{z})$ of Equation (1.7) in 
$C_{\rm loc}^{2,\delta_2}(\mathbb{R}^n)$ sense, where δ 2 is some number in 
$(0,1)$. Moreover, Equation (6.6) implies that 
$(\widetilde{u}, \widetilde{w}, \widetilde{v}, \widetilde{z})$ is non-trivial, and Equation (6.7) implies that 
$(\widetilde{u},$ 
$\widetilde{w},$ 
$\widetilde{v}$, 
$\widetilde{z})$ is a bounded solution of Equation (1.7). This contradicts the assumption of Theorem 6.1.
On the contrary, if a non-negative solution (u, v) of Equation (6.1) satisfies Equations (6.2), (6.3) and (6.4). For each 
$x_0\in \mathbb{R}^n$ and 
$R \gt 0,$ we take 
$\Omega=B(x_0, R)$. Then
\begin{equation*}u(x_0)\leq C(n,p,q)R^{-\sigma},\qquad v(x_0)\leq C(n,p,q)R^{-\tau},\end{equation*}
\begin{equation*}w(x_0),\; z(x_0)\leq C(n,p,q)R^{-\gamma}.\end{equation*} Letting 
$R\rightarrow\infty,$ we have 
$u(x_0)=w(x_0)=v(x_0)=z(x_0)=0$. Since x 0 is arbitrary, we know that Equation (1.7) has no positive solution.
Thus, we complete the proof of Theorem 6.1.
Proof of Theorem 1.4
Assume that one of the estimates (1.16) and (1.17) fails. Similar to the proof of Theorem 6.1, we write
\begin{equation*}\sigma:=\frac{(2p-1)(\beta+2)+(\alpha+2)}{(2p-1)(2q-1)-1}, \qquad \tau:=\frac{(2q-1)(\alpha+2)+(\beta+2)}{(2p-1)(2q-1)-1},\end{equation*}and
\begin{equation*}M_k:=u_k^{1/\sigma}+w_k^{1/\gamma}+v_k^{1/\tau}+z_k^{1/\eta},\quad k=1,2,\ldots,\end{equation*}where 
$w_k=|x|^{\alpha-n}\times v_k^p$, 
$z_k=|x|^{\beta-n}\times u_k^q$ and 
$\gamma=p\tau-\alpha$, 
$\eta=q\sigma-\beta$. Since 
$\sigma+2+\alpha=(2p-1)\tau$ and 
$\tau+2+\beta=(2p-1)\sigma,$ we can derive that 
$(\widetilde{u}_k, \widetilde{v}_k)$ is a solution of system (1.1) for 
$|y|\leq k.$ Similar to the proof of Theorem 6.1, there holds
\begin{equation*}
\left[\widetilde{u}_k^{1/\sigma}+\widetilde{w}_k^{1/\gamma}
+\widetilde{v}_k^{1/\tau}+\widetilde{z}_k^{1/\eta}\right](0)=1,
\end{equation*}
\begin{equation*}
\left[\widetilde{u}_k^{1/\sigma}+\widetilde{w}_k^{1/\gamma}
+\widetilde{v}_k^{1/\tau}+\widetilde{z}_k^{1/\eta}\right](y)\leq2,\quad |y|\leq k.
\end{equation*} However, by the L 2 estimates, the embedding theorem and the Schauder estimates, we can also see that the C 2-limit of some subsequence of 
$(\widetilde{u}_k, \widetilde{v}_k)$ are the bounded solution of Equation (1.1). Namely, we can also get a contradiction.
On the contrary, by the same argument in the proof of Theorem 6.1, if a non-negative solution (u, v) of Equation (1.15) satisfies Equations (1.16) and (1.17), then Equation (1.1) has no positive classical solution.
Acknowledgements
The authors thank the unknown referees very much for useful suggestions. Those suggestions have greatly improved this article.
Funding Statement
This research was supported by NNSF of China (no. 11871278) and Postgraduate Research and Practice Innovation Program of Jiangsu Province (no. KYCX23-1669).
