Electronic Journal of Differential Equations, Vol. 2025 (2025), No. 22, pp. 1–23. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu, https://ejde.math.unt.edu DOI: 10.58997/ejde.2025.22 EXISTENCE AND MULTIPLICITY OF SOLUTIONS TO QUASILINEAR DIRAC-POISSON SYSTEMS MINBO YANG, FAN ZHOU Abstract. In this article, we study the existence and multiplicity of solutions of the quasilinear Dirac-Poisson system i 3∑ k=1 αk∂ku− aβu− ωu− ϕu = h(x, |u|)u, x ∈ R3, −∆ϕ− ε4∆4ϕ = u2, x ∈ R3, where ∂k = ∂/∂xk, k = 1, 2, 3; a > 0 is a constant; α1, α2, α3 and β are 4× 4 Pauli-Dirac matrices; the operator ∆4 is the 4-Laplacian operator, defined as ∆4ϕ := div(|∇ϕ|2∇ϕ); and h(x, |u|)u describes the self-interaction. We prove the existence of the least energy solutions for the critical case and obtained that there exist finitely many critical points under certain conditions by variational methods. Additionally, we demonstrate the convergence behavior of solutions as ε tends to zero. 1. Introduction and results This study considers the Dirac system i ℏ c ∂tψ + iℏ 3∑ k=1 αk∂kψ −mcβψ − φβψ = f(x, |ψ|)ψ, x ∈ R3, −div(|∇φ| − b|∇φ|2)|∇φ| = (βψ)ψ, x ∈ R3, (1.1) where ψ denotes the wave function of the state of an electron, φ is the gauge potential of the electromagnetic field, ℏ symbolizes Planck’s constant, m > 0 means the mass of the electron, c is the speed of light, ∂k = ∂/∂xk, k = 1, 2, 3, b is a modify parameter, α1, α2, α3 and β are 4× 4 Pauli-Dirac matrices: β = ( I 0 0 −I ) , αk = ( 0 σk σk 0 ) , k = 1, 2, 3; σ1 = ( 0 1 1 0 ) , σ2 = ( 0 −i i 0 ) , σ3 = ( 1 0 0 −1 ) . 2020 Mathematics Subject Classification. 35Q40, 35J92, 49J35. Key words and phrases. Quasilinear Dirac-Poisson system; strongly indefinite problem; least energy solutions; asymptotic behavior. ©2025. This work is licensed under a CC BY 4.0 license. Submitted November 26, 2024. Published March 3, 2025. 1 2 M. YANG, F. ZHOU EJDE-2025/22 It is simple to verify that β, α1, α2 and α3 satisfy the following anticommutative relation β = β∗, αk = α∗ k, αkβ + βαk = 0, αkαl + αlαk = 2δkl, for k, l = 1, 2, 3. The Dirac-Poisson system is fundamental to relativistic quantum electrodynam- ics. It describes the complex interaction of a spin-1/2 particle with its electro- magnetic field. It plays a crucial role in quantum electrodynamics and is applied in various scientific fields, including quantum cosmology, nuclear physics, atomic physics, and gravitational physics ([29, 31]). This system plays a vital role in comprehending the quantum interactions between particles and their associated electromagnetic fields. It provides a theoretical foundation that has significantly advanced our comprehension of these phenomena. It has also been adapted to tackle specific issues in classical electrodynamics, especially those involving the infinities associated with point particles. In this regard, the Born-Infeld electromagnetic theory [10, 22] provides a non- linear alternative to Maxwell’s theory to address these infinities. A quasi-linear Dirac-Poisson system is derived by replacing the standard Maxwell’s Lagrangian density with that of Born-Infeld [23]. Born-Infeld’s theory, parameterized by b, presents a Lagrangian density in square root form. It extends classical Maxwell’s theory nonlinearly and ensures the finiteness of electric fields, thus avoiding the in- finite field issues around point particles in classical electrodynamics. Additionally, the Dirac-Born-Infeld action is used to describe D-brane dynamics in superstring theory [1], demonstrating the system’s adaptability and significance in modern the- oretical physics, bridging quantum electrodynamics with advanced string theory concepts. Many researchers have explored the solutions of Dirac equations and systems since Gross’ groundbreaking study [25] on the local existence and uniqueness of solutions for autonomous systems. Over the following decades, the nonlinear Dirac equation has garnered significant attention due to its importance in theory and application. Researchers have used variational methods to explore solutions’ ex- istence, multiplicity, and other properties based on different assumptions about potential and nonlinearity, [7, 19, 16]. When coupled with some other theories, nonlinear Dirac systems always become a nonlocal challenge. Early on, Balabane et al. [6] paved the way by transforming the Dirac equation into a planar differential system and proving the existence of a sequence of solutions. Ding and Xu [20] further analyzed stationary semi-classical solutions with general subcritical self-coupling nonlinearity. Ding and Ruf [18] also studied the multiplicity of semi-classical solutions of a nonlinear Maxwell-Dirac sys- tem, with the number of solutions described by the ratio of maximum and behaviour at infinity of the potentials. Chen et al. [11] imposed local conditions on the po- tential V , assuming it is locally Hölder continuous, ∥V ∥L∞ < a, and there exists a bounded domain Λ ∈ R3 such that ω := minΛ V < min∂Λ V.Moreover, they showed that a massive Dirac equation with critical growth has at least catMδ (M) solutions. [16] researched a nonlinear Dirac equation in space-dimension n, obtaining the ex- istence of m pairs of solutions for any ε < εm. In [21], the authors studied the multiplicity of nonlinear Dirac-Klein-Gordon systems, also describing the number EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 3 of solutions by the ratio of maximum and behaviour at infinity of the potentials. Benhassine [8] demonstrated the existence and multiplicity of stationary solutions in the asymptotically quadratic and super-quadratic cases using variational meth- ods. Recently, Alves et al. [3] complemented the results found in [17], proving that the number of global minimum points of V is directly related to the number of solutions when ε is small. The stationary wave solution of system (1.1) is a solution of the form ψ(t, x) = u(x)e −iξt ℏ , φ(x) = ϕ(x), where ξ and t are real numbers, w : R3 → C4. It is clear that (ψ,φ) solves (1.1) if and only if (u, ϕ) solves the system i 3∑ k=1 αk∂ku− aβu− ωu− ϕu = h(x, |u|)u, x ∈ R3, −∆ϕ− ε4∆4ϕ = u2, x ∈ R3, (1.2) where, for simplicity, we take b = ε4, ℏ = 1, a = mc > 0, ω = ξ c to be con- stants, functional h satisfied f(x, eiθ|u|) = eiθh(x, |u|), and the operator ∆4 is the 4-Laplacian operator, defined as ∆4ϕ := div(|∇ϕ|2∇ϕ). As a modified version of the Dirac-Maxwell system, the existence and concen- tration of minimum energy solutions for subcritical nonlinearities were recently discussed in [32]. It is a logical next step to inquire whether similar results can be achieved for quasilinear Dirac-Poisson systems. We also note several research results concerning the multiplicity and concentration phenomena of solutions for quasilin- ear problems (see [5, 4]). These works prompt us to investigate whether analogous results can be established regarding the multiplicity of solutions for quasilinear Dirac-Poisson systems. Consequently, we aim to explore the existence of minimum energy solutions and the multiplicity of solutions in quasilinear Dirac systems with critical nonlinear- ities. Here the critical exponent is 3, given by the relevant Sobolev embedding H1/2(R3,C4) ↪→ L3(R3,C4). To be specific, we consider the critical case: h(x, |u|)u = K1(x)g(|u|)u+K2(x)|u|u. (1.3) Let Ainf , Asup denote the infimum and suprermum on the whole space, respec- tively, for any function A defined in R3. Writing G(s) = ∫ s 0 g(t)tdt, we assume the nonlinear potentials satisfy the following: (A1) g(0) = 0, g ∈ C1(0,∞), g′(s) > 0 for s > 0, and there exist p ∈ (2, 3), c1 > 0 such that g(s) ≤ c1(1 + sp−2) for s ≥ 0; (A2) There exist q ≥ 2, θ > 2, and c0 > 0 such that c0s q ≤ G(s) ≤ 1 θg(s)s 2 for all s > 0; (A3) Kj ∈ C1(R3) with Kj(x) ≥ lim|x|→∞Kj(x) := kj,∞ > 0, for all x ∈ R3, j = 1, 2, and 1 < k2 := K2,sup K2,inf < R(c0, q,K1,inf ,K2,inf), 4 M. YANG, F. ZHOU EJDE-2025/22 where R(c0, q,K1,inf ,K2,inf) = (S(a2 − (ω∗)2) a2 )1/2( (c0qK1,inf) 2 q−2 6γqK2 2,inf )1/3 , ω∗ = max{ω, 0}, S is the best Sobolev constant such that S|u|26 ≤ |∇u|22, and γq is the least energy (which is attained [17]) of the equation i 3∑ k=1 αk∂ku− aβu− ωu = |u|q−2u. Our first result concerns the existence and concentration behaviour of the least energy solutions for quasilinear Dirac-Poisson systems with critical growth. We also established the existence of the limit problem and the decay properties of the solutions. Theorem 1.1. Assume ω ∈ (−a, a), h is of the form (1.3), and conditions (A1)– (A3) are satisfied. Then the following Dirac-Poisson system admits at least one least energy solution (u0, ϕ0) in ∩s≥2,r≥2W 1,r loc (R3,C4)×W 1,s loc (R3,R), iα · ∇u− aβu− ωu− ϕu = K1(x)g(|u|)u+K2(x)|u|u, x ∈ R3, −∆ϕ = u2, x ∈ R3. (1.4) If additionally ∇Kj ,j = 1, 2 are bounded, there exist C, c > 0 such that |u0(x)| ≤ C exp(−c|x|) for all x ∈ R3. Theorem 1.2. Assume ω ∈ (−a, a), h is of the form (1.3), and conditions (A1)– (A3) are satisfied. Then (1.2) admits at least one least energy solution (uε, ϕε) in ∩s≥2,r≥2W 1,r loc (R3,C4) ×W 1,s loc (R3,R) for any ε > 0. If additionally ∇Kj, j = 1, 2 are bounded, these solutions have the following properties: (1) There exist C, c > 0 such that |uε(x)| ≤ C exp(−c|x|) for any x ∈ R3; (2) The solutions (uε, ϕε) → (u0, ϕ0) in H 1 ×D1,2 as ε→ 0+. Regarding the existence of multiple solutions, we have the following results. Theorem 1.3. Assume ω ∈ (−a, a), h is of the form (1.3), and (A1)–(A3) are satisfied. For any positive integer N , there exist k∞ and m(c0, q,N,K1,inf ,K2,inf), if k∞ ≤ k2 < m(c0, q,N,K1,inf ,K2,inf), system (1.2) has at least N pairs of solutions (uε,n, ϕε,n) in ∩s≥2,r≥2W 1,r loc (R3,C4)× W 1,s loc (R3,R) for any ε > 0. If additionally ∇Kj, j = 1, 2 are bounded, there exist C, c > 0 such that |uε,n(x)| ≤ C exp(−c|x|) for any x ∈ R3. The mathematical challenges in quasi-linear Dirac-Poisson Systems are multifac- eted. Firstly, the quasi-linearity of the second equation regarding ϕ adds another layer of difficulty, as its solution lacks an explicit formula and homogeneity. Sec- ondly, the system’s strong indefiniteness means the Dirac operator’s spectrum is unbounded and contains essential spectrums, leading to a lack of a positive qua- dratic term in the energy functional of equation (1.2). Additionally, the Morse index and co-index are infinite at any critical point of this functional. Further- more, critical growth and a lack of compactness further intensify the complexity of the problem. To overcome these challenges, we will employ the reduction method EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 5 introduced by Ackermann within an appropriate variational framework. By us- ing Ding’s critical point theorems in [7], we consider the solutions of the equation as critical points of the energy functional Φε associated with system (1.2), finally proving the existence and multiplicity of solutions. The remainder part of this paper is organized as follows. In Section 2, we establish the variational framework, define the energy functionals, and recall the critical point theorems that are pivotal to our analysis. Subsequently, in Section 3, we demonstrate some preliminary results. Ultimately, in Section 4, we finish the proofs of our main results. 2. Variational framework This section aims to establish an appropriate variational setting, introduce the energy functionals, and remind the reader of the critical point theorems. We will study the ground state solutions obtained as critical points of an energy functional Φε associated with problem (1.2). Let D1,p := D1,p(R3,R) denote the Banach space defined as the completion of the test functions C∞ 0 (R3,R) with respect to the Lp-norm of the gradient provided by ∥v∥pD1,p = ∫ R3 |∇v|pdx, for p ≥ 2. Remembering the Sobolev inequality S|v|26 ≤ |∇v|22, and considering the em- beddings of D1,2(R3) and D1,4(R3) into L6(R3) and C∞ 0 (R3) respectively, we can provide equivalent characterizations as follows: D1,2(R3) = {v ∈ L6(R3) : |∇v| ∈ L2(R3)}, D1,4(R3) = {v ∈ C∞ 0 (R3) : |∇v| ∈ L4(R3)}. We define D(R3) := D1,2(R3) ∩D1,4(R3), which is a Banach space equipped with the norm ∥φ∥D := |∇φ|2 + |∇φ|4. For symbolic simplicity, let α := (α1, α2, α3) and α ·∇ := ∑3 k=1 αk∂k. Then system (1.2) can be written as iα · ∇u− aβu− ωu− ϕu = h(x, |u|)u, x ∈ R3, −∆ϕ− ε4∆4ϕ = u2, x ∈ R3. (2.1) We will write A0 := iα · ∇ − aβ, Aω := A0 − ω denote the self-adjoint operator on L2 := L2(R3,C4) with domain D(Aω) ⊂ H1 := H1(R3,C4). Let σ(Aω) and σc(Aω) signify the spectrum and continuous spectrum of Aω, respectively. Fourier analysis implies that σ(Aω) = σc(Aω) = R\(−(a+ ω), a− ω). Notice that the space L2(R3,C4) has an orthogonal decomposition: L2(R3,C4) = L+ ⊕ L−, u = u+ + u−, where A0 is positive definite on L+ and negative definite on L−. Let E := D(|Aω|1/2) = H1/2. Then E constitutes a Hilbert space equipped with a norm and inner product. For u, v ∈ E, the inner product is defined as (u, v) := ℜ(|Aω|1/2u, |Aω|1/2v)2, 6 M. YANG, F. ZHOU EJDE-2025/22 and the induced norm is ∥u∥ = (u, u)1/2, where |Aω| and |Aω|1/2 represent the absolute value of Aω and the square root of |Aω| respectively. Since σ(Aω) = R\(−(a+ ω), a− ω), we infer the inequality (a± ω)|u∓|22 ≤ ∥u∓∥2, for all u± ∈ E±. The space E can be decomposed as E = E− ⊕ E+ with E± = E ∩ L±. These subspaces are orthogonal with respect to both the (·, ·) and (·, ·)2 inner prod- ucts. This decomposition further induces a natural decomposition of Lp for every p ∈ (1,∞) [19, Proposition 2.1], thus there exists dp > 0 satisfying dp|u±|pp ≤ |u|pp for all u ∈ E ∩ Lp. (2.2) For a proof of the next lemma we refer to [7, Lemma 3.4]. Lemma 2.1. For any q ∈ [2, 3], the space E is continuously embedded in Lq(R3,C4). For any s ∈ [1, 3), E is compactly embedded in Ls loc(R3,C4). Namely, there exists constant sq > 0 such that |u|q ≤ sq∥u∥, for all u ∈ E. It is easy to see that system (2.1) is variational, and its solutions are critical points of the C2 functional Jε(u, ϕ) on E × D, defined by Jε(u, ϕ) = 1 2 (∥u+∥2 − ∥u−∥2)− 1 2 ∫ R3 ϕu2dx− F (u) + 1 4 ∫ R3 |∇ϕ|2dx+ 1 8 ε4 ∫ R3 |∇ϕ|4dx, (2.3) where u = u+ + u−, F (u) := ∫ R3 H(|u|) dx, H(|u|) = ∫ |u| 0 h(x, t)t dt. Observe that, for every u ∈ E, there exists a unique ϕuε ∈ D, which satisfies −∆ϕuε − ε4∆4ϕ u ε = u2. (2.4) For the rest of this article, ϕuε will denote the unique solution of equation (2.4), which satisfies the equation∫ R3 |∇ϕuε |2dx+ ε4 ∫ R3 |∇ϕuε |4dx = ∫ R3 ϕuεu 2dx. (2.5) For convenience, we define the operator ϕε : E → D(R3) by ϕε(u) = ϕuε for any ε > 0 fixed. By Hölder inequality and Sobolev inequality, for every u ∈ E, we have |∇ϕε(u)|22 + ε4|∇ϕε(u)|44 = ∫ R3 ϕu2dx ≤ S−1/2|∇ϕε|2|u|212 5 , which implies that |∇ϕε|2 ≤ S−1/2|u|212 5 . Thus, |∇ϕε(u)|22 + ε4|∇ϕε(u)|44 ≤ S−1|u|412 5 ≤ S−1S4 12 5 ∥u∥4. (2.6) For ease of representation, we define the functional Γε : u ∈ E 7→ 1 4 ∫ R3 |∇ϕε(u)|2dx+ 3 8 ε4 ∫ R3 |∇ϕε(u)|4dx ∈ R, EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 7 we have Γε(u) ≤ 3 8 S−1|u|412/5 ≤ 3 8 S−1S4 12/5∥u∥ 4. By inserting equation (2.5) into the functional (2.3), we can express Φε(u) := Jε(u, ϕε(u)) = 1 2 (∥u+∥2 − ∥u−∥2)− Γε(u)− F (u). In particular, we have (Φε) ′(u)v = ∂uJε(u, ϕε(u))v + ∂ϕε ϕε(u)ϕ ′ ε(u)v = ∂uJε(u, ϕε(u))v. We also deduced that (Φε) ′(u)v = (u+ − u−, v)−ℜ ∫ R3 ϕε(u)uv̄ dx−ℜ ∫ R3 h(x, |u|)uv̄ dx. We collect some valuable properties for the nonlocal term ϕε and Γε blew. Their proofs can be found in [32]. Lemma 2.2. For each ε > 0, ϕε and Γε we have the following properties: (1) ϕε maps bounded sets into bounded sets; (2) The map u 7→ Γε(u) is of class C2 in E, and its derivative satisfies Γ′ ε(u)v = ∫ R3 ϕε(u)uvdx, for all u, v ∈ E; (3) Γε is non-negative, weakly sequentially lower semi-continuous, Γ′ ε is weakly sequentially continuous; (4) If un → u in E, then, for every fixed ε > 0, Γε(un) → Γε(u) and Γ′ ε(un)un → Γ′ ε(u)u. It is easy to see that Φε ∈ C2(E,R) and critical points of Φε are weak solutions of system (2.1). To study the asymptotic behaviour of the solutions, similarly, for the limit system iα · ∇u− αβu− ωu− ϕu = h(x, |u|)u, −∆ϕ = u2, (2.7) we define functional Φ0 as Φ0(u) := 1 2 (∥u+∥2 − ∥u−∥2)− Γ0(u)− F (u), where Γ0(u) := 1 4 ∫ R3 ϕ0(u)u 2dx, and ϕ0(u) = ∫ R3 u2(y) |x−y|dy = 1 |x| ∗ u 2 is the unique solution in D1,2 such that ∫ R3 |∇ϕ0|2dx = ∫ R3 ϕ0u 2dx. Regarding the non-local term in the limit system (2.7), there are also the follow- ing properties, and their proofs can be referred to in [20]. Lemma 2.3. ϕ0 and Γ0 have the following properties: (1) ϕ0 maps bounded sets into bounded sets; (2) The map u 7→ Γ0(u) is of class C2 in E, and its derivative satisfies Γ′ 0(u)v = ∫ R3 ϕ0(u)uvdx, for all u, v ∈ E; (2.8) 8 M. YANG, F. ZHOU EJDE-2025/22 (3) Γ0 is non-negative, weakly sequentially lower semi-continuous, Γ′ 0 is weakly sequentially continuous; (4) If un → u in E, then Γ0(un) → Γ0(u) and Γ′ 0(un)un → Γ′ 0(u)u. To establish our results, we recall some abstract critical point theorems, see [7, 15]. Assume X,Y are Banach Spaces with X being separable and reflexive and set E = X ⊕ Y . Let S ⊂ X∗ be a countable dense subset. Let P be the family of semi-norms on E consisting of all semi-norms ps : E = X ⊕ Y → R, ps(x+ y) = |s(x)|+ ∥y∥, s ∈ S. Denote by TP the topology on E induced by P. Let Tw∗ be the weak∗-topology of E∗. For a functional Φ : E → R and numbers a, b ∈ R we write Φa := {u ∈ E : Φ(u) ≤ a}, Φa := {u ∈ E : Φ(u) ≥ a}, and Φa a := Φa ∪ Φa. Assume (A4) Φ ∈ C1(E,R), Φ : (E, TP) → R is upper semi-continuous, and Φ′ : (Φa, TP) → (E∗, Tw∗) is continuous for every a ∈ R; (A5) For any c > 0, there exists γ > 0, such that ∥u∥ ≤ γ∥u+∥ for all u ∈ Φc; (A6) There exists r > 0 with ρ := inf Φ(SrY ) > Φ(0) = 0 where SrY := {y ∈ Y : ∥y∥ = r}; (A7) For any e ∈ Y \{0}, there exists R with R > r > 0, such that supΦ(∂Q) ≤ ρ, where Q := {y = x+ te : x ∈ X, t > 0, ∥y∥ < R}; (A8) There exist a finite dimensional subspace Y0 ⊂ Y and R > r such that, for E0 := X × Y0 and B0 := {u ∈ E0 : ∥u∥ ≤ R}, we have supΦ(E0) <∞ and supΦ(E0\B0) < inf Φ(BrY ). We say sequence {un} ⊂ E is a (C)c sequence for Φ ∈ C1(E,R), if Φ(un) → c and (1+∥un∥)Φ′(un) ∈ 0. We say Φ satisfies the (C)c condition if any (C)c sequence for Φ has a convergent subsequence. A sequence {un} is considered a (PS)c-sequence of functional Φ if Φ(un) tends to c and Φ′(un) tends to 0. We say Φ satisfies the (PS)c condition if any (PS)c-sequence has a convergent subsequence. To prove the existence of the ground state solution, we will use the following critical point theorem. Lemma 2.4 ([15, Theorem 4.5]). Assume that conditions (A4)–(A7) are satisfied, then the functional Φ possesses a (C)c-sequence with ρ ≤ c ≤ supΦ(Q). Now we consider the set M(Φc) of maps g : Φc → E with the properties: (1) g is P-continuous and odd; (2) g(Φa) ⊂ Φa for all a ∈ [ρ, b]; (3) each u ∈ Φc has a P-open neighborhood O ⊂ E such that the set (id − g)(O ∩ Φc) is contained in a finite dimensional linear subspace. We define the pseudo-index of Φc by ψ(c) := min{genn(g(Φc) ∩ SrY ) : g ∈ M(Φc)} ∈ N0 ∪ {∞}, where gen(·) denotes the usual symmetric index. Additionally, set for d > 0 fixed M0(Φ d) := {g ∈ M(Φd) : g is a homeomorphism from Φd to g(Φd)}. We define for c ∈ [0, d], ψd(c) := min{gen(Φc ∩ SrY ) : g ∈ M0(Φ d)}. Then, by definition, we have ψ(c) ≤ ψd(c) for all c ∈ [0, d]. The following theorem plays a crucial role in proving the existence of multiple solutions. EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 9 Theorem 2.5 ([15, Theorem 4.6]). Let the assumptions (A4)–(A6), (A8) be satis- fied, and assume that Φ is even and satisfy the (C)c-condition for c ∈ [ρ, supΦ(E0)]. Then Φ has at least n := dimY0 pairs of critical points with critical values given by ci = inf{c ≥ 0 : ψ(c) ≥ i} ∈ [ρ, b], i = 1, . . . , n. If Φ has only finitely many critical points in Φ supΦ(E0) ρ , then ρ < c1 < c2 < · · · < cn ≤ supΦ(E0). We are going to use these theorems. For our purposes, we set P = X∗, thereby making TP the product topology on E = X ⊕ Y , which is defined by the weak topology on X and the strong topology on Y . 3. Preliminaries Throughout this section, we always let the hypotheses of Theorem 1.2 be satis- fied. Next, we only prove these lemmas for the problem (1.2) for critical case (1.3), the proof of the limit problem (1.4) is similar, and most of these can be checked easily in [32]. Note that (A1) and (A2) imply that for each δ > 0, there is Cδ > 0 such that g(s) ≤ δ + Cδt p−2, G(s) ≤ δs2 + Cδs p, for all s ≥ 0. Moreover, we deduce that Ĝ(s) := 1 2 g(s)s2 −G(s) ≥ 0, for all s ≥ 0. Then, we will check the assumptions in the critical theorems before. Lemma 3.1. Under assumption (A5), there exists γ > 0, satisfying ∥u∥ ≤ γ∥u+∥ for all u ∈ (Φε)c with c > 0. Proof. We argue by contradiction. Assume that there is a positive constant c and a sequence {un} ⊂ (Φε)c such that ∥un∥2 ≥ n∥u+n ∥2, for any j ∈ N. Then 0 ≥ (2− n)∥u−n ∥2 ≥ (n− 1)(∥u+n ∥2 − ∥u−n ∥2), for n ≥ 2. Hence Φε(un) = 1 2 (∥u+n ∥2 − ∥u−n ∥2)− Γε(un)− ∫ R3 K1(x)G(|un|) dx− 1 3 ∫ R3 K2(x)|un|3dx ≤ 1 2 (∥u+n ∥2 − ∥u−n ∥2) ≤ 0. But we know that Φε(un) ≥ c > 0, which is a contradiction. □ Lemma 3.2. Let Φε satisfy (A6) and (A7), that is, Φε possess the linking structure. Then (1) there exist r > 0, ρ > 0 (independent of ε), such that Φε|B+ r ≥ 0 and Φε|∂B+ r ≥ ρ where B+ r = Br ∩ E+ = {u ∈ E+ : ∥u∥ ≤ r}; (2) for each e ∈ E+\{0}, there exist R = Re > 0, C = Ce > 0 (both indepen- dent of ε), such that Φε(u) < 0 for any u ∈ Ee\BR, and supΦε(Ee) ≤ C. Proof. (1) For all u ∈ E+ and any ε > 0, we know that Γε(u) ≤ 3 8S −1S4 12 5 ∥u∥4, we have Φε(u) = 1 2 ∥u∥2 − Γε(u)− ∫ R3 K1(x)G(|u|) dx− 1 3 ∫ R3 K2(x)|u|3 dx ≥ 1 2 ∥u∥2 − 3 8 S−1S4 12 5 ∥u∥4 −K2sup(δ|u|22 + Cδ|u|pp)− K2sup 3 |u|33 10 M. YANG, F. ZHOU EJDE-2025/22 ≥ 1 4 ∥u∥2 − 3 8 S−1S4 12 5 ∥u∥4 − Cδs p pK2sup∥u∥p − s33K2sup 3 ∥u∥3. (2) For any e ∈ E+\{0}, by virtue of (2.2), for u = se+ v, we obtain Φε(u) = 1 2 (∥se∥2 − ∥v∥2)− Γε(u)− ∫ R3 K1(x)G(|u|) dx− 1 3 ∫ R3 K2(x)|u|3dx ≤ 1 2 (s2∥e∥2 − ∥v∥2)−K1,inf ∫ R3 G(|u|) dx− K2,inf 3 ∫ R3 |u|3dx ≤ s2 2 ∥e∥2 − d3s 3K2,inf 3 |e|33. The proof is complete. □ Next, we turn to study the (C)c sequence of Φε. Lemma 3.3. For all c > 0, the (C)c sequences Φε is bounded in E uniformly in ε. Proof. Given {un} satisfies Φε(un) → c and (1 + ∥un∥)Φ′ ε(un) → 0 as n → ∞. Without loss of generality, assume that ∥un∥ ≥ 1. When n is sufficiently large, we have c+ 1 ≥ Φε(un)− 1 2 (Φε) ′(un)un = 1 4 ∫ R3 |∇ϕε|2dx+ ε4 8 ∫ R3 |∇ϕε|4dx+ ∫ R3 K1(x)Ĝ(|un|) dx + 1 6 ∫ R3 K2(x)|un|3dx ≥ K1,inf( 1 2 − 1 θ )θc0|un|qq + K2,inf 6 |u|33. The sequence {un} is bounded in the spaces L2, Lq and L3 with an upper bound denoted by C1, where C1 depends only on c, K1,inf , K2,inf and q. Further deductions lead to 1 ≥ ∥un∥2 − Γ′ ε(un)(u + n − u−n )−ℜ ∫ R3 K1(x)g(un)unu + n − u−n dx −ℜ ∫ R3 K2(x)|un|2u+n − u−n dx. (3.1) Set vn = un ∥un∥ , snd recall that ϕun ε satisfies −∆ϕun ε − ε4∆4ϕ un ε = u2n. Hence, for each ψ ∈ D(R3),∫ R3 ∇ ϕun ε ∥un∥ ∇ψ dx+ ε4 ∫ R3 |∇ϕun ε |2∇ϕ un ε ∥un∥ ∇ψ dx = ∫ R3 ψunvn dx. For ∥vn∥ = 1, choose s = 6q 5q−6 such that 1 q + 1 s + 1 6 = 1. Then, since 2 < s < 3, it follows that∣∣ ∫ R3 (1 + ε4|∇ϕun ε |2)∇ ϕun ε ∥un∥ ∇ψ dx ∣∣ ≤ ∣∣ ∫ R3 unvnψ dx ∣∣ ≤ |un|q|vn|s|ψ|6 ≤ S−1/2|un|q|vn|s|∇ψ|2. EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 11 Then ∣∣∇ ϕun ε ∥un∥ ∣∣ 2 ≤ ∣∣(1 + ε4|∇ϕun ε |2)∇ ϕun ε ∥un∥ ∣∣ 2 ≤ S−1/2ss|un|q Therefore |ϕun ε |6 ≤ C1ss S ∥un∥, where S is the best Sobolev constant mentioned earlier, and ss is the constant in Lemma 2.1. Note that, by Hölder inequality, when 2 < s ≤ q < 3, we have |un|ss ≤ |un| 3(q−s) q−2 2 |un| q(s−2) q−2 q ≤ C1 3(q−s) q−2 |un| q(s−2) q−2 q , and when 2 < q < s < 3, we have |un|ss ≤ |un| 3(q−s) q−3 3 |un| q(s−3) q−3 q ≤ C1 3(q−s) q−3 |un| q(s−3) q−3 q . By Lemma 2.2, we obtain that∣∣Γ′ ε(un)(u + n − u−n ) ∣∣ = ∣∣ℜ∫ R3 ϕεunu + n − u−n dx ∣∣ ≤ |ϕε|6|un|s|u+n − u−n |q ≤ C2∥un∥1+t, (3.2) where t = q(s−2) s(q−2) when s ≤ q, and t = q(s−3) s(q−3) when s > q. It is easy to see that 0 < t < 1. Notice that by (A1), there exists r1, r2 > 0 such that g(s) ≤ a−|ω| 2K1 sup , for every s < r1, and g(s) ≤ r2s p−2, for s ≥ r1. By (A2), set δ0 := p p−2 and r3 := 2θr δ0−1 2 θ−2 , for all s ≥ r1 we have gδ0(s) ≤ rδ0−1 2 g(s)s2 ≤ r3Ĝ(s). Then, for l := pq 2q−p such that 1 δ0 + 1 q + 1 l = 1, we can estimate∣∣ℜ ∫ R3 K1(x)g(un)un(u + n − u−n ) dx ∣∣ ≤ a− |ω| 2 |un|22 +K1 sup( ∫ |u|≥r1 gδ0(|un|) dx) 1 δ0 |un|q|u+n − u−n |l ≤ a− |ω| 2a ∥un∥2 + C3∥un∥, (3.3) where C3 is independent of ε. Moreover, we have∣∣ℜ ∫ R3 K2(x)|un|unu+n − u−n dx ∣∣ ≤ K2sup|un|33 ≤ C4. (3.4) Then, the combination of estimates (3.1)-(3.4) shows that a− |ω| 2a ∥un∥2 ≤ 1 + C2∥un∥1+t + C3∥un∥+ C4. Consequently, there exists a constant Λ ≥ 1 such that ∥un∥ ≤ Λ as desired. The value of Λ is independent of ε. □ Let Kε := {u ∈ E : Φ′ ε(u) = 0} be the critical set of Φε. Due to the presence of critical terms in system (2.1), the standard bootstrap argument fails to establish the regularity of finite action weak solutions. We obtain the following regularity result using the similar argument in [26]. 12 M. YANG, F. ZHOU EJDE-2025/22 Lemma 3.4. Suppose u ∈ Kε is a critical point of Φε. Then the pair (u, ϕε) is in the space ∩s≥2,r≥2W 1,s loc (R3,C4)×W 1,r loc (R3,R). Besides this, (u, ϕε) also belongs to the space L∞(R3,C4)× L∞(R3,R). Proof. Set x ∈ R3 fixed, let ρ̄ ∈ C∞ 0 (B2(x)) be arbitrary. Choose η̄ ∈ C∞ 0 (B2(x)) such that η̄ = 1 on supp ρ̄. Define the operator D to be D = iα · ∇, we deduce D(ρ̄u) = ρ̄Du+Dρ̄ · u = η̄ · ρ̄Du+Dρ̄ · η̄. Noting that Du = aβu+ ωu+ ϕεu+K1(x)g(|u|)u+K2(x)|u|u, we have Dρ̄ · u = Aω(ρ̄u)− Tε,u(ρ̄u), (3.5) where Aω := A0 − ω. For 1 < t < 3, Tε,u :W 1,t(B2(x)) → Lt(B2(x)) is defined by w 7→ η̄ · [ϕε +K1(x)g(|u|) +K2(x)|u|]w. By applying the Gagliardo-Nirenberg inequality, it follows that |ϕε|∞ ≤ C|ϕε|1/36 |∇ϕε|2/34 . Through Sobolev embedding and the inequality (2.6), we can conclude that ϕε ∈ L∞(R3,C4). By ϕε ∈ D1,4, in value of [14], it follows that ϕε ∈ C1,α loc for 0 < α < 1. Hence, we derive that lim |x|→∞ |∇ϕε| = 0. Note that a(|∇u|) := 1 + |∇u|2 belongs to the class C1(0,∞) and satisfies the inequalities −1 < inf t>0 ta′(t) a(t) ≤ sup t>0 ta′(t) a(t) <∞, and t3 ≤ ta(t) ≤ C(t3 + 1) for t > 0. In [12, Theorem 3.1] it was shown that ∥∇ϕε∥L∞ loc(R3) ≤ C∥u2∥ 1 p−1 L3,1 loc (R3) . The embedding theorems in Lorentz space show that Lp,q loc is continuously embedded into Ln,s loc , for any 0 < n < p < ∞ and 0 < q, s < ∞. Combining this with Lq,q loc = Lq loc and u2 ∈ L∞ loc, we deduce that |Dϕε(x)| ≤ C, for all x ∈ R3, where C is independent of ε. Using the Sobolev embedding W 1,t(B2(x)) ↪→ L 3t 3−t (B2(x)) and Hölder inequal- ity, it follows that Tε,u(w) ∈ Lt(B2(x)) for w ∈ W 1,t(B2(x)) and the above map is well defined. By Minkowski and Hölder inequality, the operator norm is estimated by ∥Tε,u∥W 1,t→Lt ≤ C1(|u|L3(B) + |B|t/3) for some constant C1 (depending only on t), where B := supp η̄. Since 0 /∈ σ(Aω), Aω − Tε,u :W 1,t(B2(x)) → Lt(B2(x)) is invertible when |B| is small. EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 13 Thus, from (3.5), there exists a unique solution w ∈W 1,t(B2(x)) to the equation Aωw − Tε,u(w) = Dρ̄ · u. Let us show that there is a well-defined map Tε,u : L3(B2(x)) →W−1,3(B2(x)). In fact, by Hölder inequality, we can confirm that Tε,u(w) is belongs to L 3/2(B2(x)). Moreover, considering L3/2(B2(x)) ⊂ W−1,3(B2(x)) by the Sobolev embedding theorem, the map above remains well defined. The operator norm can be estimated as follows: ∥Tε,u∥L3→W−1,3 ≤ C2(|u|L3(B) + |B|1/3). Thus, Aω − Tε,u : L3(B2(x)) →W−1,3(B2(x)) becomes invertible if |B| is small and there exists a unique solution w̃ ∈ L3(B2(x)) solves the equation Aωw̃ − Tε,u(w̃) = Dρ̄ · u. (3.6) Consequently, we have w̃ = ρ̄ · u based on (3.5). Additionally, by W 1,s(B2(x)) ↪→ L3(B2(x)) for 3 2 ≤ s ≤ 3, we conclude that w ∈W 1,s is also a L3-solution to (3.6), given 3 2 ≤ s < 3. As a result, the uniqueness of the solution leads to w = ρ̄ · u. This implies that ρ̄ · u ∈ W 1,s(B2(x)) for any s ∈ [ 32 , 3), provided that B = supp η̄ is sufficiently small. Since ρ̄ and η̄ are arbitrary, it follows that u ∈ W 1,s(B1(x)) for any s ∈ [3/2, 3). Therefore, by Sobolev embedding, we obtain u ∈ ∩s≥2L s loc(R3) and this implies u ∈ ∩s≥2W 1,s loc (R3). Moreover, regarding equation (2.4) and u2 ∈ W 1,n loc (R3), using elliptic regularity theory [24], we deduce that ϕε ∈ W 2,n loc for any integer n ≥ 2. Using the regularity of the 4-Laplacian, we have ϕε ∈ W 1,4 loc . Consequently, ϕε ∈ ∩r≥2W 1,r loc . Finally, by elliptic estimates, we have u ∈ L∞(R3). □ Next, we state the minimax scheme and recall the mountain-pass type reduction. For each ε > 0 and e ∈ E+ \ {0}, let cε denote the minimax level of Φε deduced by the linking structure [30]: cε := inf e∈E+/{0} max u∈Ee Φε(u) = inf e∈E+/{0} max u∈Êe Φε(u), where Ee = E− ⊕ Re and Êe = E− ⊕ R+e. For u = u+ + u− ∈ E fixed, define the reduction map hε : E + → E− by Φε(u+ hε(u)) = max v∈E− Φε(u+ v). Then v ̸= hε(u) ⇔ Φε(u+ v) < Φε (u+ hε(u)) . By differentiating the functional and using the convexity of the nonlinear terms, it can be verified that hε is uniquely determined. Moreover, the following is known (see [2]): (1) hε ∈ C1(E+, E−), hε(0) = 0; (2) hε is a bounded map; (3) If un ⇀ u in E+, then hε(un)− hε(un − u) → hε(u) and hε(un)⇀ hε(u). 14 M. YANG, F. ZHOU EJDE-2025/22 We define the reduced functional Iε : E + → R by Iε(u) = Φε(u+ hε(u)), and set the Nehari-Pankov manifold Nε := { u ∈ E+\{0} : I ′ε(u)u = 0 } . Then Iε ∈ C2(E+,R) and u ∈ E+ is a critical point of Iε if and only if u + hε(u) is a critical point of Φε. We will call (hε(·), Iε(·),Nε) the Mountain-Pass reduction of system (2.1). Clearly, cε = inf u∈Nε Iε(u). For the reduction fuctional, we can verify the following result. Lemma 3.5. For any ε > 0, we have: (1) Iε possesses the mountain pass structure: Iε(0) = 0 and there exist r, ρ > 0 and e0 ∈ E+ satisfy ∥e0∥ > r such that inf Iε(S + r ) > 0 and sup Iε(e0) < 0; (2) For any u ∈ E+ \ {0}, there is a unique tε = tε(u) > 0 such that tεu ∈ Nε. Moreover, {tε(u)}ε≤1 is bounded. Proof. (1) For all u ∈ E+, by the definition of hε we have Φε(u+ hε(u)) ≥ Φε(u). Hence Iε(u) = Φε(u+ hε(u)) = 1 2 ∥u∥2 + (Φε(u+ hε(u))− Φε(u))− Γε(u)− F (u) ≥ 1 2 ∥u∥2 − Γε(u)− ∫ R3 K1(x)G(|u|) dx− 1 3 ∫ R3 K2(x)|u|3dx ≥ 1 4 ∥u∥2 − 3 8 S−1S4 12 5 ∥u∥4 − Cδs p pK2sup∥u∥p − s33K2sup 3 ∥u∥3, where δ was chosen such that δ ≤ (2s22K2sup) −1 and 2 < p < 3. For any e ∈ E+ and s > 0, we have Iε(se) = Φε(se+ hε(u)) = 1 2 (∥se∥2 − ∥hε(u)∥2)− Γε(u+ hε(u))− ∫ R3 K1(x)G(|u+ hε(u)|) dx − 1 3 ∫ R3 K2(x)|u+ hε(u)|3dx ≤ 1 2 s2∥e∥2 − K2,inf 3 ∫ R3 |u|3dx ≤ s2 2 ∥e∥2 − d3s 3K2,inf 3 |e|33. (2) We just repeat the arguments in [32] to have the results. □ When Kj(x) = kj,∞, system (2.1) becomes a autonomous problem iα · ∇u− aβu− ωu− ϕu = k1,∞g(|u|)u+ k2,∞|u|u, x ∈ R3, −∆ϕ− ε4∆4ϕ = u2, x ∈ R3. In this case we use the notation Φ∞(u) and c∞, respectively, for the functional and the least energy. For Φε(u) := 1 2 (∥u+∥2 − ∥u−∥2)− Γε(u)− k1,∞ ∫ R3 G(|un|) dx− k2,∞ 3 ∫ R3 |un|3dx, EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 15 by a similar argument as in [20], it is straightforward to verify that c∞ is attained if k22,∞ < (S(a2 − (ω∗)2) a2 )3/2 (c0qk1,∞) 2 q−2 6γq , where ω∗ = max{ω, 0}. Lemma 3.6. Φε satisfied the (C)c-condition for any 0 < c ≤ c∞. Moreover, c is attained if c < (a2 − (ω∗)2 a2 )3/2 S3/2 6k2K2 2,inf , where ω∗ = max{ω, 0}. Proof. Suppose that {un} is a (C)c-sequence of Φε. Lemma 3.3 shown that {un} is bounded. It is easy to check that {un} is relatively compact for all 0 < c ≤ c∞. Without loss of generality, let us assume that there is a uε such that un converges weakly to uε in E. Now we are going to show that uε ̸= 0 for all small ε > 0. Assume by contradiction that un is vanishing, then un → 0 in Lq for q ∈ (2, 3). Notice that c+ o(1) ≥ Φε(un)− 1 3 (Φε) ′(un)un = 1 6 (∥u+n ∥2 − ∥u−n ∥2) + 1 12 ∫ R3 |∇ϕε|2dx− ε4 24 ∫ R3 |∇ϕε|4dx − ∫ R3 K1(x)( 1 3 g(|un|)u2n −G(|un|)) dx ≥ 1 6 ||Aω|1/2un|22 + o(1), we have ||Aω|1/2un|22 ≤ 6c+ o(1). Similarly, |un|33 ≤ 6c K2,inf + o(1). Moreover, o(1) ≥(Φε) ′(un)(u + n − u−n ) = ∫ R3 ⟨Aωun, u + n − u−n ⟩ dx− Γ′ ε(un)(u + n − u−n ) − ∫ R3 K1(x)g(|un|)un(u+n − u−n ) dx− ∫ R3 K2(x)|un|un(u+n − u−n ) dx. Then ∫ R3 ⟨Aωun, u + n − u−n ⟩ dx ≤ K2sup|un|23|u+n − u−n |3 + o(1). By Calderón-Lions interpolation theorem [28], we have S1/2|un|23 ≤ ||iα · ∇|1/2un|22. We denote the Fourier transform by F : L2 → L2, recall from [19] that∫ R3 ⟨iα · ∇un, un⟩ dx = ∫ R3 |ξ||Fu(ξ)|2 dξ,∫ R3 ⟨Aωun, un⟩ dx = ∫ R3 ( (a2 + |ξ|2)1/2 − ω ) |Fu(ξ)|2 dξ. 16 M. YANG, F. ZHOU EJDE-2025/22 Taking into account that inf |ξ|>0 (a2 + |ξ|2)1/2 − ω |ξ| = {( a2−ω2 a2 )1/2 , ifω > 0, 1, ifω ≤ 0, we have ∫ R3 ⟨Aωun, un⟩ dx ≥ (a2 − (ω∗)2 a2 )1/2 ∫ R3 ⟨iα · ∇un, un⟩ dx, where ω∗ = max{ω, 0}. Finally we obtain c ≥ (a2 − (ω∗)2 a2 )3/2S3/2K2,inf 6K3 2sup = (a2 − (ω∗)2 a2 )3/2 S3/2 6k2K2 2,inf , which contradicts the hypothesis. Therefore, uε ̸= 0 and (uε, ϕε) is a solution of system (2.1). □ Similarly, for system (2.7), we have a Mountain-Pass reduction (h0(·), I0(·),N0) and the similar results as before. Denote L0 be the set of all least energy solutions. Now, we need to analyze the energy levels using the reduced functional to obtain the asymptotic behaviour of the least energy solutions. Lemma 3.7. limε→0 cε = c0. Proof. First we prove that lim infε→0+ cε ≥ c0. Arguing indirectly, assume that lim infε→0+ cε < c0. By definition and Lemma 3.5 we can choose ej ∈ Nεj and δ > 0 such that max u∈Eej Φεj (u) ≤ c0 − δ as j → ∞. By [9, Lemma 3.2], for all u ∈ E, we have Φεj (u)− Φ0(u) = Γ0(u)− Γεj (u) → 0. Note that c0 ≤ I0 (ej) ≤ max u∈Eej Φ0(u). Therefore for all j sufficiently large such that |Γ0(u)− Γεj (u)| < δ 2 , we have c0−δ ≥ max u∈Eej Φεj (u) ≥ max u∈Eej Φ0(u)+Γ0(u)−Γεj (u) ≥ c0+Γ0(u)−Γεj (u) > c0− δ 2 , which is a contradiction. Next, we turn to show that lim supε→0+ cε ≤ c0. Let u = u+ + u− ∈ L0, and set e = u+. It is evident that e ∈ N0, h0(e) = u−, and I0(e) = c0. There exist a unique tε > 0 such that tεe ∈ Nε, and we have cε ≤ Iε(tεe). (3.7) By Lemma 3.5, tε is bounded. Hence, without loss of generality we can assume tε → t0 as ε→ 0+. Setting uε = tεe+h0 (tεe), wε = tεe+hε (tεe) and vε = uε−wε, we deduce that 1 2 ∥vε∥2 + (I) = Φε (wε)− Φε (uε) = Φ0 (wε)− Φ0 (uε)− Γε (wε) + Γε (uε) + Γ0 (wε)− Γ0 (uε) , EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 17 where (I) := ∫ 1 0 (1− s) (Γ′′ ε (wε + svε)[vε, vε] + Γ′′ 0(wε + svε)[vε, vε]) ds. Considering Γε (uε)− Γε (wε) = Γ′ ε (wε) vε + ∫ 1 0 (1− s)Γ′′ ε (wε + svε) [vε, vε] ds, Γ0 (wε)− Γ0 (uε) = −Γ′ 0 (uε) vε + ∫ 1 0 (1− s)Γ′′ 0 (uε − svε) [vε, vε] ds, we have 1 2 ∥vε∥2 + (I) + (II) ≤ Γ′ ε (wε) vε + ∫ 1 0 (1− s)F ′′ (wε + svε) [vε, vε] ds− Γ′ 0 (uε) vε, where (II) := ∫ 1 0 (1− s)Γ′′ 0 (uε − svε) [vε, vε] ds. So we deduce that 1 2 ∥vε∥2 + ∫ 1 0 (1− s)Γ′′ 0 (wε + svε) [vε, vε] ds ≤ |Γ′ ε (wε) vε|+ |Γ′ 0 (us) vε| . (3.8) Since tε → t0, it is clear that {uε}, {wε} and {vε} are bounded in E. Moreover, we have Γε (zε) = o(1), ∥Γε(zε)∥ = o(1) as ε→ 0+ for zε = uε, wε, vε. By (2.8), noting that |Γ′ 0 (uε) vε| → 0. Thus from (3.8), it follows that ∥vε∥2 → 0, that is, hε (tεe) → h0 (t0e). This, jointly with the definitions,implies Φε(wε) = Φ0(wε) + o(1) = Φ0(uε) + o(1), that is Iε(tεe) = I0(t0e) + o(1) as ε→ 0+. Then, since I0(t0e) ≤ max v∈Ee Φ0(v) = I0(e) = c0, we obtain by (3.7) that lim sup ε→0+ cε ≤ lim ε→0+ Iε(tεe) ≤ c0. Above all, we have c0 ≤ lim inf ε→0+ cε ≤ lim sup ε→0+ cε ≤ c0. □ 18 M. YANG, F. ZHOU EJDE-2025/22 4. Proof of main results Given that the proof for the limit system (2.7) parallels that of the more complex system (2.1), we choose to address the existence of the least energy solutions for system (2.1) first. This approach establishes a foundation for the subsequent proof of the limit system (2.7), as outlined in Theorem 1.1. We will thus begin by proving Theorem 1.2, which pertains to system (2.1). Proof of Theorem 1.2. Part 1. Existence of least energy solutions for system (2.1). For any ε > 0, by Lemma 2.2, assumption (A4) is satisfied. And by Lemma 3.1 and Lemma 3.2, Φε satisfies all the assumptions of Lemma 2.4. Hence Φε has a (C)cε- sequence {un} with ρ ≤ cε ≤ supΦε(Ee ∩ BR). By Lemma 3.3, {un} is bounded in E. Therefore, up to a subsequence, there is a point uε such that un ⇀ uε in E. Since we have assumed k2 := K2sup K2,inf < (S(a2 − (ω∗)2) a2 )1/2( (c0qK1,inf) 2 q−2 6γqK2 2,inf )1/3 , it follows that γq < (a2 − (ω∗)2 a2 )3/2S3/2(c0qK1,inf) 2 q−2 6k32K 2 2,inf . Consider the least energy γc0qK1,inf ,q of the following equation iα · ∇u− aβu− ωu = c0qK1,inf |u|q−2u. It is easy to see that γc0qK1,inf ,q = (c0qK1,inf) −2 q−2 γq, and γc0qK1,inf ,q satisfies γc0qK1,inf ,q < (a2 − (ω∗)2 a2 )3/2 S3/2 6k32K 2 2,inf . Observe that cε ≤ γc0qK1,inf ,q, Lemma 3.6 shows that cε is attained by some point, denoted as uε. By Lemma 3.4 we see that solution (uε, ϕε) is in ∩s≥2,r≥2W 1,s loc ×W 1,r loc . Part 2. Decay estimate of solutions for system (2.1). By Lemma 3.4 it can be observed that u ∈ L∞(R3,C4). Expressing (2.1) as Du = aβu+ ωu+ ϕεu+K1(x)g(|u|)u+K2(x)|u|u. Operating the operator D on both sides and using the property that D2 = −∆, we derive a relation ∆u = a2u− (ω + ϕε +K1(x)g(|u|) +K2(x)|u|)2u −D(ϕε +K1(x)g(|u|) +K2(x)|u|)u. Let sgnu = { ū/|u|, if u ̸= 0, 0, if u = 0, and referring to Kato’s inequality [13], it can be found that ∆|u| ⩾ ℜ(△u · sgnu) = ℜ((a2u− (ω + ϕε +K1(x)g(|u|) +K2(x)|u|)2u EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 19 −D(ϕε +K1(x)g(|u|) +K2(x)|u|)u) · sgnu). Further, observing that ℜ[D(K1(x)g(|u|) +K2(x)|u|)u(sgnu)] = 0, we obtain ∆|u| ≥ a2|u| − (ω + ϕε(x) +K1(x)g(|u|) +K2(x)|u|)2|u| − |Dϕε| · |u|. (4.1) Recall that in Lemma 3.4, we have |ϕε|∞ ≤ C and |∇ϕε| ≤ C, where C is independent of ε. Subsequently, it follows from (4.1) that there exists a constant M > 0 (independent of ε) such that ∆|u| ≥ −M |u|. Then applying the maximum principle (see [27]), we can conclude that |uε(x)| ≤ C exp(−c|x|) for all x ∈ R3, and C, c is independent of ε. Part 3. Asymptotic behaviour of solutions for system (2.1). Suppose (uε, ϕε) is a pair of least energy solution for system (2.1), then uε satisfies Φε(uε) = cε, (1 + ∥uε∥)(Φε) ′(uε) → 0 and ρ ≤ cε ≤ supΦε(Ee ∩BR). With the independence of ε in Lemma 3.3, {∥uε∥} is bounded uniformly in ε, and there exists a point u0 ∈ E such that uε ⇀ u0 in E as ε→ 0+. As the proof above, we see that uε → u0 in E, then we also have uε → u0 in Lq for all q ∈ [2, 3]. Recall that |Aω(uε − u0)|2 ≤ |ϕεuε − ϕ0u0|2 + |K1(x)(g(|uε|)uε − g(|u0|)u0)|2 + |K2(x)(|uε|uε − |u0|u0)|2. Since |ϕεuε − ϕ0u0|2 ≤ |uε|3|ϕε − ϕ0|6 + |uε − u0|3|ϕε|6, and |K1(x)(g(|uε|)uε − g(|u0|)u0)|2 ≤ K1 sup|(g(uε)− g(u0))uε|2 + δK1 sup|uε − uo|2 + CδK1 sup|u0|2(p−2) ∞ |uε − u0|2, we deduce |Aω(uε − u0)|2 → 0 as ε → 0+, that is, uε → u in H1. From [9, Lemma 3.2] we have that ϕε(uε) → ϕ0(u0) in D 1,2(R3), εϕε(uε) → 0 in D1,4(R3). Supposing v ∈ C∞ 0 (R3), supp(v) ⊂ K, and K is compact. Recall that (u+ε − u−ε , v) = ℜ ∫ R3 ϕε(uε)uεv̄ dx+ ℜ ∫ R3 K1(x)g(|uε|)uεv̄ dx+ ℜ ∫ R3 K2(x)|uε|uεv̄ dx. We will pass to the limit as ε → 0+ in the above identity. Let us examine each term individually. Of course (u+ε − u−ε , v) → (u+0 − u−0 , v). Since ϕε(uε) → ϕ0(u0) in L6(R3), uε → u0 in L12/5(K) and v ∈ L12/5(K). It is easy to see that ℜ ∫ R3 ϕε(uε)uεv̄ dx→ ℜ ∫ R3 ϕ0(u0)u0v̄ dx. 20 M. YANG, F. ZHOU EJDE-2025/22 We also have ℜ ∫ R3 K1(x)g(|uε|)uεv̄ dx→ ℜ ∫ R3 K1(x)g(|u0|)u0v̄ dx, ℜ ∫ R3 K2(x)|uε|uεv̄ dx→ ℜ ∫ R3 K2(x)|u0|u0v̄ dx. As a result, we deduce that (u+0 − u−0 , v)−ℜ ∫ R3 ϕ0(u0)u0v̄ dx−ℜ ∫ R3 K1(x)g(|u0|)u0v̄ dx −ℜ ∫ R3 K2(x)|u0|u0v̄ dx = 0. This means that (u0, ϕ0) solves system (2.7) with energy Φ0(u0) = Φ0(u0)− 1 2 (Φ0) ′(u0)u0 = 1 4 ∫ R3 |∇ϕ0|2dx+ ∫ R3 K1(x)Ĝ(|u0|) dx+ 1 6 ∫ R3 K2(x)|u0|3dx. By applying Fatou’s lemma, we deduce that c0 ≤ 1 4 ∫ R3 |∇ϕ0|2dx+ ∫ R3 K1(x)Ĝ(|u0|) dx+ 1 6 ∫ R3 K2(x)|u0|3dx ≤ lim inf ε→0+ (1 4 ∫ R3 |∇ϕε|2dx+ ε4 8 ∫ R3 |∇ϕε|4dx+ ∫ R3 K1(x)Ĝ ( |uε| ) dx + 1 6 ∫ R3 K2(x) |uε|3 dx ) = lim inf ε→0+ Φε(uε) ≤ lim sup ε→0+ Φε(uε) ≤ c0. Consequently, Lemma 3.7 shows that lim ε→0+ Φε (uε) = lim ε→0+ cε = Φ0(u0) = c0. Now we can conclude that (uε, ϕε) converges in H 1 ×D1,2(R3) to (u0, ϕ0) which is a least energy solution of the system (2.7). □ Proof of Theorem 1.1. Part 1. Existence of least energy solutions for system (2.7). From the explicit expression of Γ0, in comparison with the non-local term Γε of the original problem (2.1), the non-local term Γ0 in the limit problem (2.7) possesses the same or better properties, as demonstrated in Lemma 2.4. Therefore, the series of lemmas proven for the functional Φε also hold the same conclusions for Φ0. In line with the proof of Theorem 1.2, we similarly have the functional Φ0 satisfies all the assumptions of Lemma 2.4. By assumptions (K1), (K2) and Lemma 3.6, it follows that functional Φ0 has at least one nontrivial solution that achieves the least energy. Part 2. Decay of solutions for system (2.7). Expressing (2.7) as Du = aβu+ ωu+ ϕ0u+K1(x)g(|u|)u+K2(x)|u|u. Operating the operator D on both sides, we obtain the relation ∆u = a2u− (ω + ϕ0 +K1(x)g(|u|) +K2(x)|u|)2u EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 21 −Dϕ0u−D(K1(x)g(|u|) +K2(x)|u|)u. Let sgnu = { ū/|u|, if u ̸= 0, 0, if u = 0. By Kato’s inequality [13], it can be found that ∆|u| ≥ ℜ(△u · sgnu) = ℜ((a2u− (ω + ϕ0 +K1(x)g(|u|) +K2(x)|u|)2u −D(ϕ0 +K1(x)g(|u|) +K2(x)|u|)u) · sgnu). Further, observing ℜ[D(K1(x)g(|u|) +K2(x)|u|)u(sgnu)] = 0, we obtain ∆|u| ≥ a2|u| − (ω + ϕ0(x) +K1(x)g(|u|) +K2(x)|u|)2|u| − |Dϕ0| · |u|. (4.2) Then, following the analogous arguments in [20], we can readily deduce from (4.2) that there exists a constant M > 0 such that ∆|u| ≥ −M |u|. Let Γ be a fundamental solution to −∆ + τ . We may choose Γ(x) such that |uε(x)| ≤ τΓ(x) on BR(0). Denote z = |uε| − τΓ, then ∆z = ∆|uε| − τ∆Γ ≥ τ(|uε| − τΓ) = τz for |x| ≥ R. The application of the maximum principle [27] brings us to the conclusion that |uε(x)| ≤ C exp(−c|x|), where C, c is independent of ε. □ Finally, we proceed with the proof of the existence of multiple solutions. Proof of Theorem 1.3. Obviously, Φε is even in u, and in virtue of Lemma 2.2 , 3.1 and 3.2, the assumptions (A4), (A5) and (A6) are satisfied. It remains to verify (A8) and the (C)c-conditon. For any N ∈ N+, let {en} ⊂ E+ be a standard orthogonal basis, EN := span{e1, e2, . . . , eN}. Since EN is a finite dimensional subspace, there exits cN > 0 such that for every u ∈ EN , |u|q ≥ cN∥u∥. Then for any u = u+ + u− ∈ EN ⊕E−, Φε(u) = 1 2 (∥u+∥2 − ∥u−∥2)− Γε(u)− ∫ R3 F (x, |u|) dx ≤ ∥u+∥2 − 1 2 ∥u∥2 − ∫ R3 K1(x)G(|u|) dx ≤ 1 c2N |u+|2q − c0K1,inf |u+|qq − 1 2 ∥u∥2 ≤ q − 2 q c − 2q q−2 N ( 2 qc0K1,inf ) 2 q−2 − 1 2 ∥u∥2. Clearly, there exists RN = ( 2q − 4 q )1/2c − q q−2 N ( 2 qc0K1,inf ) 1 q−2 , 22 M. YANG, F. ZHOU EJDE-2025/22 such that sup u∈EN⊕E−,∥u∥>RN Φε(u) < 0 and sup u∈EN⊕E− Φε(u) ≤ 2R2 N . We denote m(c0, q,N,K1,inf ,K2,inf) := q 6(q − 2)K2,inf (K1,infc q Nc0q 2 ) 2 q−2 (S(a2 − (ω∗)2) a2 )3/2 and T (s) := ( S(a2−(ω∗)2) a2 )3/2 1 6K2 2,infs for s > 1. Note that, T (1) > c∞ and T (s) → 0 as s → ∞. There exists k∞ > 1 such that T (k∞) = c∞. Hence for k∞ ≤ k2 < m(c0, q,N,K1,inf ,K2,inf), by Lemma 3.6 and Lemma 2.5, Φε has at least N distinct critical values. Finally, repeat the proof of Theorem 1.1, we obtain the decay estimate. □ Acknowledgments. Minbo Yang was partially supported by the NSFC (11971436, 12011530199) and by the ZJNSF (LZ22A010001, LD19A010001). References [1] H. Abe, Y. Sakamura, Y. Yamada; Matter coupled Dirac-Born-Infeld action in four- dimensional N = 1 conformal supergravity. Phys. Rev. D, 92(2015), no. 2, 025017, 8 pp. [2] N. Ackermann; A nonlinear superposition principle and multibump solutions of periodic Schrödinger equations. J. Funct. Anal., 234(2006), 277-320. [3] C. O. Alves, R. N. de Lima, A. B. Nóbrega; Existence and multiplicity of solutions for a class of Dirac equations. J. Differential Equations, 370(2023), 66-100. [4] C. O. Alves, M. Yang; Investigating the multiplicity and concentration behaviour of solutions for a quasi-linear Choquard equation via the penalization method. Proc. Roy. Soc. Edinburgh Sect. A, 146(2016), no. 1, 23-58. [5] C. O. Alves, M. Yang; Multiplicity and concentration of solutions for a quasilinear Choquard equation. J. Math. Phys., 55(2014), no. 6, 061502, 21 pp. [6] M. Balabane, T. Cazenave, A. Douady, F. Merle; Existence of excited states for a nonlinear Dirac field. Comm. Math. Phys., 119(1988), no. 1, 153-176. [7] T. Bartsch, Y. Ding; Solutions of nonlinear Dirac equations. J. Differential Equations, 226(2006), no. 1, 210-249. [8] A. Benhassine; Standing wave solutions of Maxwell-Dirac systems. Calc. Var. Partial Differ- ential Equations, 60(2021), no. 3, Paper No. 107, 20 pp. [9] K. Benmlih, O. Kavian; Existence and asymptotic behaviour of standing waves for quasilinear Schrödinger-Poisson systems in R3. Ann. Inst. H. Poincaré C Anal. Non Linéaire, 25(2008), no. 3, 449-470. [10] M. Born, L. Infeld; Foundations of the new field theory. Proc. R. Soc. Lond. A, 144(1934), 425-451. [11] Y. Chen, Y. Ding, T. Xu; Potential well and multiplicity of solutions for nonlinear Dirac equations. Commun. Pure Appl. Anal., 19(2020), no. 1, 587-607. [12] A. Cianchi, V. Maz’ya; Global gradient estimates in elliptic problems under minimal data and domain regularity. Commun. Pure Appl. Anal., 14(2015), no. 1, 285-311. [13] R. Dautray, J. Lions; Mathematical analysis and numerical methods for science and technol- ogy. vol. 3. Springer-Verlag, Berlin, 1990. [14] E. DiBenedetto; C1+α local regularity of weak solutions of degenerate elliptic equations. Nonlinear Anal., 7(1983), no. 8, 827-850. [15] Y. Ding; Variational methods for strongly indefinite problems. Interdisciplinary Mathematical Sciences, 7. World Scientific Publishing Co. Pte. Ltd., Hackensack, NJ, 2007. [16] Y. Ding, X. Dong, Q. Guo; On multiplicity of semi-classical solutions to nonlinear Dirac equations of space-dimension n. Discrete Contin. Dyn. Syst., 41(2021), no. 9, 4105-4123. [17] Y. Ding, X. Liu; Semi-classical limits of ground states of a nonlinear Dirac equation. J. Differential Equations, 252(2012), no. 9, 4962-4987. EJDE-2025/22 QUASILINEAR DIRAC-POISSON SYSTEMS 23 [18] Y. Ding, B. Ruf; On multiplicity of semi-classical solutions to a nonlinear Maxwell-Dirac system. J. Differential Equations, 260(2016), no. 7, 5565-5588. [19] Y. Ding, T. Xu, Localized concentration of semi-classical states for nonlinear Dirac equations. Arch. Ration. Mech. Anal., 216(2015), no. 2, 415-447. [20] Y. Ding, T. Xu; On semi-classical limits of ground states of a nonlinear Maxwell-Dirac system. Calc. Var. Partial Differential Equations, 51(2014), no. 1-2, 17-44. [21] Y. Ding, Y. Yu, X. Dong; Multiplicity and concentration of semi-classical solutions to non- linear Dirac-Klein-Gordon systems. Adv. Nonlinear Stud., 22(2022), no. 1, 248-272. [22] P. A. M. Dirac; An extensible model of the electron. Proc. R. Soc. Lond. A, 268(1962), 57-67. [23] D. Fortunato, L. Orsina, L. Pisani; Born-Infeld type equations for the electrostatic fields. J. Math. Phys., 43(2002), 5698-5706. [24] D. Gilbarg, N. S. Trudinger; Elliptic partial differential equations of second order. Springer, Berlin, 1998. [25] L. Gross; The Cauchy problem for the coupled Maxwell-Dirac equations. Commun. Pure Appl. Math., 19(1966), 1-5. [26] T. Isobe; Nonlinear Dirac equations with critical nonlinearities on compact Spin manifolds. J. Funct. Anal., 260(2011), no. 1, 253-307. [27] P. Rabier, C. Stuart, Exponential decay of the solutions of quasilinear second-order equations and Pohozaev identities. J. Differential Equations, 165(2000), 199-234. [28] M. Reed, B. Simon; Methods of modern mathematical physics. II. Fourier analysis, self- adjointness. Academic Press, 1975. [29] F. Schwabl; Advanced quantum mechanics. Springer, 1999. [30] A. Szulkin, T. Weth; Ground state solutions for some indefinite variational problems. J. Funct. Anal., 257(2009), no. 12, 3802-3822. [31] B. Thaller; The Dirac equation. Texts and monographs in physics. Springer, Berlin, 1992. [32] F. Zhou, Z. Shen, M. Yang; Existence and asymptotic behaviour of the least energy solutions for a quasilinear Dirac-Poisson system. Discrete Contin. Dyn. Syst. Ser. S, 16(2023), no. 11, 3427-3458. Minbo Yang School of Mathematical Sciences, Zhejiang Normal University, Jinhua, Zhejiang 321004, China Email address: mbyang@zjnu.edu.cn Fan Zhou (corresponding author) School of Mathematical Sciences, Zhejiang Normal University, Jinhua, Zhejiang 321004, China Email address: 826670487@zjnu.edu.cn 1. Introduction and results 2. Variational framework 3. Preliminaries 4. Proof of main results Acknowledgments References