Learning normalized solutions of quasilinear Schrödinger operators: a constrained deep-Ritz method for the L p -supercritical case

1 Introduction

Standing waves are among the most important structures in the analysis of nonlinear dispersive equations, and their construction with a prescribed mass is both physically and mathematically significant. For the nonlinear Schrödinger equation, a standing wave of the form Ψ(x, t) = e^−*iλt*u(x) conserves the L² mass ∫ℝN|u|2, a quantity that in Bose–Einstein condensation counts the number of particles and that is therefore the physically controlled datum rather than the temporal frequency λ. Fixing the mass and letting λ be determined by the problem leads to the notion of a normalized solution, whose systematic study was initiated in the orbital-stability analysis of Cazenave and Lions [1] and placed on a general variational footing through the concentration compactness principle of Lions [2] and the ground-state theory of Berestycki and Lions [3]. In the present study, we consider the quasilinear generalization governed by the p-Laplacian,

{−Δpu+V(x)|u|p−2u=λ|u|p−2u+|u|q−2u in ℝN,∫ℝN|u|pdx=mp, u∈W1,p(ℝN),(Pm)

where 1 < p < N, the mass m > 0 is prescribed, and the pair (u, λ) ∈ W^1,*p*(ℝ^*N*) × ℝ is unknown. The quasilinear flux |∇u|^*p*−2∇u models nonlinear and degenerate diffusion and reduces to the classical Schrödinger operator when p = 2.

The variational construction of normalized solutions turns on the position of the exponent q relative to the mass-critical threshold p + p²/N. In the mass-subcritical range the energy is bounded below on the mass constraint and solutions arise as global minimizers, whereas in the mass-supercritical range the energy is unbounded, and a minimax argument is required. The decisive step in the supercritical theory was taken by Jeanjean [4], who produced a bounded Palais–Smale sequence at a mountain-pass level and thereby obtained solutions of prescribed norm in the L²-supercritical case. This circle of ideas has since developed rapidly. Bartsch and de Valeriola [5] constructed multiple normalized solutions; Bartsch and Soave [6] introduced a natural-constraint approach based on the Pohozaev manifold that has become a standard tool, and extended it to strongly coupled systems in Bartsch and Soave [7]. Soave analyzed the delicate interplay of competing nonlinearities in the mass-subcritical and mass-critical regimes [8] and, in the technically demanding Sobolev-critical case, in Soave [9]. The mass-supercritical problem was revisited with a more flexible linking geometry by Jeanjean and Lu [10], while the existence and dynamical instability of the associated standing waves were clarified by Bellazzini and Jeanjean [11]. Related minimization and stability questions were treated by Shibata [12], and multiplicity in the whole space was established by Ikoma and Miyamoto [13]. More recently, Wei and Wu [14] obtained normalized solutions in the Sobolev-critical regime, and the effect of an external potential, which breaks the translation invariance that underlies most compactness arguments, was studied by Bartsch et al. [15], by Alves et al. [16], by Molle et al. [17], and by Pellacci et al. [18].

The quasilinear setting brings its own difficulties. The p-Laplacian is degenerate for p > 2 and singular for p < 2, the associated function space W^1,*p*(ℝ^*N*) is not a Hilbert space, and the Pohozaev identity requires the more delicate regularity theory recorded by Lindqvist [19] and the quasilinear identities of Guedda and Véron [20]. The variational analysis underlying all of this rests on the classical Sobolev-space theory of Adams and Fournier [21], the Gagliardo–Nirenberg interpolation inequality in the form due to Nirenberg [22], and the fibering method of Drábek and Pohozaev [23], all of which recur throughout the present work. Normalized solutions for p-Laplacian problems have begun to receive attention, for instance in the work of Zhang and Zhang [24]. It is against this background that Miyagaki et al. [25] established the existence of a normalized solution of Equation P_*m* in the L^*p*-supercritical range p + p²/N < q < p*, for a general class of nonnegative potentials V vanishing at infinity, by constructing a Palais–Smale–Pohozaev sequence, analyzing the fibering map t ↦ J(u_*t*), and studying the ground-state map m↦γ(m). Their theory delivers three qualitative conclusions: existence of a solution with λ < 0 for each sufficiently small mass, strict monotonicity of the map m↦γ(m), and the blow-up γ(m) → +∞ as m → 0⁺. These conclusions are powerful, but they are, by their nature, non-constructive: obtained by extracting convergent subsequences and invoking deformation and Lagrange-multiplier arguments, they neither exhibit a solution profile nor evaluate any of the quantities u, λ, or γ(m). What a normalized solution looks like, how the multiplier depends on the data, and whether the monotonicity and blow-up laws can be observed quantitatively remain open questions on the computational side. Closing this gap is the purpose of the present study.

The computation of unstable critical points of elliptic energies has a history of its own. The mountain pass algorithm of Choi and McKenna [26] and the local minimax method of Li and Zhou [27] locate saddle points of grid-based discretizations by a constrained descent along paths, and the Nehari manifold optimization of Chen et al. [28] turns index-one saddles into minimizers by restricting the energy to a natural constraint and retracting onto it by scaling; for mass-constrained minimization in the subcritical regime, the normalized gradient flow of Bao and Du [29] is the classical tool. These methods are tied to a mesh and, with the exception of Bao and Du [29], treat problems in which the spectral parameter is prescribed. The natural instrument for a mesh-free treatment is a neural network solver, and mesh-free deep learning methods for partial differential equations have matured considerably. The deep Ritz method of Yu [30] minimizes a variational energy over a neural ansatz; the physics-informed neural networks of Raissi et al. [31] and the deep Galerkin method of Sirignano and Spiliopoulos [32] instead penalize the strong residual; and high-dimensional problems have been addressed by the backward-stochastic approach ofHan et al [33]. Eigenvalue problems, which share with Equation P_*m* the presence of an unknown spectral parameter, have been treated by Han et al. [34]. Multiple and unstable solutions of nonlinear elliptic equations have been computed with neural networks through homotopy continuation by Huang et al. [35] and through Newton-informed operator learning by Hao et al. [36]. At the operator level, the DeepONet architecture of Lu et al. [37] and the Fourier neural operator of Li et al. [38] learn solution maps across families of data. The approximation power underpinning these methods rests on the classical density results of Hornik [39] and Pinkus [40], on the quantitative Sobolev-norm bounds of Gühring et al. [41] for deep ReLU networks and of De Ryck et al. [42] for tanh networks, and, for the convergence of the deep Ritz method itself, on the quasi-minimizer framework of Dondl et al. [43], while the treatment of the mass constraint relies on the augmented-Lagrangian and method-of-multipliers theory of Bertsekas [44] and Nocedal and Wright [45]. None of these methods, however, addresses the normalized-solution problem Equation P_*m*: a collocation solver for the strong form requires the multiplier λ as an input, whereas λ is unknown here; the degenerate flux destabilizes automatic differentiation of the strong residual; the prescribed-mass constraint has no place in a standard residual loss, which is minimized by the spurious solution u ≡ 0; and the problem is posed on the whole of ℝ^*N* in the supercritical regime, where the energy is a saddle rather than a minimum. A bespoke method is therefore required, and its design is the subject of this study.

Our central idea is to make the numerical method follow the variational proof. The supercritical energy J has no minimum on the mass sphere: along the L^*p*-preserving dilation ut(x)=tN/pu(tx) the sought solution is a maximum of the fibering map t ↦ J(u_*t*), and any scheme that minimizes J on the sphere, with or without a penalty on the mass, has nothing to converge to in that direction. We therefore first reduce the energy along the dilations, replacing J(u) by F(u)=maxt>0J(ut), whose infimum on the mass sphere is exactly the ground-state level γ(m) of the existence theory. The reduced functional is bounded below, the unstable direction has been eliminated, and the constrained problem becomes a minimization to which an augmented-Lagrangian scheme applies. The spectral parameter λ is then recovered as the Lagrange multiplier of the mass constraint of this reduced problem, and the normalized solution itself is the dilation of the computed minimizer onto the Pohozaev manifold. Concretely, this research makes the following contributions.

• C1. A derivation of the variational framework for Equation P_*m*: the weak formulation, the Lagrange-multiplier characterization of normalized solutions, the Pohozaev identity obtained from L^*p*-preserving dilations, and the fiber reduction, including existence, uniqueness, non-degeneracy and C¹ dependence of the fiber parameter t(u) under explicit hypotheses on the potential, the natural-constraint property of the reduced functional, and the exact scaling law of γ(m) in the potential-free case.
• C2. A constrained deep Ritz method for the reduced functional in which λ emerges as the augmented-Lagrangian multiplier of the mass constraint, with a stationarity theorem for the discrete problem that identifies the multiplier exactly via a Rayleigh quotient, explicit gradients obtained from the envelope theorem, a zero-background regularization of the degenerate flux, whole-space truncation with a quantified tail bound, and unbiased or deterministic quadrature.
• C3. An error analysis: a constraint qualification, a quantitative W^{1, *p*} approximation statement matched to the cited results, a decomposition of the error of the computable objective into truncation, regularization and quadrature parts, convergence of the discrete quasi-minimizers (energies converge to γ(m), and projected iterates converge strongly to a ground state under the compactness property of the existence theory, which we prove in the radial potential-free case), Lipschitz stability of the multiplier, and residual-based diagnostics whose scope is stated precisely.
• C4. A verification methodology based on a family of exact manufactured normalized solutions with closed-form potential, multiplier and energy, together with high-accuracy radial references with and without a potential, and an outlook on learning the solution map.
• C5. A numerical study confirming the theory: the manufactured construction to machine precision, the structural identities to 10⁻⁷–10⁻⁸, the proved monotonicity and blow-up of γ(m) with the exact power law γ ∝ m⁻², and the recovery of two L²-supercritical normalized solutions, one of them with a nonzero potential, with their multipliers and energies to relative accuracy 10⁻⁵ to 10⁻⁶.

The overall logical structure of this study is summarized in Figure 1. Section 2 fixes notation and collects the definitions and standing assumptions. Section 3 derives the variational framework, Section 4 constructs the method, Section 5 carries out the error analysis, Section 6 builds the verification benchmarks, Section 7 presents the outlook on learning the solution map, Section 8 reports the experiments, Section 9 compares the approach with existing methods, and Sections 10, 11 discuss and conclude. Throughout, we keep the epistemic boundary sharp: the numerical framework provides quantitative evidence for and exploration of the existence theory, and complements, but does not replace, its proofs.

2 Preliminaries

We work throughout in the real Sobolev space W^1,*p*(ℝ^*N*) over ℝ^*N* with 1 < p < N, and we write ∥·∥_*s* for the norm of L^*s*(ℝ^*N*), s ∈ [1, ∞], and 〈·, ·〉 for the Euclidean inner product on ℝ^*N* as well as for a duality pairing where indicated. Positive constants whose value is immaterial are denoted C, C₁, C₂, … and may change between occurrences.

Definition 2.1 [21]. The Sobolev space W^1,*p*(ℝ^*N*) is the completion of Cc∞(ℝN) under the norm ∥u∥Wp=∥∇u∥pp+∥u∥pp. It is a separable, reflexive Banach space, and the embedding W^1,*p*(ℝ^*N*)↪L^*s*(ℝ^*N*) is continuous for p ≤ s ≤ p*, where

p*: =NpN-p(1)

is the critical Sobolev exponent.

Definition 2.2 [19]. The p-Laplacian is the quasilinear operator

Δpu: =div(|∇u|p-2∇u),(2)

degenerate where ∇u = 0 if p > 2 and singular there if p < 2. It is the Euler–Lagrange operator of the p-Dirichlet energy u↦1p∥∇u∥pp.

Given a potential V ≥ 0 we use the potential-weighted gradient energy

AV(u): =∥∇u∥pp+∫ℝNV(x)|u|pdx,(3)

the positive part of the energy functional below, and we write M(x): = 〈∇V(x), x〉. Under the boundedness of V assumed below, AV(u)≤(1+∥V∥∞)∥u∥Wp, but AV is not equivalent to ∥·∥Wp when V vanishes at infinity: for 0≠φ∈Cc∞(ℝN) and uR(x)=R-N/pφ(x/R) one has ∥u_*R*∥_*p* = ∥φ∥_*p*, ∥∇uR∥pp=R-p∥∇φ∥pp→0 and ∫V|uR|p=∫V(Ry)|φ(y)|pdy→0 as R → ∞, so AV(uR)→0 while ∥u_*R*∥_*W* stays bounded away from zero. Boundedness in W^1,*p*(ℝ^*N*) will therefore always be obtained from a bound on ∥∇u∥_*p* together with the prescribed mass, through the identity ∥u∥Wp=∥∇u∥pp+mp on the mass sphere.

The energy functional and the mass sphere associated with Equation P_*m* are

J(u)=1pAV(u)-1q∫ℝN|u|qdx,(4)

S(m)={u∈W1,p(ℝN):∥u∥p=m}.

Definition 2.3 Jeanjean [4]). A pair (u, λ) ∈ W^1,*p*(ℝ^*N*) × ℝ is a weak, or normalized, solution of Equation P_*m* if u ∈ S(m) and

∫ℝN|∇u|p-2∇u·∇φ+∫ℝNV|u|p-2uφ

-∫ℝN|u|q-2uφ=λ∫ℝN|u|p-2uφ

for every φ ∈ W^1,*p*(ℝ^*N*); equivalently, u is a critical point of J|_*S*(*m*) with Lagrange multiplier λ.

The exponent range and the interpolation that controls it are recorded next; the two thresholds

p<p+p2N<past(5)

separate the mass-subcritical, mass-critical, and mass-supercritical regimes.

Definition 2.4 Nirenberg [22]. The Gagliardo–Nirenberg inequality states that for p < q < p^*ast* there is a constant C = C(N, p, q) > 0 with

∥u∥qq≤C∥u∥p(1-θq)q∥∇u∥pθqq, θq: =N(q-p)pq∈(0,1),(6)

θqq=N(q-p)p.

Definition 2.5 Guedda and Veron [20]. The Pohozaev functional and the Pohozaev manifold associated with Equation P_*m* are

P(u)=∥∇u∥pp-1p∫ℝNM(x)|u|pdx-θq∥u∥qq,(7)

Pm={u∈S(m):P(u)=0},

and the ground-state level is γ(m)=infu∈PmJ(u). Every weak solution of Equation P_*m* satisfies P(u) = 0.

Definition 2.6 Drábek and Pohozaev [23]. For u ∈ W^1,*p*(ℝ^*N*), the fibering map is ϕ_*u*(t) = J(u_*t*), t > 0, where

ut(x): =tN/pu(tx)(8)

is the L^*p*-preserving dilation; critical points of ϕ_*u* correspond to dilations u_*t* lying on the Pohozaev manifold Pm.

Definition 2.7 [4]. A sequence (u_*n*)⊂S(m) is a Palais–Smale–Pohozaev sequence for J|_*S*(*m*) at level c if J(u_*n*) → c, P(u_*n*) → 0, and (J|S(m))′(un)→0 in the dual of W^1,*p*(ℝ^*N*). Compactness of such sequences is the crux of the existence theory.

We now state the standing hypotheses on the potential. The first two coincide with the hypotheses (V₁) and (V₂) of the existence theory [25]; the third is a differential condition on M under which the fiber structure of Section 3 is unique and nondegenerate, and it is the condition that our own proofs use.

Assumption 2.8.. *The potential V satisfies the following*.

• (V*₁) V ∈ C⁰(ℝ^*N*)∩L^∞(ℝ^*N*), V(x) ≥ 0 for all x ∈ ℝ^*N*, and lim|x|→∞V(x)=0*.
• (V*₂) The map x ↦ M(x) = 〈∇V(x), x〉 is well defined, belongs to C⁰(ℝ^*N*)∩L^∞(ℝ^*N*), satisfies M(x) ≤ 0 for all x, and lim|x|→∞M(x) = 0*.
• (V*₃) The map x ↦ W(x): = 〈∇M(x), x〉 is well defined, belongs to C⁰(ℝ^*N*)∩L^∞(ℝ^*N*), and*

pM(x)+W(x)≤0 for all x∈ℝN.(9)

Remark 2.9. Condition (V*₃) is checked directly on examples. For V(x) = (1 + |x|²)^−α with α > 0, writing s = |x|², one finds M = −2αs(1 + s)^−α−1, W = −4αs(1 + s)^−α−2(1 − αs) and*

pM(x)+W(x)=-2αs(1+s)-α-2[(p+2)+(p-2α)s],

so that (V₁)–(V*₃) hold precisely when 0 < α ≤ p/2; the potential (a + b|x|²)⁻¹, a, b > 0, is the case α = 1 up to a dilation and satisfies (V₃) when p ≥ 2. The Gaussian V(x) = e^−β|x|² satisfies (V₁)–(V₂), but pM + W = −2βse^−β*s*[(p+2)−2βs] changes sign, so it violates (V₃) for large |x|. The existence theory obtains uniqueness of the fiber parameter from a scaling condition on M (its hypothesis (V₅)) rather than from Equation 9; the two conditions serve the same purpose and we use Equation 9 because it is what Proposition 3.7 below requires*.

Several statements of the existence theory are used in this study as targets and as hypotheses of convergence results. They require more than (V₁)–(V₃), namely the remaining conditions of Miyagaki et al. [25], which concern integrability and scaling monotonicity of V and M, and a smallness condition on the mass. We collect what is imported in a single hypothesis, so that every result below states exactly what it uses.

