Proof. For a function \(h\) continuous on each sphere about the origin, write \[\overline h(r)=\frac1{|\mathbb S^{n-1}|}\int_{\mathbb S^{n-1}}h(r\theta)\,\,\mathrm d\theta.\] Integration of the equations on annuli gives, for \(r>0\), \[ \begin{split} \bigl(r^{n-1}\overline u'(r)\bigr)' &=-r^{n-1+A}\overline{v^p}(r),\\ \bigl(r^{n-1}\overline v'(r)\bigr)' &=-r^{n-1+B}\overline{u^q}(r). \end{split} \tag{3}\] In particular, each flux on the left is strictly decreasing. Its extended limit as \(r\downarrow0\) exists. If that limit were nonzero or infinite, the flux would have a fixed sign and absolute value bounded below by a positive constant for all sufficiently small \(r\). Since \(\int_0^1 r^{1-n}\,\,\mathrm dr=\infty\), integration would contradict boundedness of the corresponding spherical mean at zero. Thus both limits are zero, and both fluxes are strictly negative for every \(r>0\).

If \(n=2\), then \(r\overline u'(r)\le\overline u'(1)<0\) for \(r\ge1\). Consequently \(\overline u(r)\le\overline u(1)+\overline u'(1)\log r\), contradicting positivity for large \(r\).

Continuity and \(v(0)>0\) give \(v^p\ge c>0\) near zero. Integrating the first equation in (3) over \([r/2,r]\) and using the negative flux at \(r/2\), we obtain \[-\overline u'(r) \ge c r^{1-n}\int_{r/2}^r s^{n-1+A}\,\,\mathrm ds =c_A r^{1+A},\qquad c_A>0.\] The equality remains valid when \(n+A=0\), with the corresponding constant containing \(\log2\). If \(A\le-2\), integrating this inequality towards zero contradicts continuity of \(u\) there. Interchanging the two equations proves \(B>-2\). ◻

Proof. Local boundedness of \(u,v\) gives \(g(x)\le C_{u,v}|x|^A\) and \(f(x)\le C_{u,v}|x|^B\) near zero. Since \(n+A,n+B>0\), both sources are locally integrable. Let \(\varphi\in C_c^\infty(\mathbb R^n)\) and choose \(\xi\in C_c^\infty(\mathbb R^n)\) equal to one near zero. Test the equation for \(u\) away from zero against \(\varphi(x)(1-\xi(x/\rho))\). After expanding its Laplacian, the terms involving first and second derivatives of \(\xi(x/\rho)\) have absolute values bounded respectively by \(C_{u,v,\varphi,\xi}\rho^{n-1}\) and \(C_{u,v,\varphi,\xi}\rho^{n-2}\). They tend to zero. The remaining terms pass to the limit by local integrability, proving the first distributional equation. The second follows in the same way.

Let \(\zeta_\rho\) be a nonnegative, compactly supported smooth approximate identity, and put \(u_\rho=u*\zeta_\rho\) and \(g_\rho=g*\zeta_\rho\). Then \(u_\rho>0\) and \(-\Delta u_\rho=g_\rho\ge0\) on all of \(\mathbb R^n\). Two integrations of the spherical flux formula, now about an arbitrary point \(x\), give \[u_\rho(x)-\frac1{|\mathbb S^{n-1}|}\int_{\mathbb S^{n-1}}u_\rho(x+R\theta)\,\,\mathrm d\theta =c_n\int_{|x-y|<R} \bigl(|x-y|^{2-n}-R^{2-n}\bigr)g_\rho(y)\,\,\mathrm dy.\] The spherical average is positive. Monotone convergence as \(R\to\infty\) therefore yields \(K*g_\rho\le u_\rho\). As \(\rho\downarrow0\), \(u_\rho(x)\to u(x)\) for every \(x\), while \(g_\rho\to g\) almost everywhere. Fatou’s Lemma proves \(K*g(x)\le u(x)\) at every point. The same argument proves \(K*f\le v\).

The normalization of \(K\) gives \(-\Delta K=\delta_0\) in distributions: \(K\) is harmonic off zero, \(-\int_{\partial B_r}\partial_rK=1\), and \(\int_{\partial B_r}K=O(r)\), so integration by parts outside \(B_r\) and passage to \(r=0\) gives the assertion. To apply this to \(K*g\), note that for every compactly supported smooth \(\varphi\), \[\int_{\mathbb R^n}g(y)\int_{\mathbb R^n}K(x-y)|\Delta\varphi(x)|\,\,\mathrm dx\,\,\mathrm dy =\int_{\mathbb R^n}(K*g)(x)|\Delta\varphi(x)|\,\,\mathrm dx \le\int_{\mathbb R^n}u(x)|\Delta\varphi(x)|\,\,\mathrm dx<\infty.\] Fubini’s Theorem thus gives \(-\Delta(K*g)=g\). The difference \(h=u-K*g\) is a nonnegative locally integrable harmonic distribution. Every \(h*\zeta_\rho\) is a nonnegative entire harmonic function. Such a function \(H\) is constant: for \(b=|x-z|\), the mean-value property and \(B_R(x)\subset B_{R+b}(z)\) imply \[H(x)\le\left(\frac{R+b}{R}\right)^n H(z).\] Letting \(R\to\infty\) and then interchanging \(x,z\) proves equality. Passing back to distributions shows that \(h\) itself is a nonnegative constant. The argument for \(v-K*f\) is identical.

Fix \(x_0\ne0\). The sources are locally \(C^1\) near \(x_0\), because \(u,v\) are strictly positive there and \(p,q>0\). Split \(g\) by a smooth cutoff supported near \(x_0\) and equal to one on a smaller neighborhood. For the compactly supported part, local integrability of \(K\) and \(\nabla K\) gives a continuously differentiable convolution with the asserted derivative. For the other part, \(|x-y|\) stays bounded away from zero as \(x\) varies in a sufficiently small neighborhood of \(x_0\). There both \(K(x-y)\) and \(|\nabla K(x-y)|\) are bounded by a fixed multiple of \(K(x_0-y)\). The latter is integrable against \(g\) by (5). Dominated convergence proves continuity and differentiation of this part as well. The almost-everywhere identity in (6) is therefore pointwise off zero, and (7) follows. The proof for \(f\) is the same. ◻

Proof. Set \(m_u(R)=\inf_{B_R}u\) and \(m_v(R)=\inf_{B_R}v\). Both are positive by positivity and continuity on \(\overline{B_R}\). Integrating the weights in (8) yields \[m_u(R)\ge cR^{2+A}m_v(R)^p,\qquad m_v(R)\ge cR^{2+B}m_u(R)^q.\] Combining them and dividing by \(m_u(R)^{pq}>0\) gives \[m_u(R)^{1-pq}\ge cR^{2+A+p(2+B)}.\] The exponent on the right is positive. If \(pq<1\), the left side is at most \(u(0)^{1-pq}\), and if \(pq=1\), it is one. Either case contradicts the inequality as \(R\to\infty\). ◻

Proof. Write \(M_f=\int_{B_1}f\) and \(M_g=\int_{B_1}g\). The lower bounds (8), followed by integration of \(|x|^B\) and \(|x|^A\) on \(B_1\), give \[M_f\ge c_fM_g^q,\qquad M_g\ge c_gM_f^p,\] where \(c_f,c_g>0\) are universal. In particular, \(M_f\ge c_fc_g^qM_f^{pq}\), and \(M_f>0\) and \(pq>1\) give a universal upper bound for \(M_f\). The interchanged argument bounds \(M_g\). Apply these bounds to the scaled pair (10); changing variables proves (11). The argument uses only positive powers, and requires neither \(p\ge1\) nor \(q\ge1\).

If \(c_u>0\) in (6), then \(u\ge c_u\) almost everywhere and \[\int_{B_R}f\ge c_u^q\int_{B_R}|x|^B\,\,\mathrm dx =c_u^q\frac{|\mathbb S^{n-1}|}{n+B}R^{n+B}.\] This contradicts (11) as \(R\to\infty\), since \(n+B-(n-2-\beta)=B+2+\beta=q\alpha>0\). Likewise \(c_v=0\), proving (12).

For \(x\in B_4\) and \(j\ge3\), the annulus \(2^j\le|y|<2^{j+1}\) has \(|x-y|\ge2^{j-1}\). Its contribution to the \(f\) potential is at most \[C2^{j(2-n)}\int_{B_{2^{j+1}}}f\le C2^{-j\beta},\] and its contribution to the \(g\) potential is at most \(C2^{-j\alpha}\). The two series converge because \(\alpha,\beta>0\). Thus the potentials of the sources restricted to \(\mathbb R^n\setminus B_8\) are uniformly bounded on \(B_4\). On the other hand, \[\sup_{y\in B_8}\int_{B_4}K(x-y)\,\,\mathrm dx \le\int_{B_{12}}K(z)\,\,\mathrm dz<\infty.\] Tonelli’s Theorem, the source bounds on \(B_8\), and (12) now bound \(\int_{B_4}(u+v)\). The source part of the first estimate in (13) follows directly from (11).

Finally, on \(|x-y|\ge\ell\) the kernel is at most \(c_n\ell^{2-n}\). This bounds the contribution from \(y\in B_8\) by \(C\ell^{2-n}\). The contribution from \(y\notin B_8\) has already been bounded by \(C\); since \(0<\ell\le1\) and \(n\ge3\), this proves the second estimate in (13). ◻

Proof. We first justify absolute convergence of the separate kernel-gradient integrals, before combining them. For \(r=|x|\) sufficiently small, local boundedness of \(v\) gives \(g(y)\le C_{u,v}r^A\) when \(|x-y|\le r/2\). Hence \[|X(x)|\int_{|x-y|\le r/2}|x-y|^{1-n}g(y)\,\mathrm dy \le C_{u,v}r\,r^A\int_0^{r/2}\,\mathrm ds \le C_{u,v}r^{2+A}.\] On the complementary region, \(|X(x)|\le r\le2|x-y|\), so the same expression with that region in place of the ball is at most \(C(K*g)(x)=Cu(x)\). Both bounds remain bounded as \(x\to0\). On compact subsets away from zero, the nearby source is bounded and \(|\nabla K|\) is locally integrable; the remaining part is bounded by a constant times the finite potential. Since \(X\) is supported in \(\overline B_2\) and \(f\) is locally integrable, this proves \[\int_{\mathbb R^n}f(x)|X(x)| \int_{\mathbb R^n}|\nabla K(x-y)|g(y)\,\mathrm dy\,\mathrm dx<\infty.\] Interchanging \((f,u,B)\) and \((g,v,A)\) proves the other assertion. In particular, \(fX\cdot\nabla u\) and \(gX\cdot\nabla v\) are absolutely integrable.

Since \(a=1/(q+1)\), integration by parts outside \(B_\rho\), with \(0<\rho<1\), gives \[\begin{split} -\int_{|x|>\rho} fX\cdot\nabla u\,\mathrm dx ={}&a\int_{|x|>\rho}|x|^B u^{q+1} \bigl((n+B)\eta-w\bigr)\,\mathrm dx\\ &+a\rho^{B+1}\int_{\partial B_\rho}u^{q+1}\,\mathrm dS. \end{split}\] Here \(\operatorname{div}(|x|^BX)=|x|^B((n+B)\eta-w)\), and the boundary term has the displayed positive sign because the outward normal of the punctured domain on \(\partial B_\rho\) is \(-\omega_x\). Its magnitude is \(O_{u,v}(\rho^{n+B})\), which tends to zero. The corresponding term for \(gX\cdot\nabla v\) is positive and \(O_{u,v}(\rho^{n+A})\), and also tends to zero. There is no outer boundary term because \(X\) is compactly supported. Consequently, \[-\int fX\cdot\nabla u\,\mathrm dx-\int gX\cdot\nabla v\,\mathrm dx =a(n+B)E_f+b(n+A)E_g-aW_f-bW_g.\] Finally substitute the differentiated potential identities and \(\nabla K(z)=-(n-2)zK(z)/|z|^2\). Absolute convergence permits Fubini’s Theorem and interchanging the variables in the second integral; their sum is exactly the right side of (19). ◻

Proof. Fix \(0<\delta<1/4\), and write \(\ell(z)=\delta d(z)^m\) whenever \(d(z)>0\). Call a pair \((x,y)\) close if \(|x-y|\le\ell(z)\) for at least one endpoint \(z\in\{x,y\}\) with \(d(z)>0\), and far otherwise. Since \(d\) is 1-Lipschitz, for a close pair every point \(\xi\) of the segment \([x,y]\) satisfies \[(1-\delta)d(z)\le d(\xi)\le(1+\delta)d(z).\] In particular, the entire segment is inside \(B_2\). The derivative matrix of \(X\) is \[\nabla X(x)=\eta(x)\operatorname{Id} -w(x)\omega_x\otimes\omega_x,\] with its second term zero on \(B_1\). Taylor’s formula and (18) give, on close pairs with \(x\ne y\), \[\begin{align*} |\Phi-\eta(x)+w(x)D| &\le C\delta d(z)^{N-2+m} =C\delta d(z)^{N+n+1},\\ |\eta(x)-\eta(y)| &\le C\delta d(z)^{N-1+m} =C\delta d(z)^{N+n+2}. \end{align*}\] Both are at most \(C\delta d(z)^N\). Since \(d(z)^N\le\eta(z)\le\eta(x)+\eta(y)\), the marginal identities in (16) bound the integral of their sum over close pairs by \(C\delta E\).

For a far pair, every endpoint \(z\) in \(B_2\) satisfies \(|x-y|>\ell(z)\). Using \(|\Phi|\le(|X(x)|+|X(y)|)/|x-y|\) and \(0\le D\le1\), we obtain \[|\Phi-\eta(x)+w(x)D|+|\eta(x)-\eta(y)| \le C\sum_{\substack{z\in\{x,y\}\\d(z)>0}} \left(d(z)^{N-1}+\frac{d(z)^N}{\ell(z)}\right).\] Pairs with both endpoints outside \(B_2\) contribute zero. For the terms associated with \(x\), the tail estimate in Proposition 6, (13), bounds the integral by \[\begin{split} &C\int_{B_2}f(x) \left(d(x)^{N-1}+\frac{d(x)^N}{\ell(x)}\right) \ell(x)^{2-n}\,\mathrm dx\\ &\quad\le C\int_{B_2}f(x) \left(\delta^{2-n}d(x)^{4n+7} +\delta^{1-n}d(x)^{3n+5}\right)\,\mathrm dx \le C_\delta. \end{split}\] Indeed the two exponents are \[N-1-m(n-2)=4n+7>0,\qquad N-m(n-1)=3n+5>0,\] and \(\int_{B_2}f\le C\) by (13). The terms associated with \(y\) satisfy the same bound with \(g\). All these constants are independent of the solution.

Finally, both cutoff marginals are finite, so \[|E_f-E_g| =\left|\int_{\mathbb R^{2n}}(\eta(x)-\eta(y))\,\mathrm d\pi(x,y)\right| \le\int_{\mathbb R^{2n}}|\eta(x)-\eta(y)|\,\mathrm d\pi(x,y).\] The preceding estimates bound the left side of (20) by \(C\delta E+C_\delta\). Choosing \(\delta\) sufficiently small proves the result. ◻

Proof. If at most one coefficient is positive, then \(|c_I-c_J|=c_I+c_J\), so the inequality holds pointwise before integration. Suppose now that both coefficients are positive. Both intervals are bounded and \(\lambda>0\), so all integrals in the following calculation are finite.

Integration by parts in the first variable gives \[\int_I(\mathcal K-\overline{\mathcal K}_I)\,\,\mathrm ds =-\int_I(s-m_I)\partial_s\mathcal K\,\,\mathrm ds.\] The corresponding identity in the second variable, together with \[\partial_s\mathcal K=-(n-2)\frac{h}{\lambda^2+h^2}\mathcal K, \qquad \partial_{s'}\mathcal K=(n-2)\frac{h}{\lambda^2+h^2}\mathcal K,\] yields \[ \int_{I\times J} (2\mathcal K-\overline{\mathcal K}_I-\overline{\mathcal K}_J) \,\,\mathrm ds\,\,\mathrm ds' =(n-2)\int_{I\times J} \frac{h^2-d_0h}{\lambda^2+h^2}\mathcal K\,\,\mathrm ds\,\,\mathrm ds'. \tag{22}\]

The sign of the midpoint term follows directly from interval overlap. If \(L=|I|\) and \(M=|J|\), the pushforward of \(\,\mathrm ds\,\,\mathrm ds'\) under \((s,s')\mapsto h=s-s'\) has density \[\rho(h)=F(|h-d_0|),\qquad F(r)=\max\left\{0,\min\left\{L,M,\frac{L+M}{2}-r\right\}\right\} \quad(r\ge0).\] Indeed, \(\rho(h)=|I\cap(J+h)|\), the overlap length of two intervals whose centers are separated by \(|h-d_0|\). The function \(F\) is nonincreasing. If \(d_0\ge0\) and \(h>0\), then \(|h-d_0|\le|h+d_0|\), hence \(\rho(h)\ge\rho(-h)\); for \(d_0\le0\) the inequality is reversed. Pairing positive and negative \(h\) against the odd function \(h\mathcal K(h)/(\lambda^2+h^2)\) therefore proves \[d_0\int_{I\times J}\frac{h}{\lambda^2+h^2} \mathcal K\,\,\mathrm ds\,\,\mathrm ds'\ge0.\] Consequently the left side of (22) is at most \((n-2)\int_{I\times J}\mathcal D\mathcal K\).

Multiply this last estimate by \(c=\min(c_I,c_J)\). Each individual endpoint defect is bounded above by \(\int_{I\times J}\mathcal K\), because its subtracted endpoint average is nonnegative. The remaining coefficients are \(c_I-c\) and \(c_J-c\), whose sum is \(|c_I-c_J|\). Adding these estimates and moving the finite endpoint terms to the right proves (21). ◻

Proof. Write \[T=\int w(x)D(x,y)\,\mathrm d\pi(x,y),\qquad Q=\int\bigl(\Phi(x,y)-\eta(x)+w(x)D(x,y)\bigr)\,\mathrm d\pi(x,y).\] Lemma 8 gives, for every \(\varepsilon>0\), \[ |Q|+|E_f-E_g|\le\varepsilon E+C_\varepsilon. \tag{38}\] The virial identity, Equation (19), is therefore exactly \[ \begin{gathered} c_1E_f+c_2E_g =aW_f+bW_g-(n-2)T+(n-2)Q,\\ c_1=a(n+B)-(n-2),\quad c_2=b(n+A). \end{gathered} \tag{39}\] Proposition 10 bounds the first three terms on the right from above by \(\varepsilon E+C_\varepsilon\). Consequently \[c_1E_f+c_2E_g\le(n-1)\varepsilon E+C_\varepsilon.\] Here and below the additive constant may change, but remains independent of the solution. Since \(c_1+c_2=\gamma\), \[c_1E_f+c_2E_g =\frac{\gamma}{2}E+\frac{c_1-c_2}{2}(E_f-E_g).\] Using Equation (38), we obtain \[\frac{\gamma}{2}E \le\left(n-1+\frac{|c_1-c_2|}{2}\right)\varepsilon E+C_\varepsilon.\] Choose \(\varepsilon\) so that the coefficient of \(E\) on the right is at most \(\gamma/4\). This proves the universal bound on \(E\). The first inequality in Equation (37) follows from \(\eta=1\) on \(B_1\) and the nonnegativity of the energy densities. ◻

Proof of Theorem 1. The elementary reductions in Section 2 already exclude \(n=2\), weights at most \(-2\), and \(pq\le1\). Suppose a solution exists in the remaining range. For every \(R>0\), the pair \[u_R(x)=R^\alpha u(Rx),\qquad v_R(x)=R^\beta v(Rx)\] is an admissible solution with the same equation parameters. Proposition 12 applies to every such pair with the same constant. Both energy terms have the same dilation exponent: \[\begin{align*} &\int_{B_1} \bigl(|x|^B u_R(x)^{q+1}+|x|^A v_R(x)^{p+1}\bigr)\,\mathrm dx\\ &\hspace{2em} =R^{\alpha+\beta+2-n} \int_{B_R}\bigl(|x|^B u(x)^{q+1}+|x|^A v(x)^{p+1}\bigr)\,\mathrm dx. \end{align*}\] For \(R\ge1\), the last integral is at least its fixed, strictly positive value on \(B_1\). It is finite on every bounded ball. Equation (14) and \(1-a-b>0\) give \(\alpha+\beta+2-n>0\), so the displayed expression tends to infinity as \(R\to\infty\). This contradicts the uniform bound. ◻

Proof. The strict inequality gives nonexistence by Theorem 1. Suppose instead that \[\frac{n+A}{p+1}+\frac{n+B}{q+1}\le n-2.\] Since \(n+A,n+B>n-2>0\), \[\frac1{p+1}+\frac1{q+1} <\frac1{n-2}\left(\frac{n+A}{p+1}+\frac{n+B}{q+1}\right) \le1.\] Thus \(pq>1\), because \(1-1/(p+1)-1/(q+1)=(pq-1)/((p+1)(q+1))\).

We now apply Proposition B of Phan (Phan 2012), which states the radial ground-state theorem of Bidaut-Véron and Giacomini (Bidaut-Véron and Giacomini 2010, Theorem 1.4(i)) for the same equation and solution class, with \(n\ge3\), \(A,B>-2\), and \(p,q>0\). Its radial existence region is exactly the displayed reverse weak inequality, including equality on the critical hyperbola. In the original Hamiltonian notation, the substitution is \(N=n\), \(a=A\), \(b=B\), \(\delta=p\), and \(\mu=q\), with both differential-operator exponents equal to \(2\); the preceding calculation verifies the standing superlinearity condition \(D=\delta\mu-1>0\).

The original ground states have positive radial profiles \(U,V\) on \((0,\infty)\), classical there and continuous at zero. This origin regularity gives the strict positivity at zero required here: the radial equation makes \(r^{n-1}U'(r)\) strictly decreasing. Its limit as \(r\downarrow0\) must be zero: a nonzero or infinite limit would contradict continuity of \(U\) at zero after integration against \(r^{1-n}\). Hence \(U'(r)<0\), and \(U(0)\ge U(r)>0\) for every \(r>0\). The same argument applies to \(V\). Thus \(u(x)=U(|x|)\) and \(v(x)=V(|x|)\) belong to the stated solution class and are strictly positive on all of \(\mathbb R^n\). Radial existence in this class rules out unrestricted nonexistence under the reverse weak inequality, completing the equivalence. ◻