Hypothesis 2.10. *The potential V satisfies hypotheses (V₁)–(V₆) of Miyagaki et al. [25], the exponent satisfies (Equation 10) below, and 0 < m < m₀, where m₀ = m₀(N, p, q, V) > 0 is the mass threshold of the principal theorem of Miyagaki et al. [25]. Under these conditions, the existence theory provides: (a) a normalized solution (u, λ) of Equation P_*m* with λ < 0, u∈Pm and J(u) = γ(m); (b) the strict monotonicity of m↦γ(m) on (0, m₀) and the blow-up γ(m) → +∞ as m → 0⁺. In the potential-free case V ≡ 0 these conclusions hold for every m > 0: (b) by the scaling argument of Proposition 3.11 below, and (a) because the scaling u ↦ au(a^{(*q*−*p*)/*p*}·), a > 0, maps normalized solutions to normalized solutions of every prescribed mass and multiplies λ by a^*q*−*p* (Section 8.3 uses this scaling for p = 2)*.

Remark 2.11. *We fix once and for all the L^*p*-supercritical exponent range*

p+p2N<q<p*=NpN-p.(10)

In terms of 2.4, the two thresholds in Equation 5 correspond to θ_qq* = p (mass-critical) and θ_*q*q = p* (Sobolev-critical), so the supercritical condition q > p + p²/N is equivalent to θ_*q*q > p, that is to N(q − p)/p > p. This inequality drives the supercritical geometry: it makes the exponent N(q − p)/p of the nonlinear term in the fibering map exceed the exponent p of the gradient term, so that J|_*S*(*m*) is unbounded and the relevant critical level is of mountain-pass type on the natural constraint Pm*.

3 The variational framework

This section derives, from first principles, the variational objects on which both the existence theory and our numerical method rest. Nothing is assumed here; every identity is obtained by an explicit computation.

3.1 Energy functional and weak formulation

We first show that Equation P_*m* is the Euler–Lagrange eqnarray of the energy J constrained to the sphere S(m), and that λ is the associated Lagrange multiplier.

Proposition 3.1. Let

J(u)=1pAV(u)-1q∥u∥qq, Ψ(u)=1p∥u∥pp,(11)

*with AV(u) as in Equation 3, and let S(m) = {u ∈ W^1,*p*(ℝ^*N*):Ψ(u) = m^*p*/p}. Then J, Ψ ∈ C¹(W^1,*p*(ℝ^*N*);ℝ), and for u ∈ W^1,*p*(ℝ^*N*) and φ ∈ W^1,*p*(ℝ^*N*)*,

〈J′(u),φ〉=∫ℝN|∇u|p-2∇u·∇φdx(12)

+∫ℝNV(x)|u|p-2uφdx-∫ℝN|u|q-2uφdx,

〈Ψ′(u),φ〉=∫ℝN|u|p-2uφdx.(13)

*Moreover, u ∈ S(m) is a critical point of J|_*S*(*m*) if and only if there exists λ ∈ ℝ such that*

〈J′(u),φ〉=λ〈Ψ′(u),φ〉 for allφ∈W1,p(ℝN),(14)

*which is precisely the weak formulation of Equation P_*m**.

Proof: Step 1: differentiability and the derivative of the gradient term. Fix u, φ ∈ W^1,*p*(ℝ^*N*) and consider g(s): =1p∥∇(u+sφ)∥pp=1p∫ℝN|∇u+s∇φ|pdx for s ∈ ℝ. The map ξ↦1p|ξ|p on ℝ^*N* is C¹ with gradient |ξ|^*p*−2ξ. Hence, pointwise,

dds1p|∇u+s∇φ|p=|∇u+s∇φ|p-2(∇u+s∇φ)·∇φ.

For |s| ≤ 1 the right-hand side is dominated in L¹(ℝ^*N*) by C(|∇u|^*p*−1+|∇φ|^*p*−1)|∇φ|, which is integrable by Hölder's inequality since |∇u|^*p*−1 ∈ L^{*p*/(*p*−1)} and |∇φ| ∈ L^*p*. By dominated convergence, we may differentiate under the integral and set s = 0:

g′(0)=∫ℝN|∇u|p-2∇u·∇φdx.

Step 2: the potential and nonlinear terms. An identical argument applied to s↦1p∫V|u+sφ|p and to s↦1q∫|u+sφ|q, using that r↦1p|r|p and r↦1q|r|q are C¹ on ℝ with derivatives |r|^*p*−2r and |r|^*q*−2r, together with V ∈ L^∞ and the Sobolev embedding W^1,*p*(ℝ^*N*)↪L^*q*(ℝ^*N*) valid for p ≤ q ≤ p*, yields

dds|s=01p∫V|u+sφ|p=∫V|u|p-2uφ,

dds|s=01q∫|u+sφ|q=∫|u|q-2uφ.

Collecting these three limits gives the Gateaux derivative (Equation 12); continuity of u ↦ J′(u) in the operator norm follows from the continuity of the Nemytskii operators involved, so J ∈ C¹. The computation of Ψ′ in Equation 13 is the special case with V ≡ 0, q = p, and the factor 1/p.

Step 3: the multiplier rule. Since 〈Ψ′(u),u〉=∥u∥pp=mp≠0 for u ∈ S(m), the constraint functional Ψ has a surjective derivative at every point of S(m), so S(m) is a C¹ Banach sub-manifold of W^1,*p*(ℝ^*N*) of codimension one and Ψ′(u) ≠ 0. By the Lagrange multiplier theorem in Banach spaces, a point u ∈ S(m) is critical for J|_*S*(*m*) if and only if J′(u) annihilates the tangent space TuS(m)=kerΨ′(u), equivalently if and only if J′(u) is proportional to Ψ′(u), i.e., Equation 14 holds for some λ ∈ ℝ. Substituting Equations 12–13 into Equation 14 gives, for all φ ∈ W^1,*p*(ℝ^*N*),

∫ℝN|∇u|p-2∇u·∇φ+∫ℝNV|u|p-2uφ

-∫ℝN|u|q-2uφ=λ∫ℝN|u|p-2uφ,

which is the weak form of −Δpu + V|u|^*p*−2u−|u|^*q*−2u = λ|u|^*p*−2u, i.e., of Equation P_*m*.

Remark 3.2. *Testing Equation 14 with φ = u gives the identity*

AV(u)-∥u∥qq=λ∥u∥pp=λmp.(15)

*In the supercritical regime, the mountain-pass solution satisfies AV(u)<∥u∥qq, whence λ < 0; this is consistent with the standing-wave interpretation e^−*iλt*u and is reproduced by the numerical multiplier in Section 4. Equation 15 provides. Furthermore, a cheap consistency diagnostic for computed solutions*.

3.2 The Pohozaev identity via L^*p*-preserving dilations

The second-order (mountain-pass) nature of J|_*S*(*m*) in the supercritical regime forces the introduction of a further constraint, the Pohozaev manifold. We derive it as a necessary condition satisfied by every solution, using the scaling (Equation 8). The key virtue of that particular scaling is that it preserves the constraint.

Lemma 3.3. *Let u ∈ W^1,*p*(ℝ^*N*) and, for t > 0, let ut(x)=tN/pu(tx). Then*

∥ut∥pp=∥u∥pp,(16)

∥∇ut∥pp=tp∥∇u∥pp,(17)

∥ut∥qq=tN(q-p)/p∥u∥qq,(18)

∫ℝNV(x)|ut|pdx=∫ℝNV(y/t)|u(y)|pdy.(19)

*In particular u_*t* ∈ S(m) whenever u ∈ S(m), for every t > 0. Moreover, (_{u_*s*)t} = u_*st* for all s, t > 0, the map u ↦ u_*t* is a linear isometry of L^*p*(ℝ^*N*) and a bounded linear bijection of W^1,*p*(ℝ^*N*), and the L^*p* pairing is dilation invariant*,

∫ℝN|ut|p-2utφtdx=∫ℝN|u|p-2uφdx(20)

for all u,φ∈W1,p(ℝN), t>0.

Proof: All four identities follow from the substitution y = tx, dy = t^*N*dx.

For Equation 16,

∥ut∥pp=∫ℝNtN|u(tx)|pdx=∫ℝN|u(y)|pdy=∥u∥pp.

For Equation 17, the chain rule gives ∇ut(x)=tN/pt∇u(tx)=tN/p+1∇u(tx), hence

∥∇ut∥pp=∫ℝNtp(N/p+1)|∇u(tx)|pdx=∫ℝNtN+p|∇u(tx)|pdx

=tp∫ℝN|∇u(y)|pdy=tp∥∇u∥pp.

For Equation 18,

∥ut∥qq=∫ℝNtNq/p|u(tx)|qdx=tNq/pt-N∫ℝN|u(y)|qdy

=tN(q-p)/p∥u∥qq.

For Equation 19,

∫ℝNV(x)|ut(x)|pdx=∫ℝNV(x)tN|u(tx)|pdx

=∫ℝNV(y/t)|u(y)|pdy,

using x = y/t. Finally Equation 16 shows the constraint ∥u∥_*p* = m is invariant under t ↦ u_*t*. The group property (us)t(x)=tN/psN/pu(stx)=ust(x) is read off from Equation 8; linearity of u ↦ u_*t* is evident, and Equations 16–17 show that it is an isometry of L^*p* and bounded on W^1,*p*(ℝ^*N*) with inverse u ↦ u_1/*t*. For Equation 20 the same substitution gives ∫t^{*N*(*p*−1)/*p*}|u(tx)|^*p*−2u(tx)t^*N*/*p*φ(tx)dx = ∫|u(y)|^*p*−2u(y)φ(y)dy.

Proposition 3.4. *For u ∈ W^1,*p*(ℝ^*N*) define the fibering map ϕ_*u*(t): = J(u_*t*), t > 0. Then*

ϕu(t)=tpp∥∇u∥pp(21)

+1p∫ℝNV(y/t)|u|pdy-tN(q-p)/pq∥u∥qq,

and its derivative satisfies

ϕu′(t)=1tP(ut), P(w): =∥∇w∥pp(22)

-1p∫ℝNM(x)|w|pdx-N(q-p)pq∥w∥qq,

*where M(x) = 〈∇V(x), x〉. In particular ϕu′(1)=P(u), so any critical point of J|_*S*(*m*) satisfies the Pohozaev identity P(u) = 0*.

Proof: Formula Equation 21 is immediate from Equation 11, Lemma 3.3. We differentiate Equation 21 in t. The gradient term gives ddttpp∥∇u∥pp=tp-1∥∇u∥pp. The nonlinear term gives

ddt(-tN(q-p)/pq∥u∥qq)=-N(q-p)pqtN(q-p)/p-1∥u∥qq.

For the potential term we differentiate under the integral (justified by V ∈ C⁰∩L^∞ and |u|^*p* ∈ L¹): since ddtV(y/t)=∇V(y/t)·(-y/t2)=-1t〈∇V(y/t),y/t〉=-1tM(y/t), we obtain

ddt1p∫ℝNV(y/t)|u|pdy=-1pt∫ℝNM(y/t)|u(y)|pdy.

Collecting the three contributions,

ϕu′(t)=tp-1∥∇u∥pp-1pt∫ℝNM(y/t)|u|pdy

-N(q-p)pqtN(q-p)/p-1∥u∥qq.

We now show this equals 1tP(ut). Using Equation 17, ∥∇ut∥pp=tp∥∇u∥pp, so 1t∥∇ut∥pp=tp-1∥∇u∥pp. Using Equation 18, 1t·N(q-p)pq∥ut∥qq=N(q-p)pqtN(q-p)/p-1∥u∥qq. For the middle term, the substitution x = y/t gives

∫ℝNM(x)|ut(x)|pdx=∫ℝNM(x)tN|u(tx)|pdx

=∫ℝNM(y/t)|u(y)|pdy,

so 1t·1p∫M(x)|ut|p=1pt∫M(y/t)|u|p. Hence the displayed expression for ϕu′(t) is exactly 1tP(ut), proving Equation 22. Setting t = 1 gives ϕu′(1)=P(u). Finally, if u ∈ S(m) is critical for J|_*S*(*m*), then t = 1 is a critical point of the scalar map t↦ϕ_*u*(t) = J(u_*t*) because u_*t* ∈ S(m) for all t by Lemma 3.3 and t ↦ u_*t* is a curve in S(m) through u; therefore ϕu′(1)=P(u)=0.

Remark 3.5. *The functional P in Equation 22 coincides with the Pohozaev functional obtained by the classical route, i.e., by testing (Equation P_*m*) against the field x·∇u and integrating by parts. Indeed, combining the test-function identity (Equation 15) with the translation-Pohozaev relation N-pp∥∇u∥pp=λNp∥u∥pp+Nq∥u∥qq-Np∫(M(x)N+V)|u|p eliminates λ and returns P(u) = 0 with the same P. Our dilation derivation has the advantage of being constructive and of exposing the scaling exponents that drive the numerical fiber reduction and continuation in Section 4*.

3.3 The Pohozaev manifold, the fiber parameter, and the reduced functional

Definition 3.6. The Pohozaev manifold and the associated ground-state level are

Pm: ={u∈S(m):P(u)=0}, γ(m): =infu∈PmJ(u).(23)

The following proposition is the analytical core of the method. It shows that, under (V₁)–(V₂), every non-zero u can be dilated onto the Pohozaev manifold at a maximum of its fibering map, that under (V₃) this dilation is unique and non-degenerate, and that the fiber parameter depends continuously differentiably on u. For u ≠ 0 and t > 0 we introduce the auxiliary function

ψu(t): =t1-pϕu′(t)=∥∇u∥pp-θqtθqq-p∥u∥qq(24)

-1ptp∫ℝNM(y/t)|u(y)|pdy,

whose zeros are the critical points of ϕ_*u*.

Proposition 3.7. *Let q satisfy (Equation 10) and let u ∈ W^1,*p*(ℝ^*N*)\ {0}*.

• (i) Under (V₁)–(V₂), ψ_u(t*) > 0 for all sufficiently small t > 0, ψ_*u*(t) → −∞ as t → ∞, and ϕ_*u* attains its supremum at some t(u) ∈ (0, ∞), which satisfies ut(u)∈Pm when u ∈ S(m)*.
• (ii) Under (V₁)–(V₃), ψ_u is strictly decreasing on (0, ∞). Consequently ϕ_u* has exactly one critical point t(u), which is its strict global maximizer and is non-degenerate:*
• (iii) Under (V₁)–(V*₃), the map t:W^1,*p*(ℝ^*N*)\{0} → (0, ∞), u ↦ t(u), is of class C¹, and t(u_*s*) = t(u)/s for every s > 0*.
• *(iv) If V ≡ 0, then*

In particular, under (V₁)–(V₃) the ground-state level admits the minimax characterization

γ(m)=infu∈S(m)maxt>0J(ut),(27)

so that it is a mountain-pass level along the dilations.

Proof: (i) By Equation 10, θ_*q*q − p > 0. Since M ≤ 0 by (V₂), the last term of Equation 24 is non-negative, so ψu(t)≥∥∇u∥pp-θqtθqq-p∥u∥qq, which is positive for small t because ∥∇u∥_*p* > 0 for u ≠ 0. For large t, |1ptp∫M(y/t)|u|p|≤∥M∥∞∥u∥pp/(ptp)→0 while tθqq-p→∞, whence ψ_*u*(t) → −∞. Hence ϕu′=tp-1ψu is positive near 0 and negative near ∞; since ϕu(t)→1plimt→0∫V(y/t)|u|pdy=0 as t → 0⁺ by dominated convergence and (V₁), and ϕ_*u*(t) → −∞ as t → ∞ by Equation 21, the continuous function ϕ_*u* attains its supremum at an interior point t(u), where ϕu′(t(u))=0, that is P(u_*t*(*u*)) = 0 by Equation 22; if u ∈ S(m) then u_*t*(*u*) ∈ S(m) by Equation 16, so ut(u)∈Pm.

(ii) Set h(t): = t^−*p*∫M(y/t)|u(y)|^*p*dy. Under (V₃) we may differentiate under the integral sign, and since ddtM(y/t)=-1t〈∇M(y/t),y/t〉=-1tW(y/t),

h′(t)=-1tp+1∫ℝN[pM+W](y/t)|u(y)|pdy ≥0

by Equation 9. Therefore

ψu′(t)=-θq(θqq-p)tθqq-p-1∥u∥qq-1ph′(t)

≤-θq(θqq-p)tθqq-p-1∥u∥qq<0,

so ψ_*u* is strictly decreasing, and by (i) it has exactly one zero t(u), with ψ_*u* > 0 on (0, t(u)) and ψ_*u* < 0 on (t(u), ∞). Thus ϕ_*u* increases on (0, t(u)) and decreases on (t(u), ∞), so t(u) is its unique critical point and strict global maximizer. Differentiating ϕu′(t)=tp-1ψu(t) and using ψ_*u*(t(u)) = 0 gives Equation 25.

(iii) Consider Θ(u, t): = ψ_*u*(t) on (W^1,*p*(ℝ^*N*)\{0}) × (0, ∞). Each of the three terms in Equation 24 is C¹ in u (the first two by Proposition 3.1, the third because v↦∫M(y/t)|v|^*p* is C¹ on L^*p* for bounded continuous M) and C¹ in t (by the computation in (ii), which uses W ∈ C⁰∩L^∞), with partial derivatives continuous in (u, t); hence Θ ∈ C¹. Since Θ(u, t(u)) = 0 and ∂tΘ(u,t(u))=ψu′(t(u))<0, the implicit function theorem yields that u ↦ t(u) is C¹ near every point. The identity t(u_*s*) = t(u)/s follows from (_{u_*s*)t} = u_*st*: the map t↦ϕ_{u_*s*}(t) = ϕ_*u*(st) is maximized at t = t(u)/s.

(iv) For V ≡ 0, Equation 24 reads ψu(t)=∥∇u∥pp-θqtθqq-p∥u∥qq, whose zero is Equation 26.

For Equation 27, let u ∈ S(m). By (i) the maximum of ϕ_*u* is attained at some t(u) with ut(u)∈Pm, so maxt>0J(ut)=J(ut(u))≥γ(m), and taking the infimum over S(m) gives infS(m)maxtJ(ut)≥γ(m). Conversely, for w∈Pm one has ϕw′(1)=P(w)=0, so t = 1 is a critical point of ϕ_*w*, hence its unique maximizer by (ii): t(w) = 1 and maxtJ(wt)=J(w). Therefore infS(m)maxtJ(ut)≤infPmJ=γ(m).

Definition 3.8. Under (V₁)–(V₃), the reduced functional is

F(u): =maxt>0J(ut)=J(ut(u)), u∈W1,p(ℝN)\{0},(28)

and π(u*): = u_*t*(*u*) is the fiber projection of u onto the Pohozaev manifold [π(u)∈Pm whenever u ∈ S(m)]*.

Proposition 3.9. Assume (V₁)–(V₃) and Equation 10. Then:

• *(i) F is invariant under dilations, F(us)=F(u) for all s > 0, and F(u)>0 for all u ≠ 0*.
• *1. F∈C1(W1,p(ℝN)\{0}), and its derivative is the derivative of J at the projected point, tested against dilated directions*:
• *(ii) γ(m)=infu∈S(m)F(u), and a point u ∈ S(m) is a critical point of F|S(m) with multiplier Λ, that is F′(u)=ΛΨ′(u), if and only if w=π(u)∈Pm is a critical point of J|_*S*(*m*) with the same multiplier, J′(w) = ΛΨ′(w). In particular every minimizer of F on S(m) projects onto a normalized solution of Equation P_*m* at the ground-state level*.
• *(iii) If V ≡ 0, then with c0: =1p-1θqq>0 and β: = pq(1 − θ_*q*)/(θ_*q*q − p) > 0*,

Proof: (i) The orbit {u_*t*:t > 0} of u_*s* coincides with that of u because (_{u_*s*)t} = u_*st*, so the maximum of J over the orbit is the same; positivity follows from F(u)≥ϕu(t)>0 for small t > 0, since ϕu(t)=tpp∥∇u∥pp(1+o(1))+1p∫V(y/t)|u|p≥tp2p∥∇u∥pp for small t by Equation 21 and V ≥ 0.

(ii) By Proposition 3.7 (iii), u ↦ t(u) is C¹, and (u, t)↦ϕ_*u*(t) = J(u_*t*) is C¹ jointly, with ∂uϕu(t)[φ]=〈J′(ut),φt〉 because u ↦ u_*t* is bounded and linear. The chain rule gives 〈F′(u),φ〉=∂uϕu(t)[φ]|t=t(u)+ϕu′(t(u))〈t′(u),φ〉, and the second term vanishes since ϕu′(t(u))=0; this is Equation 29.

(iii) The identity γ(m)=infS(m)F is Equation 27. Let u ∈ S(m), τ = t(u) and w=π(u)=uτ∈Pm. By Equations 29, 20, for every φ ∈ W^1,*p*(ℝ^*N*),

〈F′(u)-ΛΨ′(u),φ〉=〈J′(w),φτ〉-Λ〈Ψ′(w),φτ〉

=〈J′(w)-ΛΨ′(w),φτ〉.

Since φ↦φ_τ is a bijection of W^1,*p*(ℝ^*N*) (Lemma 3.3), the left-hand side vanishes for all φ if and only if the right-hand side does, which is the claim. If u minimizes F on S(m), it is a critical point of F|S(m) (Proposition 3.1, Step 3, applied to F), so w is a critical point of J|_*S*(*m*), hence a normalized solution by Proposition 3.1, with J(w)=F(u)=γ(m).

(iv) Inserting Equation 26 into Equation 21 with V ≡ 0 and using θqt(u)θqq∥u∥qq=t(u)p∥∇u∥pp gives F(u)=c0t(u)p∥∇u∥pp, which is the first formula in Equation 30. Replacing u by au multiplies ∥∇u∥pp by a^*p* and ∥u∥qq by a^*q*, so F(au)=a[pθqq-pq]/(θqq-p)F(u)=a-βF(u).

Remark 3.10. *Item (iii) is the natural-constraint property in the form used by the numerical method: the mountain-pass level of J is the minimum of F on the sphere, and the fiber direction, along which J has a maximum, has been absorbed into F, which is constant along it by (i). The price is that the minimizers of F are never isolated: the whole orbit {u_*t*} of a minimizer consists of minimizers, and only its projection π(u)∈Pm is a solution. This is harmless for a minimization algorithm, since the orbit direction is a flat direction of both F and the mass constraint, and the projection is applied once, after the optimization*.

The following facts are the qualitative laws our method is designed to reproduce numerically. In the potential-free case they follow from Equation 30 by scaling; in general they are imported from the existence theory under 2.10.

Proposition 3.11. *Let V ≡ 0 and let q satisfy Equation 10. Then for all m, a > 0*,

γ(am)=a-βγ(m), β=pq(1-θq)θqq-p>0,(31)

and γ(m) = γ(1)m^−β > 0 is strictly decreasing with γ(m*) → +∞ as m → 0⁺. For N = 3, p = 2, q = 4 one has θ_*q* = 3/4 and β = 2*.

Proof: The map u ↦ au is a bijection of S(m) onto S(am), and Equation 30 gives F(au)=a-βF(u); taking the infimum over S(m) and using Equation 27 yields (Equation 31). Positivity of γ(1) follows from the Gagliardo–Nirenberg inequality (Equation 6), which bounds ∥u∥qq by a power of ∥∇u∥_*p* on S(1) and so bounds F from below on S(1) by a positive constant through Equation 30. The values for N = 3, p = 2, q = 4 are read off from Equations 6, 31.

Proposition 3.12. *Under 2.10, the map m↦γ(m) is strictly decreasing on (0, m₀), γ(m) → +∞ as m → 0⁺, and for every m ∈ (0, m₀) there is a normalized solution (u, λ) of Equation P_*m* with u∈Pm, J(u) = γ(m) and λ < 0*.

This is a restatement of what 2.10 imports from Miyagaki et al. [25]; no proof is given here. These are the statements Section 4 is built to verify: the computed energy along the ground-state branch should decrease in m and diverge as m → 0⁺, with the exact exponent (Equation 31) when V ≡ 0. We now turn to the method itself.

4 The constrained deep-Ritz method

We construct a mesh-free approximation of the normalized solution and its multiplier by minimizing the reduced functional F of Definition 3.8 over a neural ansatz on the mass sphere, prove that the stationary points of the discrete problem satisfy the weak form of Equation P_*m* on the dilated tangent space of the ansatz with a multiplier identified exactly, and derive the ingredients of a stable implementation: the computation of the fiber parameter and of the gradient of F, the emergence of λ through the augmented Lagrangian, a regularized degenerate flux, whole-space truncation with a quantified tail bound, and quadrature. Throughout this section and the next we assume (V₁)–(V₃) and Equation 10, so that t(u), π(u), and F are defined and C¹ on W^1,*p*(ℝ^*N*)\{0}.

4.1 Neural ansatz and the reduced constrained program

Let net_ϑ:ℝ^*N* → ℝ denote a neural network with parameters ϑ, e.g., a multilayer perceptron with a smooth activation. To encode positivity of the ground state and decay at infinity, we use the ansatz

uθ(x)=cσ+(netϑ(x))χ(x), θ=(c,ϑ)∈ℝD,(32)

χ(x)=e-a|x|2 or χ(x)=e-a〈x〉,

where σ+(s)=log(1+es) is the softplus, 〈x〉=1+|x|2, a > 0 is a fixed decay rate (or a trainable parameter), and c > 0 is an explicit output scale, so that ∂_*c*uθ = uθ/c; this scale is what makes the constraint qualification of Section 5.1 automatic and what identifies the multiplier in Theorem 4.1 below. The choice (Equation 32) guarantees uθ > 0, uθ ∈ W^1,*p*(ℝ^*N*)∩L^*q*, and, for the smooth activations used here, uθ ∈ C^∞(ℝ^*N*) with all derivatives decaying like χ, so that x·∇uθ ∈ W^1,*p*(ℝ^*N*) and the curve t↦(uθ)_*t* is C¹ into W^1,*p*(ℝ^*N*).

Motivated by Proposition 3.1 and Proposition 3.9, we do not discretize the strong form of Equation P_*m* using collocation-based physics-informed networks [31], nor do we minimize the energy J itself. The reason for the second choice is the geometry of the supercritical problem: for u ∈ S(m) every dilation u_*t* lies in S(m) by Equation 16, so the mass constraint, with or without a penalty, is inactive along the dilation orbit, on which J has a strict maximum at u_*t*(*u*) by Proposition 3.7. A descent method for J on the sphere therefore runs away from the sought solution in the fiber direction, and no minimization structure supports it. Instead, in the spirit of the deep-Ritz method [30] and of the Nehari manifold optimization of Chen et al. [28], we minimize the reduced functional, which has already absorbed the fiber direction. Define

F(θ): =F(uθ)=maxt>0J((uθ)t), g(θ): =Ψ(uθ)-mpp(33)

=1p(∫ℝN|uθ|pdx-mp),

with the same symbol F for the functional and for its composition with the ansatz. The normalized-solution problem (Equation P_*m*) is then modeled by the equality-constrained program

minθ∈ℝD F(θ) subject to g(θ)=0,(34)

and the computed normalized solution is the fiber projection of the minimizer,

wθ: =π(uθ)=(uθ)t(uθ)∈Pm.(35)

By Proposition 3.9(iii), the infimum of Equation 34 over the whole space S(m) is the ground-state level γ(m), so Equation 34 is a minimization problem bounded below by γ(m) > 0, and by Proposition 3.9(i) its objective is constant along the dilation orbit of each uθ, so the position of the raw iterate on its orbit is immaterial and is fixed only by the projection Equation 35. The normalization of g in Equation 33 is chosen so that ∇_θg is the derivative of Ψ composed with the ansatz; with the alternative constraint ∫|uθ|^*p*−m^*p* = pg every multiplier below would be divided by p.

4.1.1 Computing the fiber parameter

For V ≡ 0 the fiber parameter is the closed-form expression (Equation 26), evaluated from the two integrals ∥∇uθ∥pp and ∥uθ∥qq. For a general potential, t(uθ) is the unique zero of the strictly decreasing function ψ_*uθ* of Equation 24,

ψuθ(t)=∥∇uθ∥pp-θqtθqq-p∥uθ∥qq-1ptp∫ℝNM(y/t)|uθ(y)|pdy,

which is positive for small t and negative for large t (Proposition 3.7). A bracket is obtained by doubling and halving, and bisection then converges unconditionally with one evaluation of the integral ∫M(y/t)|uθ|^*p* per step; Newton's method may be used instead, its derivative being available from the formula for ψu′ in the proof of Proposition 3.7. The integrals of the dilated potential terms are evaluated on the same quadrature nodes as the other integrals, with V and M sampled at y/t, so that no interpolation of uθ is ever needed.

4.1.2 The gradient of the reduced functional

By Proposition 3.9(ii) the derivative of F does not involve the derivative of the fiber parameter: with t = t(uθ), w = w_θ, and the dilation identities of Lemma 3.3,

∂θiF(θ)=〈J′(w),(∂θiuθ)t〉=tp∫ℝN|∇uθ|p-2∇uθ·∇∂θiuθdy+∫ℝNV(y/t)|uθ|p-2uθ∂θiuθdy

-tθqq∫ℝN|uθ|q-2uθ∂θiuθdy,(36)

all integrals being taken over the undilated variable y with ∂_θuθ and ∂_θ∇uθ supplied by reverse-mode automatic differentiation or by an explicit backward pass. Thus one evaluation of ∇θF costs one root-find for t and one gradient of the energy density, exactly as for the unreduced energy.

4.2 Emergence of the spectral parameter

We now show that the multiplier λ of Equation P_*m* arises exactly as the Lagrange multiplier of Equation 34, and that the standard method of multipliers, introduced by Bertsekas [44] and developed in the monograph of Nocedal and Wright [45], produces it. Introduce the augmented Lagrangian

L(θ;λ,μ): =F(θ)-λg(θ)+μ2g(θ)2, μ>0.(37)

Theorem 4.1. Assume (V₁)–(V₃) and Equation 10, let the ansatz be Equation 32, and let (θ^⋆, λ^⋆) be a Karush–Kuhn–Tucker (KKT) point of Equation 34, i.e.,

∇θF(θ⋆)=λ⋆∇θg(θ⋆), g(θ⋆)=0.(38)

*Put u⋆=uθ⋆, t^⋆ = t(u^⋆), w⋆=π(u⋆)=ut⋆⋆∈Pm, let*

T⋆: =span{(∂θiu⋆)t⋆: i=1,…,D}⊂W1,p(ℝN)

be the dilated tangent space of the ansatz at θ^⋆, and let ζ⋆: =ddtwt⋆|t=1=Npw⋆+x·∇w⋆ be the dilation direction. Then:

• *w^⋆ satisfies the weak form (Equation 14) of Equation P_*m* on T^⋆⊕span{ζ^⋆} with multiplier λ^⋆*: *and P(w^⋆) = 0*.
• The multiplier is given exactly by the Rayleigh quotient
• *w^⋆ is a weak solution of Equation P_*m* if and only if the residual functional J′(w^⋆)−λ^⋆Ψ′(w^⋆), which annihilates T^⋆⊕span{ζ^⋆} by (i), vanishes on all of W^1,*p*(ℝ^*N*)*.

*Moreover, any minimizer θ^⋆ of the unconstrained problem minθL(θ;λ⋆,μ) with g(θ^⋆) = 0 is such a KKT point*.

Proof: By Equation 36, ∂θiF(θ⋆)=〈J′(w⋆),(∂θiu⋆)t⋆〉, and by Equation 13 together with the dilation invariance Equation 20 of the L^*p* pairing,

∂θig(θ⋆)=〈Ψ′(u⋆),∂θiu⋆〉=〈Ψ′(w⋆),(∂θiu⋆)t⋆〉.

Substituting into Equation 38 gives 〈J′(w⋆)-λ⋆Ψ′(w⋆),(∂θiu⋆)t⋆〉=0 for every i, i.e., Equation 39 on T^⋆. For the dilation direction, the curve t↦wt⋆ is C¹ into W^1,*p*(ℝ^*N*) (the ansatz is smooth with decaying derivatives), so by the chain rule 〈J′(w⋆),ζ⋆〉=ddtJ(wt⋆)|t=1=ϕw⋆′(1)=P(w⋆) by Equation 22, and 〈Ψ′(w⋆),ζ⋆〉=ddtΨ(wt⋆)|t=1=0 by Equation 16. Since w⋆=ut⋆⋆ with t^⋆ = t(u^⋆), Proposition 3.7 gives P(w⋆)=ϕu⋆′(t⋆)t⋆=0, which proves (Equation 39) on ζ^⋆ and completes (i).

For (ii), the scale parameter gives ∂cu⋆=u⋆/c, hence (∂cu⋆)t⋆=w⋆/c∈T⋆, and testing (Equation 39) with φ = w^⋆ yields, by Equations 12–13 and ∥w⋆∥p=∥u⋆∥p=m (feasibility and Equation 16),

AV(w⋆)-∥w⋆∥qq=〈J′(w⋆),w⋆〉=λ⋆〈Ψ′(w⋆),w⋆〉=λ⋆mp,

which is Equation 40. Statement (iii) is a restatement of Proposition 3.1: w^⋆ ∈ S(m) is a weak solution with multiplier λ^⋆ if and only if Equation 14 holds for every φ ∈ W^1,*p*(ℝ^*N*). The final assertion is the first-order optimality condition for minθL(·;λ⋆,μ): differentiating Equation 37,

∇θL(θ;λ⋆,μ)=∇θF(θ)-(λ⋆-μg(θ))∇θg(θ),

which vanishes together with g(θ^⋆) = 0 exactly when Equation 38 holds.

Theorem 4.1 is a statement about the finite-dimensional problem, and it is deliberately no more than that: the KKT conditions are D scalar equations, and they yield the weak equation tested against the D dilated tangent directions and the dilation direction, a Galerkin-type stationarity. Density of a union of ansatz classes in W^1,*p*(ℝ^*N*) does not make the tangent space at one parameter value dense, so no claim about arbitrary test functions is made here; the passage to the weak equation for every test function is the object of Theorem 5.6, which works with quasi-minimizers rather than with stationarity. What the theorem does certify is the identity of the number λ^⋆ delivered by the optimizer: at a KKT point it coincides with the Rayleigh quotient Equation 40 of the projected iterate, and if the limit of the projected iterates is a solution, it is the spectral parameter of the equation.

The identity

∇θL(θ;λ,μ)=∇θF(θ)-(λ-μg(θ))∇θg(θ)(41)

also dictates the classical method-of-multipliers update: after each inner minimization producing θ_*k*, set

λk+1=λk-μkg(θk), μk+1=min(βμk,μmax), β>1.(42)

Comparing Equation 41 at successive iterates shows that when the constraint residual g(θ_*k*) → 0, the effective multiplier λk-μkg(θk)→λ⋆, so Equation 42 converges to the KKT multiplier without μ_*k* → ∞, the decisive advantage of the augmented Lagrangian over a pure quadratic penalty, which would require μ → ∞ and cause ill-conditioning. The inner minimizations are supported by a minimum structure: on the feasible set F≥γ(m)>0, and for V ≡ 0 the homogeneity (Equation 30) gives F(uθ)≥γ(1)∥uθ∥p-β, so L(·;λ,μ) is bounded below and grows both as ∥uθ∥_*p* → 0 and, through the quadratic penalty, as ∥uθ∥_*p* → ∞. Its only flat directions near a minimizer are the dilation orbit and, for V ≡ 0, translations, along which neither F nor g varies; they do not obstruct descent and are resolved by the projection (Equation 35).

4.3 Gradient formulae and the regularized degenerate flux

Two features distinguish the present setting from a standard second-order Dirichlet energy: only first derivatives of uθ enter F, and the integrand |∇uθ|^*p* is not twice differentiable where ∇uθ = 0 when p < 2. We record the exact gradient and its regularization.

The energy density in Equation 33 depends on θ only through uθ and ∇_*x*uθ. Writing u = uθ and G=∇xuθ∈ℝN, and using ∂θ(1p|G|p)=|G|p-2G·∂θG and ∂θ(1pV|u|p)=V|u|p-2u∂θu, formula Equation 36 takes the form

∇θF(θ)=∫ℝN[tp|G|p-2G·∇θG+V(y/t)|u|p-2u∇θu(43)

-tθqq|u|q-2u∇θu]dy, t=t(uθ),

where ∇_θu and ∇_θG are supplied by reverse-mode automatic differentiation. Only the first spatial derivative G = ∇_*x*uθ is differentiated with respect to x; the potentially unstable second x-derivative of the degenerate flux, which would appear in the strong form −Δpu, never occurs. This is the central numerical reason for adopting the energy (deep-Ritz) formulation.

To control the singular/degenerate factor |G|^*p*−2 we replace |G|^*p* by the smooth surrogate

Φε(G): =(|G|2+ε2)p/2-εp,(44)

∇GΦε(G)=p(|G|2+ε2)(p-2)/2G,

with ε > 0 annealed to 0 during training. The subtraction of the constant ε^*p* is essential on the whole space: without it the surrogate tends to ε^*p* > 0 as G → 0 and is not integrable over ℝ^*N*. With it, Φ_ε(0) = 0, Φ_ε ≥ 0, G↦Φ_ε(G) is C^∞, all gradients are bounded on bounded sets, and Φε→|G|p uniformly on compacta. For p = 2 the surrogate is exact, Φε(G)=|G|2.

Lemma 4.2. *Let 1 < p < ∞, ε > 0, and let Φ_ε be as in Equation 44. Then |∇GΦε(G)|≤p(|G|+ε)p-1, and the following pointwise bounds hold for all G ∈ ℝ^*N**:

1<p≤2: 0≤|G|p-Φε(G)≤min{εp,|G|p},(45)

p>2: 0≤Φε(G)-|G|p(46)

≤p2ε2(|G|2+ε2)(p-2)/2≤Cpε2(|G|p-2+εp-2),

*and, when 2 < p ≤ 4, also Φε(G)-|G|p≤p2εp-2|G|2. Consequently, if Eε denotes an energy in which |G|^*p* is replaced by Φ_ε(G) over a bounded domain Ω, then*

{εpp|Ω|,1<p≤2,Cpε2(∫Ω|G|p-2dx+εp-2|Ω|),p>2,

and over ℝ^N* the difference Φε(G)-|G|p is integrable and its integral tends to 0 as ε → 0 whenever ∫ℝN|G|rdx<∞ for all r > 0, which holds for the ansatz Equation 47. In all cases Eε→E and ∇θEε→∇θE pointwise in θ as ε → 0*.

Proof: The gradient bound follows from (|G|² + ε²)^1/2 ≤ |G|+ε. Write s = |G|² and f(s): = (s + ε²)^*p*/2−ε^*p*−s^*p*/2, so that f(0) = 0 and f′(s)=p2[(s+ε2)(p-2)/2-s(p-2)/2]. For p ≤ 2 the exponent (p−2)/2 is non-positive, so f′ ≤ 0 and f ≤ 0, i.e., Φε≤|G|p; moreover −f(s) = s^*p*/2 + ε^*p*−(s + ε²)^*p*/2 ≤ ε^*p* because (s + ε²)^*p*/2 ≥ s^*p*/2, and −f(s) ≤ s^*p*/2 because (s + ε²)^*p*/2 ≥ ε^*p*; this is Equation 45. For p > 2 the exponent is positive, so f′ ≥ 0 and f ≥ 0. To bound f, consider instead F(η): = (s+η)^*p*/2−η^*p*/2 for η ∈ [0, ε²], so that f(s) = F(ε²)−F(0) and F′(η)=p2[(s+η)(p-2)/2-η(p-2)/2]≤p2(s+ε2)(p-2)/2; integrating over [0, ε²] gives the first bound in Equation 46, and the second follows from (a+b)(p-2)/2≤Cp(a(p-2)/2+b(p-2)/2). When 2 < p ≤ 4, the mean value theorem applied to σ↦σ^(*p*−2)/2 gives (s+η)(p-2)/2-η(p-2)/2≤p-22sη(p-4)/2, whose integral over η ∈ [0, ε²] is sε^*p*−2; hence f(s)≤p2εp-2s. Integrating the pointwise bounds over Ω gives Equation 47. On ℝ^*N*, for p ≤ 2 the bound min{ε^*p*, |G|^*p*} ≤ |G|^*p* ∈ L¹ tends to zero pointwise, and dominated convergence applies; for 2 < p ≤ 4 we use p2ε2|G|p-2·2(p-2)/2 on {|G| ≥ ε} and p2εp-2|G|2≤p2ε2|G|p-2 on {|G| < ε}, so the integrand is bounded by Cε²|G|^*p*−2 with |G|^*p*−2 integrable; for p > 4 one uses instead, on {|G| < ε}, the monotonicity of F′ in η to get f(s)≤p2ε2[(s+ε2)(p-2)/2-εp-2]≤Cε2(s(p-2)/2+sεp-4), which is integrable when |G| ∈ L^*p*−2∩L². Convergence of the energies and, differentiating (Equation 44) and using the gradient bound, of the parameter gradients follow from dominated convergence.

Remark 4.3 (Choosing and annealing ε). *The regularization error (Equation 47) is deterministic, monotone in ε and, by Theorem 5.5 below, adds to a quadrature error of order n^−1/2 for n Monte-Carlo samples (or of order h⁴ for the radial trapezoidal rule of Section 4.5) and to the approximation error of the ansatz. Three practical rules follow. First, there is no benefit in driving ε below the level at which εmin(p,2)|ΩR| falls under the quadrature error, so the annealing floor should be εmin~(n-1/2/|ΩR|)1/min(p,2) for Monte-Carlo quadrature. Second, for p < 2 the Hessian of Φ_ε has norm of order ε^*p*−2 near G = 0, so the Lipschitz constant of the gradient of the objective grows like ε^*p*−2 as ε↓0; a fixed-step optimizer must reduce its step in proportion, and annealing faster than the optimizer resolves the profile only produces stiffness. A geometric schedule ε_*k*+1 = ρε_*k* with ρ ∈ [1/2, 4/5] per outer iteration of Algorithm 1, started at ε₀ comparable to the typical magnitude of |∇uθ|, and combined with a quasi-Newton inner solver with line search, is a safe default. Third, for p > 2 the regularization only smooths the density near critical points of uθ and may be kept at a fixed small value, while for p = 2 it is not needed at all. In the experiments of Section 8, where p = 2, no regularization is used*.

Algorithm 1 Constrained deep-Ritz solver for the normalized solution of Equation P_*m*. Algorithm pseudocode shown with initialization of network parameters, penalty, and multiplier, followed by an iterative optimization loop. The loop includes quadrature node selection, minimization of a penalized objective, multiplier and penalty updates, and diagnostic checks. The process stops based on convergence criteria, then returns the projected solution, multiplier, and energy. Equations referenced in blue include 26, 42, 43, and 40.

4.4 Whole-space truncation and a tail estimate

Problem (Equation P_*m*) is posed on ℝ^*N*. We solve on a ball Ω_*R* = {x:|x| < R} and control the discarded tail. The decay envelope χ in Equation 32 makes the ansatz integrable on ℝ^*N*; the truncation error is then governed by the tail of χ.

Proposition 4.4. *Let uθ be of the form (Equation 32) with |net_ϑ| ≤ B on ℝ^*N* and let s ∈ {p, q}. Then*

0≤∫ℝN|uθ|sdx-∫ΩR|uθ|sdx(48)

=∫|x|≥R|uθ|sdx≤(cσ+(B))sωN-1∫R∞rN-1χrad(r)sdr={O(RN-2e-asR2),χ=e-a|x|2,O(RN-1e-asR),χ=e-a〈x〉,

where ω_N−1 is the surface measure of the unit sphere and χ_rad(r*) denotes χ as a function of r = |x|. An identical bound holds for ∫V|uθ|^*p* (with V ∈ L^∞) and, using |∇uθ| ≤ C(1 + |x|)cσ₊(B)χ(x), for ∫|∇uθ|^*p*. Consequently all energy and constraint integrals computed on Ω_*R* differ from their ℝ^*N* values by*

τχ(R): ={RNe-amin(p,q)R2,χ=e-a|x|2,RNe-amin(p,q)R,χ=e-a〈x〉,(49)

*up to a constant factor, a quantity which is driven below any prescribed tolerance by a moderate R*.

Proof: The equality is the definition of the tail. Bounding |uθ| ≤ cσ₊(B)χ(x) (since σ₊ is increasing and |net_ϑ| ≤ B) and passing to polar coordinates gives the integral over r ≥ R; the asymptotics ∫R∞rN-1e-asr2dr~12asRN-2e-asR2 and ∫R∞rN-1e-asrdr~1asRN-1e-asR (using e^{−*as*〈*r*〉} ≤ e^−*asr*) as R → ∞ (integration by parts) yield the stated rates. The gradient bound follows from ∇uθ=cσ+′(netϑ)∇netϑχ+cσ+(netϑ)∇χ and |∇χ| = 2a|x|e^−a|x|² or |∇χ| = a|x| 〈x〉⁻¹e^{−*a*〈*x*〉} ≤ aχ, together with a bound on ∇net_ϑ, which is uniformly bounded for bounded activations such as tanh and bounded weights. For tanh networks the hypothesis |net_ϑ| ≤ B holds automatically with B the sum of the absolute values of the output weights and bias.

The bound (Equation 48) justifies choosing R so that the tail term is below the target accuracy, and verifying R-independence a posteriori by refinement.

4.5 Quadrature

All integrals over Ω_*R* are estimated by Monte-Carlo in the general, non-radial setting. Let ρ be a probability density supported on Ω_*R* (uniform, or an importance density concentrated where |uθ| is large), and draw x1,…,xn~iidρ. For any integrand f,

I^n(f): =1n∑j=1nf(xj)ρ(xj), 𝔼[I^n(f)]=∫ΩRfdx,(50)

Var[I^n(f)]=1nVarρ[fρ].

Thus I^n is unbiased and converges at the dimension-independent rate O(n^−1/2); importance sampling with ρ ∝ |uθ|^*p* minimizes the variance of the mass functional and is essential for the sharply peaked profiles that occur near the small-m threshold. Replacing each integral in F and g by its estimator (Equation 50) yields a stochastic objective whose gradient, by Equation 49 and dominated convergence, is an unbiased estimator of ∇θF at fixed t; this is precisely what stochastic optimizers (Adam, then L-BFGS on a frozen batch) consume.

When the potential is radial and the ansatz depends on x only through r = |x|, as in the experiments of Section 8, every integral reduces to ωN-1∫0Rf(r)rN-1dr and a deterministic rule is preferable. We use the composite trapezoidal rule with n equispaced nodes and spacing h = R/(n−1). By the Euler–Maclaurin formula its error is -h212[F′(R)-F′(0)]+O(h4) for F(r) = f(r)r^*N*−1; here F′(0) = 0 when N ≥ 3 and F′(R) is exponentially small, so the error is O(h⁴) up to exponentially small terms, and the quadrature term in the error decomposition of Theorem 5.5 is replaced by Ch⁴.

4.6 Continuation in the mass and branch tracing

Proposition 3.9 shows that the sought solution is the fiber projection π(u) of a minimizer of F on S(m); Proposition 3.12 predicts a smooth ground-state branch m↦(γ(m), λ(m)). We exploit both. Given a converged solution θ^⋆(m) at mass m, we initialize the solve at a nearby mass m′ by the analytic predictor obtained from the scaling: rescaling uθ↦(m′/m)uθ adjusts the mass to m′ exactly (by homogeneity of ∥·∥pp under scalar multiples), after which a short corrector solve of Equation 34 at m′ restores stationarity. Sweeping m over a grid in this predictor–corrector fashion traces the branch and yields the curves γ(m) and λ(m) whose monotonicity and blow-up are the numerical targets of Proposition 3.12; for V ≡ 0 the predictor is exact up to the fiber projection, by Equation 30.

4.7 Algorithm

Algorithm 1 realizes the theory of Section 3: the inner objective is the reduced energy F of Definition 3.8, whose minimum on the sphere is the mountain-pass level by Proposition 3.9; the returned λ is the multiplier identified by Theorem 4.1; the Rayleigh identity Equation 40 and the strong-form residual serve as correctness diagnostics that the network was never trained to satisfy; the projected output w_θ lies on Pm by construction; and the outer sweep over m produces the branch of Proposition 3.12. The method is thus a computational counterpart to the existence theory, closing the gap identified in Section 1.

4.8 Validation protocol and scope

The following independent checks are applied to a computed pair (w_θ, λ) without reference to the training objective: (i) the constraint residual |g(θ)|/m^*p*; (ii) the gap between the augmented-Lagrangian multiplier and the Rayleigh quotient Equation 40 of w_θ, which vanishes at an exact KKT point by Theorem 4.1 and therefore measures the optimization residual; (iii) the strong-form residual of Equation P_*m* at (w_θ, λ), evaluated a posteriori, which bounds the dual norm of the constrained Euler–Lagrange residual (Proposition 5.11); (iv) the Pohozaev functional P(w_θ), which vanishes by construction up to the tolerance of the fiber computation and thereby checks that computation, together with the fiber parameter t(uθ) of the raw iterate, which records the immaterial position of the iterate on its orbit; (v) agreement with a high-accuracy radial reference where one is available; and (vi) refinement independence in R, the number of quadrature nodes, the network size, and ε. We emphasize that these are consistency checks and quantitative evidence: the method computes and explores the solutions whose existence is proved by the variational theory, and does not purport to replace that proof.

5 Error analysis and convergence

Theorem 4.1 identifies the multiplier at a stationary point of the discrete program Equation 34 and gives the weak equation on the dilated tangent space of the ansatz. In this section we (i) verify the constraint qualification that guarantees a well-defined, unique multiplier; (ii) state the approximation properties of the ansatz in W^1,*p*(ℝ^*N*) as available in the literature, and derive from them the approximability of a positive solution by the ansatz Equation 32; (iii) bound the error of the computable objective, which is a truncated, regularized and quadrature-based version of F; (iv) prove that quasi-minimizers of the discrete problem have energies converging to γ(m), and that their fiber projections converge strongly to a ground state under the compactness property of the existence theory, which we verify in the radial potential-free case; (v) establish Lipschitz stability of the multiplier; and (vi) derive residual-based diagnostics and state precisely what each of them measures. Throughout, the ansatz is Equation 32 with its explicit output scale c > 0, and (V₁)–(V₃) and Equation 10 are in force.

5.1 Constraint qualification and the multiplier

Lemma 5.1. *Let uθ be given by Equation 32 with uθ ≢ 0. Then ∇_θg(θ) ≠ 0; in fact*

〈∇θg(θ),ec〉=∂cg(θ)=∫ℝN|uθ|p-2uθ∂cuθdx(51)

=1c∫ℝN|uθ|pdx=pg(θ)+mpc>0 on{g=0},

*where e_*c* is the unit vector in the c-direction. Consequently the feasible set {g = 0} is a smooth manifold near θ and the KKT multiplier λ^⋆ in Equation 8 is unique*.

Proof: By Equation 33 and Equation 13, ∂θig=∫|uθ|p-2uθ∂θiuθ. Taking θ_*i* = c and using ∂_*c*uθ = uθ/c gives ∂cg=1c∫|uθ|p=1c(pg+mp), which is strictly positive since ∫|uθ|^*p* = pg + m^*p* = m^*p* > 0 on the feasible set. Hence ∇_θg ≠ 0, i.e., LICQ holds; uniqueness of the multiplier is the standard consequence of LICQ for the stationarity system ∇θF=λ∇θg.

A second, mesh-independent characterization of the multiplier follows from testing the equation against the solution itself and is the estimator we use in practice.

Proposition 5.2. *For w ∈ S(m) define the Rayleigh multiplier*

λR(w): =AV(w)-∥w∥qqmp(52)

=1mp(∥∇w∥pp+∫ℝNV|w|pdx-∥w∥qq).

Then: (i) if (w*, λ) is a weak solution of Equation P_*m*, then λ = λ_*R*(w); (ii) if (θ^⋆, λ^⋆) is a KKT point of Equation 34, then λ⋆=λR(wθ⋆); (iii) if V ≡ 0, then λ_*R*(w) < 0 for every w∈Pm*.

Proof: (i) Take φ = w in the weak form Equation 14: the left side is 〈J′(w),w〉=AV(w)-∥w∥qq by Equation 12, and the right side is λ〈Ψ′(w),w〉=λ∥w∥pp=λmp by Equation 13. Solving for λ gives Equation 52. This is exactly the identity (Equation 15) of Equation 32. (ii) is Theorem 4.1(ii). (iii) For V ≡ 0 and w∈Pm, P(w) = 0 reads ∥∇w∥pp=θq∥w∥qq, so AV(w)-∥w∥qq=(θq-1)∥w∥qq<0 because θ_*q* < 1.

The estimator (Equation 52) depends only on integrals of the current iterate, is insensitive to the penalty parameter μ, and converges at the same rate as the underlying solution error (Proposition 5.9). We report it alongside the method-of-multipliers value (Equation 42); by (ii) their gap vanishes at an exact KKT point, so it measures the optimization residual.

5.2 Neural approximation in W^1,*p*(ℝ^*N*)

Because the constraint and energy are evaluated on the truncated ball Ω_*R* (Section 4.4), it suffices to approximate the target in W1,p(ΩR). We state the approximation results in the form in which are proved in the literature, separating the qualitative density statement, which is all that the convergence theory below uses, from the quantitative rates, which are available for specific architectures.

Theorem 5.3. Let Ω_R⊂ℝ^N* be the ball of radius R, 1 ≤ p < ∞, and let Q⊃ΩR¯ be a cube*.

• *(a) (Density.) Let the activation σ ∈ C^∞(ℝ) be non-polynomial. For every w⋆∈W1,p(ΩR) and every δ > 0 there is a single-hidden-layer network net_ϑ with activation σ such that*
• *(b) (ReLU rates.) Let k ≥ 2, B > 0. There is a constant c = c(k, N, p, B) such that for every δ ∈ (0, 1/2) and every w^⋆ ∈ W^{*k, p*}(Q) with ∥w⋆∥Wk,p(Q)≤B there is a ReLU network Φ with at most clog2(δ-k/(k-1)) layers and at most cδ-N/(k-1)log22(δ-k/(k-1)) non-zero weights such that ∥Φ-w⋆∥W1,p(Q)≤δ. Equivalently, in terms of the number W of nonzero weights, the W^{1, *p*} error is O(W-(k-1)/N) up to logarithmic factors*.
• *(c) (tanh rates.) Let k ≥ 2 and w^⋆ ∈ W^{*k*, ∞}(Q). For every N∈ℕ there is a tanh network with two hidden layers, of widths O(N) and O(NN), such that ∥Φ-w⋆∥W1,∞(Q)≤C(k,N,w⋆)(1+logN)N-(k-1), and hence the same bound holds in W^{1, *p*}(Q) up to the factor |Q|^1/*p**.

Statement (a) is the classical density theorem: single-hidden-layer networks with a smooth non-polynomial activation are dense in C¹(K) for every compact K⊂ℝ^*N*, a result going back to Hornik [39] and surveyed and sharpened by Pinkus [40]; since C1(ΩR¯) is dense in W1,p(ΩR) for the Lipschitz domain Ω_*R* and embeds continuously into it, Equation 53 follows. Statement (b) is the theorem of Gühring et al. [41] on approximation in W^{*s, p*} norms by deep ReLU networks, stated there on the unit cube with s = 1 and transferred to Q by scaling; a target given on Ω_*R* is first extended to Q by a bounded extension operator Wk,p(ΩR)→Wk,p(Q), which exists for Lipschitz domains [21]. Statement (c) is the theorem of De Ryck et al. [42] on tanh networks with two hidden layers. The rates in (b) and (c) concern specific architectures whose depth or second width grows with the accuracy; they are not asserted for a fixed-depth network with an arbitrary smooth activation, and the convergence theory below uses only (a).

Corollary 5.4. Let (u*^⋆, λ) be a weak solution of Equation P_*m* with u^⋆ ≥ 0, u^⋆≢0, and fix R > 0, c > 0 and an envelope χ as in Equation 32. Then u⋆∈Cloc1,α(ℝN) for some α ∈ (0, 1), u^⋆ > 0 on ℝ^*N*, and the function*

v⋆: =σ+-1(u⋆cχ)=log(exp(u⋆cχ)-1)(54)

belongs to C1(ΩR¯)⊂W1,∞(ΩR). For every δ > 0 there are parameters θ = (c, ϑ) with

∥uθ-u⋆∥W1,p(ΩR)≤L∥netϑ-v⋆∥W1,p(ΩR)<δ,(55)

L=c(∥χ∥L∞(ΩR)+∥∇χ∥L∞(ΩR))(1+14∥∇v⋆∥L∞(ΩR)).

*If, u⋆∈Wk,∞(ΩR) with k ≥ 2 (as is the case for p = 2 and V ∈ C^∞ by elliptic regularity), then v⋆∈Wk,∞(ΩR) and the rates of Theorem 5.3 b–c transfer to ∥uθ-u⋆∥W1,p(ΩR) with the same constant L*.

Proof: Local C^{1, α} regularity of weak solutions of −Δpu = f(x, u) with f bounded on bounded sets of u is the theorem of DiBenedetto [46] and Tolksdorf [47]; here f(x, u) = (λ−V(x))|u|^*p*−2u+|u|^*q*−2u and u⋆∈Lloc∞ by the standard Moser iteration for subcritical q. Strict positivity of the nonnegative, nontrivial solution follows from the strong maximum principle of Vázquez [48]: u^⋆ satisfies −Δpu^⋆ + β(u^⋆) ≥ 0 with β(s)=(∥V∥∞+|λ|)sp-1, which is nondecreasing with β(0) = 0 and ∫0(sβ(s))-1/pds=∞. Consequently u^⋆/(cχ) is a C¹ function on the compact set ΩR¯ with values in a compact subinterval [δ₀, K₀] of (0, ∞), on which σ+-1 is smooth; hence v⋆∈C1(ΩR¯). Given ϑ with ∥netϑ-v⋆∥W1,p(ΩR)<δ/L, write v = net_ϑ. Since σ+′=(1+e-s)-1∈(0,1) and |σ+″|≤1/4,

∥σ+(v)-σ+(v⋆)∥Lp≤∥v-v⋆∥Lp,∥∇[σ+(v)-σ+(v⋆)]∥Lp≤∥∇v-∇v⋆∥Lp+14∥∇v⋆∥L∞∥v-v⋆∥Lp,

where the second estimate uses ∇[σ+(v)-σ+(v⋆)]=σ+′(v)(∇v-∇v⋆)+(σ+′(v)-σ+′(v⋆))∇v⋆. Multiplication by cχ is bounded on W1,p(ΩR) with norm at most c(∥χ∥_∞+∥∇χ∥_∞), and cσ+(v⋆)χ=u⋆; combining gives Equation 55. The final statement follows because σ+-1 is smooth on [δ₀, K₀], so v^⋆ inherits the W^{*k*, ∞} regularity of u^⋆/(cχ), and the Lipschitz estimate Equation 55 does not depend on how net_ϑ was obtained.

5.3 Error of the computable objective

The quantity actually optimized is neither J nor F but a truncated, regularized and quadrature-based counterpart of F. Write u = uθ, G = ∇uθ, and consider the five integrals that enter the fibering map Equation 21, the constraint and the fiber Equation 24:

I1=∫ℝN|G|pdx,(56)

I2(t)=∫ℝNV(y/t)|u|pdy, I3=∫ℝN|u|qdx,I4=∫ℝN|u|pdx, I5(t)=∫ℝNM(y/t)|u|pdy.

Their computable counterparts I^j replace ℝ^*N* by Ω_*R*, |G|^*p* by Φ_ε(G), and each integral by the estimator Equation 50 with n samples. The computable fibering map is ϕ^θ(t)=tppI^1+1pI^2(t)-tθqqqI^3, the computable reduced functional is F^ε,R,n(θ)=maxt∈Kϕ^θ(t) over a fixed compact interval K⊂(0, ∞) of admissible fiber parameters, and the computable constraint is g^=I^4/p-mp/p.

Theorem 5.5. *Let uθ be of the form (Equation 32) with |net_ϑ| ≤ B and |∇net_ϑ| ≤ B on ℝ^*N*, let K = [t_min, t_max]⊂(0, ∞), and assume ε ≤ 1. There is a constant C = C(p, q, N, B, c, a, V, K) such that, with probability at least 1 − η over the Monte-Carlo sample, for j = 1, …, 5*,

supt∈K|I^j(t)-Ij(t)|≤Cεmin(p,2)(1+RN)⊮j=1︸regularization+(57)

Cτχ(R)︸truncation+Clog(1/η)n︸quadrature,

where τ_χ(R*) is the truncation rate (Equation 49) of the envelope. If the maximizer t(uθ) of ϕ_θ lies in the interior of K and ∥uθ∥_*p* ≥ m/2, then*

≤C(εmin(p,2)(1+RN)+τχ(R)+log(1/η)n),

*where λ^R is the Rayleigh quotient (Equation 52) evaluated with the computable integrals at the computable fiber maximizer. The same bounds, with a θ-uniform constant on the parameter set {|net_ϑ| ≤ B, |∇net_ϑ| ≤ B}, hold for the parameter gradients. When the radial trapezoidal rule of Section 4.5 is used instead of Equation 50, the quadrature term is replaced by Ch⁴ and the statement is deterministic*.

Proof: Fix j and t ∈ K and insert two intermediate quantities: the exact integrand over Ω_*R*, and the regularized integrand over Ω_*R*. The triangle inequality gives

where IjR,ε is the regularized integral over Ω_*R*.

(I) Truncation. Each tail integral is bounded by Proposition 4.4; for j = 2, 5 the bounded factors V(y/t), M(y/t) only contribute ∥V∥_∞, ∥M∥_∞, uniformly in t. This gives (I) ≤ Cτ_χ(R).

(II) Regularization. By Equation 47, (II)≤εp|ΩR|/p for p ≤ 2 and (II)≤Cε2(∫ΩR|G|p-2+εp-2|ΩR|) for p > 2; since |G| ≤ C(1 + |x|)χ by the gradient bound of Proposition 4.4, ∫ℝN|G|p-2≤C, and (II) ≤ Cε^{min(*p*, 2)}(1 + R^*N*).

(III) Quadrature. With f_θ the integrand of IjR,ε, which is bounded on Ω_*R* by a constant M₀ = M₀(B, c, a, R, V) because |net_ϑ|, |∇net_ϑ| and χ are bounded and Φε(G)≤(|G|+1)p, the estimator (Equation 50) is unbiased and f_θ/ρ is bounded by M=M0supΩRρ-1. Hoeffding's inequality gives Pr{|I^j-IjR,ε|>s}≤2e-2ns2/M2, i.e., (III)≤Mlog(2/η)/(2n) with probability ≥ 1 − η. For j = 2, 5 the sample is fixed while t ranges over K; since t ↦ V(y/t) and t ↦ M(y/t) are uniformly continuous on ΩR¯×K, a union bound over a finite δ-net of K [whose cardinality enters log(1/η) only through a constant depending on K and the moduli of continuity] makes the bound uniform in t ∈ K. Collecting (I)–(III) gives Equation 57.

For Equation 58, note that ϕ_θ and ϕ^θ are affine in (I₁, I₂, I₃) with coefficients bounded on K, so supK|ϕ^θ-ϕθ| is bounded by the right-hand side of Equation 57; since the maximum of ϕ_θ over K is attained at t(uθ), the elementary inequality |maxKϕ^θ-maxKϕθ|≤supK|ϕ^θ-ϕθ| gives the bound for F^ε,R,n-F. The constraint is affine in I₄. For the Rayleigh quotient, λ_*R*(w_θ) is the value at t = t(uθ) of the function Λθ(t): =(tpI1+I2(t)-tθqqI3)/I4, which is Lipschitz in t on K and in the integrals on the set {I4≥(m/2)p}. The computable maximizer t^ satisfies ϕ^θ′(t^)=0, while ϕθ′(t(uθ))=0 and ϕθ″≤-c1<0 on a neighborhood of t(uθ) by Proposition 3.7; once the right-hand side of Equation 57 is small enough for t^ to lie in that neighborhood, |t^-t(uθ)|≤|ϕθ′(t^)|/c1=|ϕθ′(t^)-ϕ^θ′(t^)|/c1≤supK|ϕθ′-ϕ^θ′|/c1, and ϕθ′ is affine in I₁, I₃ and I₅(t) with coefficients bounded on K, so this is again bounded by the right-hand side of Equation 57. Combining these facts gives the last term of Equation 58. The gradient statement follows by applying the same three estimates to Equation 43, whose integrand is again bounded and tail-controlled; uniformity over the parameter set is standard, the integrands being equi-bounded and equi-Lipschitz in θ there. For the radial trapezoidal rule the term (III) is the deterministic bound of Section 4.5.

Balancing the three terms in Equation 57 yields the practical prescription ε~n^{−1/(2min(*p*, 2))} and R~1amin(p,q)logn (Gaussian envelope) or R~1amin(p,q)logn (exponential envelope), so that all three contributions are O(n-1/2logn) up to powers of logn; Remark 4.3 translates this balance into an annealing schedule.

5.4 Convergence of the discrete quasi-minimizers

We now pass to the limit. The natural convergence notion for a deep Ritz method is that of quasi-minimizers over a sequence of ansatz classes, in the spirit of Dondl et al. [43]: it accounts for the approximation, regularization, truncation, quadrature and optimization errors through a single tolerance, and it requires no statement about stationarity on tangent spaces. Let (Uh)h∈ℕ be a sequence of ansatz classes, Uh⊂W1,p(ℝN)\{0}, each closed under multiplication by positive scalars (as Equation 32 is, through the output scale c), and let uh∈Uh be the output of the discrete solver at stage h. Set

ũh: =m∥uh∥puh∈S(m), wh: =π(ũh)=(ũh)t(ũh)∈Pm,(59)

λh: =λR(wh).

We consider the following hypotheses.

• (D) (Approximability.) There is a ground state w⋆∈Pm, J(w^⋆) = γ(m), and elements vh∈Uh∩S(m) with ∥vh-w⋆∥W→0.
• (Q) (Quasi-minimality.) F(ũh)≤inf{F(v):v∈Uh∩S(m)}+δh with δ_*h* → 0.
• (H_*c*) (Compactness.) Every minimizing sequence for J on Pm has a subsequence converging strongly in W^1,*p*(ℝ^*N*).

Hypothesis (D) is supplied by Corollary 5.4 together with the tail bound of Proposition 4.4 for growing widths and radii; (Q) is supplied by Theorem 5.5, which shows that a τ_*h*-quasi-minimizer of the computable objective is a (τ_*h*+2η_*h*)-quasi-minimizer of F on the parameter set, where η_*h* is the right-hand side of Equation 58, and by the rescaling to the sphere, which changes F by O(|∥u_*h*∥_*p* − m|) by local Lipschitz continuity of F; (H_*c*) is the compactness property of the existence theory, on which Proposition 5.7, Remark 5.8 below comment.

Theorem 5.6. Assume (V₁)–(V₃), Equation 10, (D) and (Q). Then:

• *(i) F(ũh)=J(wh)→γ(m), so (wh)⊂Pm is a minimizing sequence for J on Pm, and it is bounded in W^1,*p*(ℝ^*N*)*.
• (ii) If (H_c*) holds, then along a subsequence w_*h*→w_∞ strongly in W^1,*p*(ℝ^*N*), where w∞∈Pm is a ground state, J(w_∞) = γ(m), and a weak solution of Equation P_*m* for every test function φ ∈ W^1,*p*(ℝ^*N*), with multiplier λ_∞ = λ_*R*(w_∞); moreover λ_*h* → λ_∞ and the constrained Euler–Lagrange residual satisfies ∥J′(wh)-λhΨ′(wh)∥(W1,p(ℝN))→0.

Proof: (i) Since w⋆∈Pm, Proposition 3.7 (ii) gives t(w^⋆) = 1 and F(w⋆)=J(w⋆)=γ(m). By Proposition 3.9 (ii), F is continuous at w^⋆, so (D) gives F(vh)→γ(m), and (Q) gives F(ũh)≤F(vh)+δh. On the other hand F(ũh)=J(wh)≥γ(m) because wh∈Pm. Hence γ(m)≤J(wh)≤F(vh)+δh→γ(m). For boundedness, use P(w_*h*) = 0 in the form θq∥wh∥qq=∥∇wh∥pp-1p∫M|wh|p≤∥∇wh∥pp+1p∥M∥∞mp together with ∫V|wh|p≥0:

J(wh)≥1p∥∇wh∥pp-1qθq(∥∇wh∥pp+1p∥M∥∞mp)

=c0∥∇wh∥pp-∥M∥∞mppqθq,

with c0=1p-1qθq>0 by Equation 10. Since J(w_*h*) is bounded, so is ∥∇w_*h*∥_*p*, and ∥wh∥Wp=∥∇wh∥pp+mp is bounded.

(ii) By (i) and (H_*c*), a subsequence converges strongly, w_*h*→w_∞ in W^1,*p*(ℝ^*N*). The set Pm is closed in W^1,*p*(ℝ^*N*) (P and Ψ are continuous) and J is continuous, so w∞∈Pm and J(w_∞) = limJ(w_*h*) = γ(m): w_∞ is a ground state. Since t(w_∞) = 1, F(w∞)=J(w∞)=γ(m)=infS(m)F, so w_∞ minimizes F on S(m) and Proposition 3.9 (iii) shows that w_∞ = π(w_∞) is a critical point of J|_*S*(*m*), i.e., a weak solution of Equation P_*m* for every φ ∈ W^1,*p*(ℝ^*N*), with a multiplier that equals λ_*R*(w_∞) by Proposition 5.2 (i). The maps w↦λ_*R*(w), w ↦ J′(w) and w↦Ψ′(w) are continuous from W^1,*p*(ℝ^*N*) to ℝ and to (W^1,*p*(ℝ^*N*))^* (Proposition 5.9 and the continuity of the Nemytskii operators in Proposition 3.1), so λ_*h* → λ_∞ and J′(wh)-λhΨ′(wh)→J′(w∞)-λ∞Ψ′(w∞)=0 in (W^1,*p*(ℝ^*N*))^*.

Proposition 5.7. *Let V ≡ 0, N ≥ 2, q as in Equation 10, and replace W^1,*p*(ℝ^*N*) by the subspace W^1,*p*rad(ℝ^*N*) of radial functions in Theorem 5.6 (all ansatz classes, S(m), Pm, and γ(m) are understood in the radial class). Then (H_*c*) holds, so Theorem 5.6 (ii) applies without further assumption. The limit w_∞ is a critical point of J|_*S*(*m*) on the whole space W^1,*p*(ℝ^*N*) by the principle of symmetric criticality [49], hence a normalized solution of Equation 23. For p = 2, w_∞ is the unique positive radial solution with λ = λ_∞ [50], and γ(m) in the radial class coincides with the ground-state level of Equation 17, because a minimizer of J on Pm may be replaced by its modulus, which is a positive solution, and positive decaying solutions are radial about some point [51]*.

Proof: Let (wh)⊂Pm∩W1,prad(ℝN) satisfy J(w_*h*) → γ(m). For V ≡ 0, P(w) = 0 reads θq∥w∥qq=∥∇w∥pp, so J(w)=c0∥∇w∥pp on Pm with c0=1p-1qθq>0; hence ∥∇wh∥pp→γ(m)/c0 and (w_*h*) is bounded in W^1,*p*(ℝ^*N*). Passing to a subsequence, w_*h*harpoonupw_∞ weakly in W^1,*p*rad(ℝ^*N*), and by the compactness of the embedding W^1,*p*rad(ℝ^*N*)↪L^*q*(ℝ^*N*) for p < q < p* and N ≥ 2 [52, 53], w_*h*→w_∞ strongly in L^*q*. Therefore ∥w∞∥qq=lim∥wh∥qq=γ(m)/(c0θq)>0, so w_∞ ≠ 0, while weak lower semicontinuity gives ∥∇w∞∥pp≤γ(m)/c0=θq∥w∞∥qq and m′: =∥w∞∥p≤m. Thus ψw∞(1)=∥∇w∞∥pp-θq∥w∞∥qq≤0, and since ψ_{w_∞} is strictly decreasing with zero τ: = t(w_∞) (Proposition 3.7), τ ≤ 1. The dilation (_{w_∞)τ} lies in Pm′ and, by Equation 17, J((w∞)τ)=c0τp∥∇w∞∥pp≤c0∥∇w∞∥pp≤γ(m). Hence γ(m′) ≤ γ(m); but γ is strictly decreasing by Proposition 3.11, whose proof applies verbatim in the radial class since the scaling u ↦ au preserves radial symmetry, so m′ ≥ m, and therefore m′ = m. Now all the inequalities are equalities: γ(m)≤J((w∞)τ)≤c0∥∇w∞∥pp≤γ(m) forces τ = 1 (otherwise τ^*p* < 1 would give a strict inequality) and ∥∇w∞∥pp=lim∥∇wh∥pp. Weak convergence together with convergence of the norms in the uniformly convex spaces L^*p*(ℝ^*N*) and L^*p*(ℝ^*N*; ℝ^*N*) yields w_*h*→w_∞ and ∇w_*h* → ∇w_∞ strongly, i.e., w_*h*→w_∞ in W^1,*p*(ℝ^*N*), and w∞∈Pm with J(w_∞) = γ(m). The remaining assertions are the cited results: a critical point of J restricted to the radial functions of S(m) is a critical point on S(m) because J and Ψ are invariant under the orthogonal group [49]; uniqueness of positive radial solutions of −Δw+κ²w = w^*q*−1 is Kwong [50]; and radial symmetry of positive decaying solutions is Gidas et al. [51].

Remark 5.8. Theorem 5.6 is a consistency statement, not an independent existence proof. It presupposes (D), which requires a ground state to exist, and (Q), which requires the optimizer to return a quasi-minimizer, that is, to have reached the basin of the minimum of the reduced functional rather than a spurious local minimizer of the non-convex parametric problem; the continuation scheme of Section 4.6 is designed to maintain (Q) along the branch, and Theorem 5.5 converts the tolerance of the computable objective into the tolerance δ_h. The compactness property (H_c*) is exactly the mechanism by which the existence theory produces a ground state: under 2.10 the minimizing sequences constructed in Miyagaki et al. [25] on Pm converge strongly. For the sequence generated by the numerical scheme, this is an assumption in general, and a theorem in the radial potential-free case (Proposition 5.7). For V ≡ 0 on the whole space, (H_*c*) can hold only up to translations, which is why the statement is formulated in the radial class in Proposition 5.7; a non-radial ansatz for V ≡ 0 would require recentering the iterates. Finally, only the projected iterates w_*h* converge: the raw iterates u_*h* may sit anywhere on their dilation orbits, as Proposition 3.9 (i) predicts and Section 8 confirms*.

5.5 Stability of the multiplier

Proposition 5.9. *Let u, ũ ∈ W^1,*p*(ℝ^*N*) lie in a bounded set B={v:∥v∥W≤Θ, ∥v∥p=m}, and let λR(v)=(AV(v)-∥v∥qq)/mp be the Rayleigh multiplier (Equation 52). Then*

so the computed multiplier converges at least as fast as the solution error.

Proof: The maps v↦AV(v) and v↦∥v∥qq are C¹ with derivatives pJ0′(v) and q|v|^*q*−2v (cf. Equation 12), which are bounded operators on B: indeed |〈AV′(v),φ〉|≤p(∥∇v∥pp-1∥∇φ∥p+∥V∥∞∥v∥pp-1∥φ∥p)≤C(Θ)∥φ∥W, and similarly |〈(∥·∥qq)′(v),φ〉|≤q∥v∥qq-1∥φ∥q≤C(Θ)∥φ∥W using the Sobolev embedding W^1,*p*(ℝ^*N*)↪L^*q*. By the mean value inequality along the segment [ũ, u] (which lies in the ball {∥v∥_*W* ≤ Θ} by convexity, and along which the same bounds hold with ∥v∥_*p* ≤ Θ in place of ∥v∥_*p* = m), |AV(u)-AV(ũ)|≤C(Θ)∥u-ũ∥W and likewise for the L^*q* term. Dividing by m^*p* gives Equation 60.

5.6 Residual-based diagnostics and their scope

The Pohozaev functional P and the constrained Euler–Lagrange residual provide computable quantities that require no knowledge of the exact solution. It is important to clarify what each measures. The condition P(u) = 0 is necessary for a normalized solution but not sufficient: a point of Pm must also be a critical point of J restricted to the natural constraint. Thus, a small value of |P(u)| does not, by itself, imply that u is close to a solution. What |P(u)| controls, under nondegeneracy, is the distance from u to Pm along its own dilation orbit.

Proposition 5.10. *Let u ∈ S(m) satisfy the non-degeneracy condition ϕu″(1)≠0, where ϕ_*u*(t) = J(u_*t*) is the fibering map of Proposition 3.4, and assume that t ↦ u_*t* is C¹ into W^1,*p*(ℝ^*N*) near t = 1 (as it is for the ansatz Equation 32). Let t^⋆ = t(u) be the projection time with ut⋆=π(u)∈Pm (Proposition 3.7). Then*

∥u-π(u)∥W≤C|t⋆-1|=C|ϕu″(1)||P(u)|+O(|P(u)|2).

*Thus, the Pohozaev residual is, to leading order, proportional to the distance from u to the Pohozaev manifold along its dilation orbit. It carries no information on the distance from π(u) to a solution; in particular, for the projected output w_θ = π(uθ) of Algorithm 1 one has P(w_θ) = 0 by construction, and the numerical value of P(w_θ) only checks the accuracy of the fiber computation*.

Proof: By Equation 22, ϕu′(t)=t-1P(ut); in particular ϕu′(1)=P(u) and ϕu′(t⋆)=0. Since ϕu∈C2 near t = 1 (differentiate Equation 21 twice, using W ∈ C⁰∩L^∞), Taylor's theorem with the Lagrange remainder gives, for some ξ between 1 and t^⋆,

0=ϕu′(t⋆)=ϕu′(1)+ϕu″(ξ)(t⋆-1)

=P(u)+ϕu″(1)(t⋆-1)+[ϕu″(ξ)-ϕu″(1)](t⋆-1).

Solving the linear part for t^⋆ − 1,

t⋆-1=-P(u)ϕu″(1)-ϕu″(ξ)-ϕu″(1)ϕu″(1)(t⋆-1),

and since ϕu″ is continuous, the last coefficient is O(|t^⋆−1|), so t⋆-1=-P(u)/ϕu″(1)+O(|t⋆-1|2); as |t^⋆−1| = O(|P(u)|) by the same relation, the remainder is O(|P(u)|²), giving the first identity in Equation 61. For the second, the curve t ↦ u_*t* is C¹ into W^1,*p*(ℝ^*N*) with ddtut|t=1=Npu+x·∇u∈W1,p(ℝN), so ∥u-ut⋆∥W≤supt∈[1,t⋆]∥ddtut∥W|t⋆-1|≤C|t⋆-1|. Combining the two displays yields (Equation 61).

A statement about the distance to a solution requires the constrained Euler–Lagrange residual and a stability property of the solution. The residual is bounded by a computable quantity.

Proposition 5.11. *Let w ∈ S(m) and λ ∈ ℝ, and suppose that the strong-form residual*

ϱ(w,λ): =-Δpw+V|w|p-2w-λ|w|p-2w-|w|q-2w(62)

*is well defined as a function in L^*p*′(ℝ^*N*), p′ = p/(p−1), with the flux |∇w|^*p*−2∇w in Wloc1,p′(ℝN) and ∫|∇w|^*p*−2∇w·∇φ = −∫div(|∇w|^*p*−2∇w)φ for all φ ∈ W^1,*p*(ℝ^*N*) [for p = 2 this holds for every smooth w with Δw ∈ L², in particular for the ansatz (Equation 32)]. Then*

∥J′(w)-λΨ′(w)∥(W1,p(ℝN))*≤∥ϱ(w,λ)∥p′.(63)

Proof: For φ ∈ W^1,*p*(ℝ^*N*), Equations 12–13 and the integration by parts give 〈J′(w)−λΨ′(w), φ〉 = ∫ϱ(w, λ)φ, and Hölder's inequality yields |∫ϱφ|≤∥ϱ∥p′∥φ∥p≤∥ϱ∥p′∥φ∥W.

Proposition 5.12. *Let p = 2, 2 < q < 2^*, V satisfy (V₁), and let (w^⋆, λ^⋆) ∈ H¹(ℝ^*N*) × ℝ be a weak solution of Equation P_*m*. Consider*

T:H1(ℝN)×ℝ→H-1(ℝN)×ℝ,

T(w,λ)=(J′(w)-λΨ′(w), Ψ(w)-m2/2),

*which is of class C¹, and assume that T′(w^⋆, λ^⋆) is an isomorphism (a non-degeneracy, or inf–sup, condition). Then there are δ > 0 and C > 0 such that for all (w, λ) with ∥w-w⋆∥H1+|λ-λ⋆|<δ*,

∥w-w⋆∥H1+|λ-λ⋆|(64)

≤C(∥J′(w)-λΨ′(w)∥H-1+|Ψ(w)-m2/2|)≤C(∥ϱ(w,λ)∥2+|Ψ(w)-m2/2|),

the second inequality holding whenever Proposition 5.11 applies.

Proof: For p = 2 and 2 < q < 2^* the functional J is C² on H¹(ℝ^*N*) (the second derivative of 1q∫|w|q is (q−1)∫|w|^*q*−2φψ, continuous in w by Hölder's inequality and the Sobolev embedding), and Ψ is quadratic; hence T ∈ C¹. Since T(w^⋆, λ^⋆) = 0 and T′(w^⋆, λ^⋆) is an isomorphism, the inverse function theorem provides a neighborhood on which T is a C¹ diffeomorphism onto a neighborhood of 0 with Lipschitz inverse; the first inequality in Equation 64 is the Lipschitz bound ∥(w, λ)−(w^⋆, λ^⋆)∥ ≤ C∥T(w, λ)−T(w^⋆, λ^⋆)∥. The second follows from Equation 63 with p′ = 2. □

Remark 5.13. *The constant C in Equation 64 is not computable from the residual alone, so Equation 64 is a qualitative a posteriori statement: a small strong-form residual together with a small constraint residual implies a small error provided the solution is non-degenerate, and the iterate is already in its neighborhood. Both provisos are verifiable in the setting of Section 8. There V ≡ 0, N = 3, q = 4 and the computation is radial; the linearized operator L = −Δ−λ^⋆−(q−1)|w^⋆|^*q*−2 at the positive radial solution has a kernel spanned by the partial derivatives ∂jw⋆ [50, 54], which are not radial, so L is invertible on radial functions; and T′(w⋆,λ⋆)(v,ν)=(Lv-νw⋆,〈w⋆,v〉L2) is then invertible on Hrad1×ℝ if and only if 〈w⋆,L-1w⋆〉L2≠0. Differentiating the equation -Δwκ+κ2wκ=wκq-1 of the branch of Section 8.3 with respect to κ gives L(∂_κw_κ) = −2κw_κ, whence 〈wκ,L-1wκ〉=-14κddκ∥wκ∥22, which is nonzero because the mass along the branch is a strictly monotone power of κ (see Equation 82). Hence the non-degeneracy condition holds in the radial class for the computed example, and the strong-form residual reported in Section 8.4 is an upper bound for the error up to a local constant. For p ≠ 2 the same argument requires the finer regularity theory of the p-Laplacian and is not pursued here*.

Propositions 5.10 to 5.12 delimit the role of the diagnostics of Section 4.8: the constraint residual and the multiplier gap measure feasibility and optimization accuracy; the strong-form residual bounds the dual norm of the Euler–Lagrange residual and, under non-degeneracy, the error, while the Pohozaev functional of the projected iterate is a check on the fiber computation rather than an error estimator.

6 Verification

A numerical method for a problem whose exact solutions are unknown must be validated against cases with known answers. We construct two such cases. First (Section 6.1), we reduce Equation P_*m* to a radial ordinary differential equation, derive the exact decay rate of its ground state, and formulate the shooting problem that provides a high-accuracy reference for radial potentials, including the potential-free case V ≡ 0. Second (Section 6.3), we construct a family of manufactured normalized solutions, explicit profiles that solve Equation P_*m* exactly for an associated potential, with a closed-form multiplier and energy, which serve as ground truth for the construction of the variational identities with a nonzero potential.

6.1 Radial reduction and the shooting problem

Lemma 6.1. *If u(x) = u(r), r = |x|, is radial, then*

Δpu=(|u′|p-2u′)′+N-1r|u′|p-2u′(65)

=1rN-1(rN-1|u′|p-2u′)′.

*Consequently Equation P_*m* for a positive radial u and a radial potential V = V(r) is equivalent to the singular two-point problem*

-1rN-1(rN-1|u′|p-2u′)′+V(r)up-1=λup-1+uq-1,(66)

u′(0)=0, u(r)→0 (r→∞),

with the mass constraint ωN-1∫0∞u(r)prN-1dr=mp, where ωN-1=|𝕊N-1|.

Proof: For radial u, ∇u = u′(r)x/r and |∇u| = |u′|, so the flux is the radial field |∇u|^*p*−2∇u = |u′|^*p*−2u′(x/r). For any radial field F(x) = f(r)x/r one has divF=f′(r)+N-1rf(r)=r-(N-1)(rN-1f)′; applying this with f = |u′|^*p*−2u′ gives (Equation 65). Substituting into Equation P_*m* and using |u|^*p*−2u = u^*p*−1, |u|^*q*−2u = u^*q*−1 for u > 0 yields (Equation 66). The mass integral is the radial form of ∫ℝN|u|p=mp.

Proposition 6.2. *Let p = 2, let V be radial with V(r) → 0 as r → ∞, and let u > 0 solve Equation 66 with λ = −κ² < 0. Then, as r → ∞*,

u(r)=Cr-(N-1)/2e-κr(1+o(1)),(67)

u′(r)=-κCr-(N-1)/2e-κr(1+o(1)),

*provided that V(r) = O(r⁻²) as r → ∞; for merely vanishing V the leading exponential rate e^−κ*r* persists, with ∫1rV entering the sub-exponential factor*.

Proof: For p = 2 Equation 68 reads -u″-N-1ru′+V(r)u=λu+uq-1. Since u → 0, the term u^*q*−1 = o(u) is negligible at leading order, giving the linear equation

u″+N-1ru′-(κ2+V(r))u=0.(68)

Substitute u(r) = r^{−(*N*−1)/2}v(r). Then u′=r-(N-1)/2(v′-N-12rv) and a direct computation gives

u″+N-1ru′=r-(N-1)/2(v″-[(N-1)(N-3)4r2]v).

Hence Equation 68 becomes v″-(κ2+V(r)+(N-1)(N-3)4r2)v=0, a perturbation of v″ = κ²v by an O(r⁻²) potential when V(r) = O(r⁻²). By the standard asymptotic theory of such equations (a regular-singular analysis at infinity, or comparison with the modified Bessel equation, whose decaying solution is Kν(κr)~π/2κre-κr with ν = (N−2)/2), the recessive solution satisfies v(r) = Ce^−κ*r*(1 + o(1)). Therefore u(r) = Cr^{−(*N*−1)/2}e^−κ*r*(1 + o(1)), and differentiating the leading term [the O(r⁻²) correction does not affect the leading exponential rate] gives the stated expression for u′. For a potential that vanishes more slowly, the Liouville–Green approximation v~(κ2+V)-1/4exp(-∫rκ2+V) shows that the exponential rate is unchanged while the prefactor acquires the factor exp(-∫rV/(2κ)+… ).

Proposition 6.2 justifies the exponential envelope in the ansatz and fixes the far-field behavior that the shooting reference must reproduce. The reference itself is generated as follows. Recast (Equation 66) with p = 2 as the first-order system in (u, w) with w = u′,

u′=w, w′=-N-1rw+V(r)u-λu-uq-1,(69)

u(0)=α>0, w(0)=0,

and, for fixed λ = −κ², view the terminal behavior as a function of the shooting amplitude α. For a supercritical exponent there is a threshold α^⋆(λ) separating solutions that cross zero in finite r (oscillatory/nodal) from those that blow up; α^⋆ is the positive, monotone ground state whose existence and uniqueness for the model semilinear problem is classical [3, 50], isolated by bisection on α using the sign of the first zero as the objective. The pair (λ, α^⋆) is then tuned in an outer loop so that the computed mass ωN-1∫0∞u2rN-1dr equals m²; monotonicity of the mass in κ makes this a well-posed scalar root-find (for V ≡ 0 the mass is an explicit power of κ, see Equation 82, and no outer loop is needed; for the potential of Section 8.2 we fix κ and take the resulting mass as the target of the neural solver). Integrating Equation 69 with a high-order adaptive Runge–Kutta scheme [regularizing the 1/r term at r = 0 via the Taylor data w(r)=1N((V(0)-λ)α-αq-1)r+O(r3) obtained from Equation 69] yields a reference solution of controllable accuracy.

6.2 Radial Pohozaev identity as an ODE-side check

Specializing P(u) = 0 of Proposition 3.4 to radial functions and using the radial forms of the three integrals,

∥∇u∥pp=ωN-1∫0∞|u′|prN-1dr,∥u∥qq=ωN-1∫0∞uqrN-1dr,∫ℝNM|u|pdx=ωN-1∫0∞rV′(r)uprN-1dr

[the last one because M(x) = rV′(r) for a radial potential], gives the scalar identity

∫0∞|u′|prN-1dr-1p∫0∞rV′(r)uprN-1dr(70)

=N(q-p)pq∫0∞uqrN-1dr,

a relation the shooting solution must satisfy independently of how it was computed; its residual measures the quality of the reference. The test-function identity Equation 15 provides a second, independent check in the form ∫(|u′|^*p* + Vu^*p*)r^*N*−1−∫u^*q*r^*N*−1 = λ∫u^*p*r^*N*−1.

6.3 Manufactured normalized solutions

To validate the variational identities with a potential, we prescribe the profile and construct the data that make it an exact solution.

Theorem 6.3. *Fix λ ∈ ℝ and a profile w ∈ C²(ℝ^*N*) with w > 0, w, ∇w → 0 at infinity, and ∇w(x) ≠ 0 for x ≠ 0. Define*

Vw(x): =λ+w(x)q-p+Δpw(x)w(x)p-1.(71)

*Then w is an exact classical solution of the equation in Equation P_*m* with potential V = V_*w* and multiplier λ, i.e.*,

-Δpw+Vwwp-1=λwp-1+wq-1 inℝN,(72)

*and w is normalized to mass m_*w*: = ∥w∥_*p*. Its energy is given in closed form by*

J(w)=(1p-1q)∥w∥qq+λpmwp.(73)

*Moreover, P(w) = 0: the manufactured solution lies on the Pohozaev manifold*.

Proof: By construction Equation 71, multiply by w^*p*−1 > 0:

Vwwp-1=λwp-1+wq-pwp-1+Δpw=λwp-1+wq-1+Δpw,

which rearranges to Equation 72. Thus w solves the PDE with V = V_*w*; positivity and the decay of w make w ∈ W^1,*p*(ℝ^*N*), and its mass is ∥w∥_*p* = m_*w* by definition. For the energy, test Equation 72 against w and integrate by parts: ∫(-Δpw)w=∥∇w∥pp, so ∥∇w∥pp+∫Vwwp=λmwp+∥w∥qq, i.e., AVw(w)=λmwp+∥w∥qq [recall AVw(w)=∥∇w∥pp+∫Vw|w|p]. Substituting into J(w)=1pAVw(w)-1q∥w∥qq gives

J(w)=1p(λmwp+∥w∥qq)-1q∥w∥qq=(1p-1q)∥w∥qq+λpmwp,

which is Equation 73. Finally P(w) = 0 is automatic: w is a weak solution, hence by Proposition 3.4 it satisfies the Pohozaev identity (the required regularity of M_*w* = 〈∇V_*w*, x〉 holds because w ∈ C² and w > 0 make Vw∈C0).

Corollary 6.4. *Let p = 2, let 〈r〉: =1+r2, and take the smooth radial profile w(x) = e^{−κ〈*r*〉} with κ > 0 and λ = −κ². Then*

Vw(x)=-κ2︸λ+e-(q-2)κ〈r〉︸wq-p→0+Δww︸→ κ2,(74)

Δww=κ2r2〈r〉2-κ〈r〉3-(N-1)κ〈r〉,

*so that Vw∈C∞(ℝN)∩L∞(ℝN) and V_*w*(x) → 0 as |x| → ∞. Hence w is an exact normalized solution of Equation P_*m* for a smooth, bounded potential that vanishes at infinity. The sign of V_*w*, however, places it outside 2.8 : one has Vw(0)=-κ2+e-(q-2)κ-Nκ<0 and Vw=-(N-1)κ〈r〉-1+O(〈r〉-2)<0 in the tail, and for the parameters used in Section 8 the function r ↦ V_*w*(r) is negative for every r and strictly increasing, so that Mw=rVw′(r)>0 for r > 0. Thus V_*w* satisfies the regularity and decay requirements of (V₁) and (V₂) but violates both sign conditions, V ≥ 0 and M ≤ 0: it is a potential well rather than a nonnegative potential. This does not diminish its role as an exact benchmark for the identities of Section 3, since Theorem 6.3 requires no sign condition on V_*w*; it does mean that the fiber structure of Proposition 3.7, and hence the reduced method, is not asserted for V_*w*, which is why the solver is tested in Section 8 against radial references with admissible potentials rather than against w. A manufactured example inside the class 2.8 would require a profile whose reconstructed potential is nonnegative with M ≤ 0, which we do not pursue here*.

Proof: For p = 2, Δp = Δ and Vw=λ+wq-2+Δw/w by Equation 71. With w = e^{−κ〈*r*〉} and 〈r〉=1+r2, we have 〈r〉′ = r/〈r〉, so w′=-κr〈r〉w and, using Δw=w″+N-1rw′ with w″=(κ2r2〈r〉2-κ〈r〉-r〈r〉′〈r〉2)w=(κ2r2〈r〉2-κ1〈r〉3)w,

Δww=κ2r2〈r〉2-κ1〈r〉3-N-1r·κr〈r〉=κ2r2〈r〉2-κ〈r〉3-(N-1)κ〈r〉.

As r → ∞, r²/〈r〉² → 1 and the remaining terms vanish, so Δw/w → κ²; combined with w^*q*−2 → 0 and λ = −κ² this gives V_*w* → 0, in agreement with Equation 74. Smoothness and boundedness of V_*w* follow from smoothness and boundedness of each term (all denominators 〈r〉 ≥ 1). The value at the origin follows from r²/〈r〉² = 0 and 〈r〉 = 1 there, and the tail expansion from r²/〈r〉² = 1 − 〈r〉⁻². The monotonicity statement for the parameters of Section 8 is verified numerically in Section 8.1.

6.4 Acceptance criteria for the benchmark

Given a network output (w_θ, λ_net) and a reference (u^⋆, λ^⋆) from the shooting solution of Section 6.1, we report the relative errors

ELp=∥wθ-u⋆∥p∥u⋆∥p, EW1,p=∥wθ-u⋆∥W∥u⋆∥W,(75)

Eλ=|λnet-λ⋆||λ⋆|, EJ=|J(wθ)-J(u⋆)||J(u⋆)|,

together with the reference-free diagnostics of Section 4.8: the constraint residual |g|/m^*p*, the gap between the augmented-Lagrangian and Rayleigh multipliers, the strong-form residual of Equation P_*m*, and the Pohozaev functional of the projected iterate. A run is accepted when Equation 75 and the residuals fall below prescribed tolerances and exhibit the refinement behavior predicted by Theorem 5.5. The manufactured family of Theorem 6.3 is used to verify the construction itself and the closed-form identities (Equations 72, 73) and P(w) = 0 to machine precision, so that the computer implementation of the operators and functionals is checked against an exact object before it is applied to a problem without closed-form solution.

7 Outlook: learning the solution map

The existence theory is stated for a class of potentials, and the per-instance solver of Section 4 treats one triple (V, m, q) at a time. It is natural to ask for the map from the data to the solution. This section offers a formal outlook: it records what such an operator-learning extension would require, states the one approximation fact available under explicit hypotheses, and formulates a training objective. Nothing in this section is validated numerically in the present study, whose experiments (Section 8) concern the per-instance solver only.

7.1 Obstructions and the required ingredients

Three issues stand between the existence theory and a well-defined solution operator. First, the theory does not define a single-valued map (V, m, q)↦(u, λ): even for V ≡ 0 the translates u(·−x₀) of a solution are solutions, and in general several critical points may share the level γ(m). A canonical selection is therefore needed, for instance the radially symmetric decreasing representative when the potential is radial [for p = 2 and V ≡ 0 this selection is unique by Kwong [50]], or the representative whose center of mass is at the origin, or the class modulo translations. Second, continuity of the scalar level m↦γ(m) (Proposition 3.12) says nothing about continuity of a selected profile u_ξ in ξ = (V, m, q); continuity of the selected branch in the W^1,*p*(ℝ^*N*) topology is a separate statement, which for p = 2, V ≡ 0 follows from the explicit scaling of the branch in Section 8.3 but is not available in general. Third, an approximation theorem for operators applies to continuous single-valued maps on compact input sets with outputs on a bounded physical domain [37]; it does not by itself provide uniform approximation in W^1,*p*(ℝ^*N*) × ℝ on the whole space, and the admissible set of potentials must be made compact, for instance by restricting to a finite-dimensional family V = V_α, α ∈ A⊂ℝ^*d*, with A compact.

Under these three provisos, the following finite-rank approximation statement holds; its proof is elementary and is included only to delimit what is and what is not established.

Lemma 7.1. *Let A be a compact metric space of data, R > 0, and let G:A→W1,p(ΩR)×ℝ, ξ↦(u_ξ|_{Ω_*R*}, λ_ξ), be continuous. Then for every δ > 0 there are K ∈ ℕ, continuous functions β0,…,βK:A→ℝ and functions τ1,…,τK∈W1,p(ΩR) such that*

supξ∈A(‖uξ-∑k=1Kβk(ξ)τk‖W1,p(ΩR)+|λξ-β0(ξ)|)<δ.(76)

If A is a compact subset of ℝ^d, the β_k may be taken to be networks with a smooth non-polynomial activation, and each τ_k may be taken to be a network as in Theorem 5.3 (a).

Proof: The image G(A) is compact in the Banach space X=W1,p(ΩR)×ℝ, hence covered by finitely many balls of radius δ/2 centered at G(ξ1),…,G(ξK). Let (η_*k*) be a continuous partition of unity on A subordinate to the open cover {ξ:∥G(ξ)-G(ξk)∥X<δ/2}. Then ∑kηk(ξ)G(ξk) is within δ/2 of G(ξ) for every ξ, which is Equation 76 with β_*k* = η_*k*, τ_*k* = u_{ξ_*k*} and β0=∑kηkλξk, up to δ/2. Replacing the η_*k* by network approximations uniformly within δ/(4Kmaxk∥G(ξk)∥X) on the compact set A⊂ℝd, and each τ_*k* by a network within δ/(4K) in W1,p(ΩR) (Theorem 5.3a), keeps the total error below δ.

The form Equation 76 is the branch–trunk structure of a DeepONet [37]; a Fourier neural operator [38] provides an alternative parametrization. The lemma concerns the existence of an approximant, not about its computation, and its hypotheses (single-valued selection, continuity of the selected branch, compact data set) are precisely the items listed above, which remain to be established for Equation P_*m* beyond the potential-free radial case.

7.2 A physics-informed operator objective

Writing J_ξ and g_ξ for the energy and constraint of Section 4 at datum ξ, the natural training objective combines a data loss on solver-generated pairs with the variational identities of Section 3,

minΘ 𝔼ξ~ν[‖uΘ(ξ)-uξ‖W2︸data+ω1|gξ(uΘ(ξ))|2︸mass(77)

+ω2|Pξ(uΘ(ξ))|2︸Pohozaev+ω3Rξ(uΘ(ξ),λΘ(ξ))2︸PDE residual],

where ν is a sampling distribution over A and Rξ is the strong-form residual of Equation P_*m*. The last three terms are reference-free and could allow the operator to extrapolate to data for which no pre-computed solution is available; the data pairs (ξ, u_ξ) would be generated by Algorithm 1. Whether such a surrogate can be trained to match the per-instance solver's accuracy, and how the selection and continuity issues above manifest in practice, are open computational questions that we leave for future work. In any case, the surrogate would be checked against the variational solver and could not replace it, and it would not certify existence.

8 Numerical experiments

We now report numerical results. Every number in this section is computed; none is quoted from the literature. We work in the setting p = 2, N = 3 and q = 4, which lies in the L²-supercritical window q > 2+4/N = 10/3, with two potentials: V ≡ 0, which is covered by the existence theory for every mass by Proposition 3.11 and admits the exact scaling of Section 8.3, and V(x) = (1 + |x|²)⁻¹, which satisfies (V₁)–(V₃) by Remark 2.9 (the case α = 1 = p/2) and for which the fiber parameter has no closed form. All integrals use |𝕊²| = 4π. The computations are organized to verify, in turn, the manufactured-solution construction of Theorem 6.3, the radial references and their structural identities, the qualitative laws of Proposition 3.11, and finally the constrained deep-Ritz solver of Section 4.

8.1 Verification of manufactured solutions

Theorem 6.3 asserts that, for the profile w(x) = e^{−κ〈*r*〉} and multiplier λ = −κ², the potential Vw=λ+wq-p+Δpw/wp-1 makes w an exact solution. We verify this directly: on a fine radial grid we evaluate the analytic Δw/w of Corollary 6.4, assemble V_*w*, and measure the residual of -Δw+Vww-λw-wq-1 relative to ∥wq-1∥∞. Table 1 shows the residual is at the level of double-precision round-off across three (N, q) pairs, confirming the construction and its closed-form potential; the analytic and numerically differentiated Δw/w agree to 10⁻¹⁵.

Table 1. Verification of the manufactured-solution family (Theorem 6.3 and Corollary 6.4), κ = 0.8.

Second column: maximal discrepancy between the closed-form Δw/w of Equation 74 and its numerical evaluation; third column: relative residual of the Equation 72 with the reconstructed potential; last column: value of V_*w* at r = 8, negative in all cases. The PDE residual is at machine precision, confirming that w solves Equation P_*m* exactly with the constructed potential V_*w*.

Figure 2 displays the manufactured solution and its constructed potential for (N, q) = (3, 4); the potential is smooth, bounded, and decays to zero. As stated in Corollary 6.4, it is a potential well: for the three pairs of Table 1 we find numerically that V_*w*(r) < 0 for all r ∈ [0, 40], with V_*w*(0) = −2.838, −2.149 and −3.638 respectively, and that V_*w* is strictly increasing, so that Mw=rVw′>0 for r > 0; V_*w* thus violates the sign conditions of (V₁) and (V₂) while satisfying their regularity and decay requirements. This does not affect its use as an exact benchmark for the construction, which is all that this subsection uses it for.

8.2 Radial references and structural identities

For V ≡ 0 we compute the positive ground state of the fixed-frequency form -U″-N-1rU′+κ2U=Uq-1 (with λ = −κ²) by the shooting method of Section 6.1. Bisection on the amplitude α = U(0) across the separatrix between the divergent and nodal regimes yields, for κ = 1,

α⋆=U(0)=4.3373876800,(78)

after which the profile is integrated with an adaptive high-order scheme and truncated at the radius r_cut = 8.11 where U falls below 10⁻⁴ (beyond which the growing mode of the linearization would contaminate the tail; cf. Proposition 6.2); the neglected tail contributes less than 10⁻⁸ in relative terms to every integral below. The resulting integrals are

A=∥∇U∥22=56.6917, B=∥U∥qq=75.5890,(79)

D=∥U∥22=18.8972,

so the mass at κ = 1 is m=D=4.347096. The reference passes the two structural identities of the theory to the accuracy of the quadrature: the radial Pohozaev identity (Equation 70) and the test-function identity (Equation 15) hold with relative residuals

These identities are independent of how the profile was computed, so their near-vanishing certifies the reference. The profile U(r) appears (solid) in Figure 3. The exponential decay predicted by Proposition 6.2, namely U(r)~Cr^{−(*N*−1)/2}e^−κ*r*, is confirmed in Figure 4, where the computed ground state is indistinguishable from the theoretical envelope Cr⁻¹e^−*r* over four decades.

For the potential V(r) = (1 + r²)⁻¹ the same shooting procedure applied to Equation 69, with κ = 1 fixed, gives the amplitude α^⋆ = U(0) = 5.4310130317, the truncation radius r_cut = 7.73, and the integrals

A=70.9575, B=98.3279, D=16.3424,(81)

IV=∫ℝ3VU2dx=11.0280,IM=∫ℝ3MU2dx=-5.5769,

so that the reference has mass m=D=4.042578, multiplier λ = −1 and energy γ=12A+12IV-14B=16.410753. The Pohozaev identity with potential, A-12IM-34B=0, and the test-function identity A+IV+κ2D-B=0 hold with relative residuals 6.8 × 10⁻⁸ and 8.7 × 10⁻⁸. The fibering map ϕU(t)=t22A+12∫V(y/t)U2-t34B, evaluated from the reference integrals on a grid of 2001 values of t ∈ [0.5, 1.5], attains its maximum at t = 1.0000 with value γ, as Proposition 3.7 requires (Figure 5). This reference is used in Section 8.4 as the target for the neural solver with a potential.

8.3 The ground-state branch

We now verify the qualitative laws of Proposition 3.11. Using the exact dilation covariance uκ(x)=κ2/(q-2)U(κx), which solves Equation P_*m* for V ≡ 0 with λ = −κ², we find that the mass and energy along the ground-state branch are pure powers of κ,

m(κ)=Dκ12(4q-2-N), γ(κ)=12Aκ4q-2+2-N(82)

-1qBκ2qq-2-N,

where the two energy exponents coincide (both equal 4q-2+2-N) because uκ∈Pm. For N = 3, q = 4 the exponent of κ in m is 12(4q-2-N)=-12, so m decreases monotonically in κ, while the energy exponent is +1 > 0, so γ → +∞ as κ → ∞, i.e. as m → 0⁺. Eliminating κ gives the exact power law

γ(m)∝m2(4q-2+2-N)/(4q-2-N)=m-2,(83)

in agreement with the exponent β = 2 of Proposition 3.11, and a least-squares fit to the computed branch over κ ∈ [0.25, 4] recovers the slope −2.0000 to four digits. Table 2 lists representative points and Figure 6 displays the branch: γ(m) is strictly decreasing and blows up as m → 0⁺, and λ(m) = −κ² increases to 0⁻ as m grows. This is a direct numerical confirmation of the monotonicity and blow-up statements of Proposition 3.11.

Table 2. Computed ground-state branch for N = 3, q = 4, V ≡ 0 (m=Dκ-1/2, λ = −κ²).

Columns: dilation parameter κ of the branch u_κ = κU(κx); its mass m; its energy γ(m) = J(u_κ); its multiplier λ(m); and the product γm², which is constant to the displayed precision, confirming the exact law Equation 83.

8.4 The constrained deep-Ritz solver

Finally, we test the method of Section 4. The reduced functional Equation 28 is minimized on the mass sphere; for p = 2, V ≡ 0, N = 3, q = 4, the fiber maximum has the closed form

F(u)=maxt>0J(ut)=827Au3Bu2,(84)

Au=∥∇u∥22, Bu=∥u∥qq, t(u)=4Au3Bu,

obtained by substituting the maximizer of ϕu(t)=t22Au-t34Bu, in agreement with Equation 26 and Equation 30 (here θ_*q* = 3/4, c₀ = 1/6). For the potential V = (1 + r²)⁻¹ the fiber parameter is computed by bisection on ψu(t)=Au-34tBu-12t2∫M(y/t)u2, which takes 51 steps to reach relative accuracy 10⁻¹⁵, and F(u)=t22Au+12∫V(y/t)u2-t34Bu at t = t(u). The fiber structure is illustrated in Figure 5: for both references the map ϕ_*U*(t) attains its maximum at t(U) = 1 with value ϕ_*U*(1) = γ(m), so that U lies on the Pohozaev manifold, as it must.

8.4.1 Setup

We parametrize uθ(r)=cσ+(∑k=1Hvktanh(wkr+bk)+b0)e-〈r〉, a single-hidden-layer network with the explicit output scale c of Equation 32 (3H+2 parameters in all, c initialized to 1) and the exponential envelope of Equation 32 with a = 1, which matches the decay rate e^−κ*r* of Proposition 6.2. All integrals are evaluated by the composite trapezoidal rule of Section 4.5 on [0, R] with R = 12 and 4, 000 nodes, so that the truncation error of Proposition 4.4 is below e^−2*R*R³≈10⁻⁷ and the quadrature error is of order h⁴≈10⁻¹⁰; no flux regularization is needed for p = 2. The constraint is g(θ)=12(∥uθ∥22-m2), as in Equation 33, so that the augmented-Lagrangian multiplier estimates λ directly. The minimization follows Algorithm 1: 1200 Adam steps [55] with learning rate 3 × 10⁻³ and a multiplier update every 400 steps (μ₀ = 20, growth factor 2), then twelve outer method-of-multipliers iterations, each with an L-BFGS inner solve [56] of at most 400 iterations (μ from 200, doubled at each outer iteration up to 4 × 10⁶). The computation is carried out in double precision in NumPy with an explicit backward pass for the gradient Equation 43. The analytic gradient of the augmented Lagrangian agrees with complex-step differentiation [57] to a maximal relative deviation of at most 1.3 × 10⁻¹⁵ over the four widths for V ≡ 0, and with central differences to at most 2 × 10⁻¹⁰ for the potential case, where the complex-step route is not available through the root-finder. The projected solution ŵ = π(u_θ) is compared with the shooting reference on [0, r_cut] using the relative L², H¹ and sup norms. The strong-form residual of Proposition 5.11 is evaluated on 40000 nodes as ∥-Δŵ+Vŵ-λRŵ-ŵ3∥L2(ℝ3)/∥ŵ3∥L2(ℝ3), with the second derivative obtained by central differences of the analytic first derivative.

8.4.2 Results

Tables 3, 4 report the two problems for H ∈ {8, 16, 24, 32}; Figure 3 overlays the recovered profiles on the references, and Figures 7, 8 show the optimizer history and the width dependence. Four observations summarize them. First, the energy and the multiplier are recovered to relative accuracy between 7 × 10⁻⁷ and 1.2 × 10⁻⁵ (energy) and between 2.7 × 10⁻⁶ and 1.4 × 10⁻⁵ (Rayleigh multiplier) over the eight runs of the two width studies. The discrete energies lie above the reference level in every run, as they must, since F≥γ(m) on the sphere by Proposition 3.9 (iii). For V ≡ 0 the relative errors of γ and λ_*R* coincide to three digits: on Pm one has Aw=34Bw, hence γ=16Aw and λR=-Aw3m2, so γ=-m22λR identically, and the two relative errors are the same number.

Table 3. Constrained deep-Ritz solver vs. network width H for N = 3, q = 4, V ≡ 0, m = 4.347096 (reference γ = 9.448623, λ = −1).

Columns: γ_net = J(ŵ), the energy of the projected iterate; λ_R, its Rayleigh multiplier Equation 52; λ_AL, the multiplier returned by the augmented-Lagrangian update Equation 42; the relative errors of γ_net and λ_R against the reference; the relative L² and H¹ profile errors EL2, EH1 of Equation 75 on [0, r_cut]; the relative strong-form residual ϱ of Proposition 5.11; and the fiber parameter t(u_θ) of the raw iterate before projection. In all runs |g|/m² ≤ 2 × 10⁻¹¹ and |P(ŵ)|/Aŵ≤2×10-16.

Table 4. Constrained deep-Ritz solver vs. network width H for N = 3, q = 4, V = (1 + r²)⁻¹, m = 4.042578 (reference γ = 16.410753, λ = −1), with the fiber parameter computed by bisection at every evaluation.

The columns have the same meaning as in Table 3. In all runs |g|/m² ≤ 4 × 10⁻¹² and |P(ŵ)|/Aŵ≤6×10-16.

Second, the augmented-Lagrangian multiplier λ_AL delivered by the update (Equation 42) agrees with the Rayleigh quotient (Equation 52) to between 2 × 10⁻⁶ and 9 × 10⁻⁵; by Theorem 4.1 (ii) this gap is an optimization residual, and its size is consistent with the level at which the L-BFGS inner solves terminate. The mass constraint is satisfied to |g|/m² ≤ 2 × 10⁻¹¹.

Third, the profiles are recovered to a relative L² error between 8 × 10⁻⁴ and 1.6 × 10⁻³ (V ≡ 0) and between 7 × 10⁻⁴ and 2.2 × 10⁻³ (potential case), with H¹ and sup-norm errors of the same order. The profile errors are two to three orders of magnitude larger than the energy errors, as expected at a minimizer, where the energy error is of second order in the profile error. The strong-form residual, which involves second derivatives of the network, lies between 1.4 × 10⁻³ and 2 × 10⁻²; by Remark 5.13 it bounds the error up to a local constant, and its observed ratio to the L² profile error is between 2 and 10. The fiber parameter of the raw iterate lies between 1.27 and 1.52 in all runs (Tables 3–5), i.e., the optimizer leaves the iterate well away from the Pohozaev manifold along its orbit, as Proposition 3.9 (i) allows, and the projection π brings it back exactly: the Pohozaev functional of the projected iterate is at round-off level, |P(ŵ)|/Aŵ≤10-15, which checks the fiber computation and nothing else (Proposition 5.10).

Table 5. Spread of the results over five random initializations of the network parameters at fixed width H = 32 (seeds 0 to 4; the runs of Tables 3, 4 are seed 0).

Each entry gives the interval between the smallest and the largest value observed; the columns have the same meaning as in Table 3.

Fourth, the width dependence is mild. All quantities improve markedly from H = 8 to H = 16 and then fluctuate without a clear trend: the energy, multiplier and profile errors within a factor of about two for V ≡ 0 and, in the energy error of the potential case, by up to an order of magnitude (H = 24 being the least accurate run of that sweep, H = 32 the most accurate), and the strong-form residual by up to a factor of four. Table 5 shows that this fluctuation is of the same nature as the spread over five random initializations at fixed H = 32: the energy error varies between 3.4 × 10⁻⁶ and 5.3 × 10⁻⁶ and the profile error between 7 × 10⁻⁴ and 1.2 × 10⁻³ for V ≡ 0, and between 7 × 10⁻⁷ and 3.8 × 10⁻⁶ and between 7 × 10⁻⁴ and 1.3 × 10⁻³ for the potential case. The accuracy at these widths is therefore set by the termination of the nonconvex optimization on the flat orbit directions of Remark 3.10, not by the approximation capacity of the ansatz, and the width study should be read as evidence of robustness rather than as a convergence rate. The optimizer history in Figure 7 shows the mechanism of Section 4.2: the Adam warm-up brings the energy error to 21% with the constraint only approximately satisfied; the first L-BFGS outer iteration reduces it below 1%, and the subsequent multiplier updates take it below 10⁻⁴ after about 490 checkpoints, below 10⁻⁵ after about 950, and to 4.4 × 10⁻⁶ at the end.

Taken together, the experiments substantiate the computational claims of the study. The manufactured-solution family is exact to machine precision (Table 1); the radial references satisfy the Pohozaev and test-function identities to 10⁻⁷–10⁻⁸ (Equation 80, Section 8.2); the ground-state branch exhibits the proved monotonicity of γ(m) and the blow-up γ(m) → +∞ as m → 0⁺, with the exact power law γ ∝ m⁻² (Table 6, Figure 6); and the constrained deep-Ritz solver recovers the solution, multiplier, and energy of two L²-supercritical normalized-solution problems, one of which requires the root-finding fiber projection for the potential (Tables 3 to 5, Figure 3), with the multiplier emerging from the mass constraint as Theorem 4.1 predicts. The theory of Sections 3 to 5 is thereby confirmed in practice.

Table 6. Comparison of candidate approaches for the normalized-solution problem Equation P_*m* in the L^*p*-supercritical regime.

A check mark denotes native support, a dash denotes its absence or the need for external machinery.

9 Comparison with existing approaches

It is useful to situate the present method among the numerical and deep learning approaches that could, in principle, be applied to Equation P_*m*, and to clarify why the normalized-solution problem in the supercritical regime falls outside their scope. Table 6 summarizes the comparison across the five features that make Equation P_*m* difficult: enforcement of the prescribed mass, recovery of the unknown multiplier λ, treatment of the degenerate p-Laplacian, the whole-space domain, and the mountain-pass (saddle) nature of the supercritical energy.

Classical grid-based solvers, such as finite-element or radial shooting methods, are accurate for a fixed frequency and provide the reference used in Section 8. However, they are tied to a mesh, scale poorly with dimension, and do not by themselves impose a prescribed mass or return the branch m↦(γ, λ) without an outer continuation. The normalized gradient flow of Bao and Du [29] imposes the mass by projection at every step and recovers λ as a Rayleigh quotient, but it is a descent method for the energy on the sphere and therefore targets mass-subcritical ground states, which are minimizers; in the supercritical regime the sought state is a saddle along the dilations and the flow does not converge to it. The mountain pass algorithm of Choi and McKenna [26], the local minimax method of Li and Zhou [27], and the Nehari manifold optimization of Chen et al. [28] compute saddle points of a given energy on a mesh. The latter minimizes the Nehari manifold after a scaling retraction, which is the same device as our fiber reduction with the Pohozaev manifold in place of the Nehari manifold. None of them treats a prescribed mass with an unknown multiplier, and they are formulated on bounded domains. Among neural solvers, the physics-informed networks of Raissi et al. [31] and the deep Galerkin method of Sirignano and Spiliopoulos [32] minimize the strong residual and therefore require λ as an input, differentiate the degenerate flux twice, and admit the spurious minimizer u ≡ 0 unless the mass is separately enforced. The deep Ritz method of Yu [30] minimizes an energy and shares our first-order treatment of the flux, but in its standard form it targets coercive, unconstrained problems rather than a mass-constrained saddle. The deep eigenvalue solver of Han et al. [34] addresses an unknown spectral parameter, but for linear self-adjoint operators via a stochastic representation that does not extend to the quasilinear, mass-constrained setting. The high-dimensional method of Han et al. [33] is designed for parabolic equations via backward stochastic differential equations and is structurally different. The homotopy networks of Huang et al. [35] and the Newton-informed operator of Hao et al. [36] compute multiple solutions of nonlinear elliptic equations with prescribed coefficients on bounded domains, and are complementary in that they could enumerate further critical points once the multiplier is known. Operator-learning architectures, the DeepONet of Lu et al. [37] and the Fourier neural operator of Li et al. [38], learn solution maps across data and are complementary to our per-instance solver, as discussed in the outlook of Section 7; on their own they are surrogates that neither enforce the constraint nor certify existence. Against this landscape, the distinguishing feature of the present method is that it is built on the variational structure of the problem itself: the mass constraint is native, the multiplier is recovered as the constraint's Lagrange multiplier with an exact identification at stationary points, the degenerate flux is handled via the first-order energy with regularization, the whole space is treated with a decaying ansatz with a quantified truncation bound, and the supercritical saddle is reached by the fiber reduction onto the Pohozaev manifold. To our knowledge, no previously reported method combines these five capabilities for a quasilinear normalized-solution problem.

10 Discussion, limitations, and outlook

We have built and analyzed a mesh-free method that computes normalized solutions of the quasilinear Schrödinger problem (Equation P_*m*) in the L^*p*-supercritical regime, and we have validated it against exact and high-accuracy references. The method follows the variational proof at every stage: the functional it minimizes is the fiber-reduced energy of Definition 3.8, whose minimum on the mass sphere is the mountain-pass level by Proposition 3.9;; the constraint it enforces is the mass sphere S(m); the multiplier it returns is identified by Theorem 4.1 with the Rayleigh quotient of the projected iterate; and the output is the fiber projection onto the Pohozaev manifold of Proposition 3.4. The error analysis of Section 5 bounds the error of the computable objective, proves convergence of the energies of quasi-minimizers to γ(m) and, under the compactness property of the existence theory, strong convergence of the projected iterates to a ground state, and delimits what the residual diagnostics measure; the experiments of Section 8 confirm the theory's qualitative predictions quantitatively.

10.1 Limitations

Four limitations should be stated plainly. First, the numerical study is confined to p = 2, N = 3, q = 4, radial ansatz functions and two radial potentials, where independent radial references exist; the degenerate case p ≠ 2, non-radial potentials and non-radial ansatz functions, though covered by the analysis, are not exercised numerically here, and for V ≡ 0 a non-radial ansatz would require recentering the iterates to compensate for translation invariance. Second, the convergence theorem assumes the compactness property (H_*c*) for the sequence generated by the scheme; we prove it in the radial potential-free case and otherwise rely on the mechanism of the existence theory, which applies to the sequences constructed there. Third, the accuracy of the prototype solver, of order 10⁻³ in the profile and 10⁻⁵ to 10⁻⁶ in the energy and the multiplier, is set by the termination of the nonconvex optimization along the flat orbit directions and not by the approximation theory; reaching the refinement-limited regime of Theorem 5.5 would require a globalization strategy such as multiple restarts and a second-order inner solver with a projection onto the Pohozaev manifold at every iteration. Fourth, the operator-learning extension of Section 7 is an outlook: the selection, continuity and compactness questions listed there are open, and no operator has been trained. None of these limitations affects the central contribution, a variationally consistent, analyzed and verified numerical realization of the existence theory.

10.2 Outlook

The natural continuations are the degenerate p-Laplacian (using the regularized flux of Section 4.3 and the manufactured family of Theorem 6.3, stated for general p); a systematic study across the admissible potential class of 1 with the root-finding fiber projection of Section 4.1; non-radial computations with Monte-Carlo quadrature; and the questions raised in the outlook on the solution map. Each is a concrete computational project for which the present study supplies the mathematical foundation and the verified building block.

11 Conclusion

Normalized solutions of quasilinear Schrödinger operators in the L^*p*-supercritical regime are known to exist via variational arguments that are, by their nature, non-constructive. This study has supplied a computational counterpart. We derived the variational framework (weak formulation, Lagrange multiplier rule, Pohozaev identity via L^*p*-preserving dilations, and the fiber reduction with its existence, uniqueness and regularity theory); we constructed a constrained deep-Ritz method for the reduced energy, in which the spectral parameter emerges as the multiplier of the mass constraint and is identified exactly at stationary points of the discrete problem; we carried out an error analysis of the computable objective and proved convergence of the discrete quasi-minimizers, together with residual-based diagnostics whose meaning is stated precisely; and we verified the framework numerically, reproducing the theory's monotonicity and blow-up laws exactly and recovering two supercritical normalized solutions, their multipliers and energies, including one for a nonzero potential. Throughout, we have kept the epistemic boundary sharp: the method provides quantitative evidence for and exploration of the existence theory, and complements.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Symbol glossary

Table 7.

Footnotes

Edited by: Nenad Manojlović, University of Algarve, Portugal

Reviewed by: Feng Liu, China Aerodynamics Research and Development Center, China

Fanfei Meng, China Agricultural University, China