UNIFORM ERROR ESTIMATES OF A FINITE DIFFERENCE METHOD FOR THE KLEIN-GORDON-SCHR¨ODINGER SYSTEM IN THE NONRELATIVISTIC AND MASSLESS LIMIT REGIMES

. We establish a uniform error estimate of a ﬁnite diﬀerence method for the Klein-Gordon-Schr¨odinger (KGS) equations with two dimensionless pa- rameters 0 < γ ≤ 1 and 0 < ε ≤ 1, which are the mass ratio and inversely proportional to the speed of light, respectively. In the simultaneously nonrelativis- tic and massless limit regimes, i.e., γ ∼ ε and ε → 0 + , the KGS equations con-verge singularly to the Schr¨odinger-Yukawa (SY) equations. When 0 < ε ≪ 1, due to the perturbation of the wave operator and/or the incompatibility of the initial data, which is described by two parameters α ≥ 0 and β ≥ − 1, the solution of the KGS equations oscillates in time with O ( ε )-wavelength, which requires harsh meshing strategy for classical numerical methods. We propose a uniformly accurate method based on two key points: (i) reformulating KGS system into an asymptotic consistent formulation, and (ii) applying an integral approximation of the oscillatory term. Using the energy method and the limiting equation via the SY equations with an oscillatory potential, we establish two independent error bounds at O ( h 2 + τ 2 /ε ) and O ( h 2 + τ 2 + τε α ∗ + ε 1+ α ∗ ) with h mesh size, τ time step and α ∗ = min { 1 ,α, 1+ β } . This implies that the method converges uniformly and optimally with quadratic convergence rate in space and uniformly in time at O ( τ 4 / 3 ) and O ( τ 1+ α ∗ 2+ α ∗ ) for well-prepared ( α ∗ = 1) and ill-prepared (0 ≤ α ∗ < 1) initial data, respectively. Thus the ε -scalability of the method is τ = O (1) and h = O (1) for 0 < ε ≤ 1, which is signiﬁcantly better than classical methods. Numerical results are reported to conﬁrm our error bounds. Finally, the method is applied to study the con- vergence rates of KGS equations to its limiting models in the simultaneously nonrelativistic and massless limit regimes.

1. Introduction. We consider the coupled Klein-Gordon-Schrödinger (KGS) equations which describe a system of conserved scalar nucleons interacting with neutral scalar mesons coupled through the Yukawa interaction [12,27,40]: ∆ψ(x, t) + ηφ(x, t)ψ(x, t) = 0, x ∈ R d , d = 1, 2, 3, (1.1) Here ψ represents a complex scalar nucleon field and φ is a real scalar meson field, is the Planck constant, c is the speed of light, m 1 > 0 is the mass of a nucleon, m 2 > 0 is the mass of a meson and η > 0 is the coupling constant. The KGS system describes a classical model of the Yukawa interaction between conservative complex nucleon field and neutral meson in quantum field theory [40]. It is widely applied in many physical fields, such as many-body physics [11], nonlinear plasmas and complex geophysical flows [18], nonlinear optics and optical communications [33] and nonlinear quantum electrodynamics [31].
For scaling the KGS system (1.1), introducẽ where x s , t s = 2m1x 2 s and φ s = x −d/2 s √ 2m1 are length unit, time unit and meson field unit, respectively, to be taken for the nondimensionalization of the KGS (1.1) via (1.2). Plugging (1.2) into (1.1) and removing all '∼', we get the following dimensionless KGS system as    i∂ t ψ(x, t) + ∆ψ(x, t) + λφ(x, t)ψ(x, t) = 0, x ∈ R d , t > 0, ε 2 ∂ tt φ(x, t) − ∆φ(x, t) + γ 2 ε 2 φ(x, t) − λ|ψ(x, t)| 2 = 0, (1.3) where γ := m2 2m1 is the mass ratio, λ = If one sets the dimensionless length unit x s = 2cm1 , then ε = 1, which corresponds to the classical regime. This choice of x s is appropriate when the wave speed is at the same order of the speed of light. However, a different choice of x s is more suitable when the wave speed is much smaller than the speed of light. We remark here that the choice of x s determines the observation scale of time evolution of the system and decides which phenomena can be resolved by discretization on specified spatial/temporal grids and which phenomena is 'visible' by asymptotic analysis.
Different parameter regimes could be considered for the KGS system(1.3) which is displayed in Figure 1.1: • Standard regime, i.e., ε = 1 and γ = 1 (⇐⇒ x s = 2cm1 and m 2 = 2m 1 ): there were extensive analytical and numerical studies for the KGS equations  (1.3) with ε = γ = 1 in the last two decades. For the well-posedness, we refer to [12][13][14][15]19]; for the attractors and asymptotic behavior of the system, we refer to [7,16,17,25,26,29]; and for plane, solitary, and periodic wave solutions, we refer to [9,21,36] as well as the references therein. For the numerical part, many efficient numerical methods have been proposed for the KGS system, such as the finite difference method [30,38,41], the conservative spectral method [39], the time-splitting spectral method [5], trigonometric spectral method [20], the discrete-time orthogonal cubic spline collocation method [37], the Chebyshev pseudospectral multidomain method [10], and the symplectic and multi-symplectic methods [22][23][24]. • Massless limit regime, i.e., ε = 1 and 0 < γ ≪ 1 (⇐⇒ x s = 2cm1 and m 2 ≪ m 1 ): the KGS system (1.3) converges -regularly -to the Schrödingerwave equations i∂ t ψ(x, t) + ∆ψ(x, t) + λφ(x, t)ψ(x, t) = 0, with quadratic convergence rate in terms of γ. Any numerical methods for the KGS equations (1.3) in the standard regime can be applied in this regime. • Nonrelativistic limit regime, i.e., γ = 1 and 0 < ε ≪ 1: by taking the ansatz wherez denotes the complex conjugate of a complex-valued function z, the KGS (1.3) converges -singularly -to the Schrödinger equations [6], i.e., (ψ, z) satisfies either the Schrödinger equations with wave operator [6]    i∂ t ψ(x, t) + ∆ψ(x, t) + λ e it/ε 2 z(x, t) + e −it/ε 2 z(x, t) ψ(x, t) = 0, (1.8) or the Schrödinger equations [6]   (1.9) In addition, a multiscale time integrator Fourier pseudospectral method was proposed in [6] and it was proved that the method converges in space and time with exponential and linear convergence rates, respectively, which are uniformly for 0 < ε ≤ 1. • Simultaneously nonrelativistic and massless limit regimes, i.e., γ ∼ ε and 0 < ε ≪ 1, the KGS system (1.3) converges -singularly -to the Schrödinger-Yukawa (SY) equations, which was rigourously analyzed in [2]. To our best knowledge, there is no rigorous numerical analysis for different numerical methods for the KGS system (1.3) in this regime, especially on how the error bound depends on the small parameter ε ∈ (0, 1]. In this paper we consider the KGS equations (1.3) in the simultaneously nonrelativistic and massless limit regimes, i.e., 0 < ε ≪ 1 and γ = δε with δ > 0 a fixed constant which is independent of ε. For simplicity of notation, we choose δ = 1 and λ = 1, in which case we denote the functions as (ψ ε , φ ε ) in (1.3) and the system reads as Similar to the properties of the Zakharov system [3,28,32], the solution of the KGS equations (1.10) propagates highly oscillatory waves at wavelength O(ε) and O(1) in time and space, respectively, and rapid outgoing initial layers at speed O(1/ε) in space. This highly temporal oscillatory nature in the solution of the KGS equations (1.10) brings significant numerical difficulties, especially when 0 < ε ≪ 1 [3,28,32]. For example, classical methods may request harsh meshing strategy (or ε-scalability) in order to get 'correct' oscillatory solutions when ε ≪ 1 [8,34]. Recently, we proposed and analyzed uniform accurate finite difference methods for the Zakharov system [3] and Klein-Gordon-Zakharov system [4] in the subsonic limit regime by adopting an asymptotic consistent formulation. The main aim of this paper is to propose and analyze a finite difference method for the KGS equations, which is uniformly accurate in both space and time for 0 < ε ≪ 1.
The key ingredients rely on (i) reformulating the KGS system into an asymptotic consistent formulation and (ii) using an integral approximation of the oscillatory term. Other techniques include the energy method, cut-off technique for treating the nonlinearity and the inverse inequalities to bound the numerical solution, and the limiting equation via a Schrödinger-Yukawa system with an oscillatory potential.
The rest of the paper is organized as follows. In Section 2, we recall the singular limit of the KGS system in the nonrelativistic and massless limit regimes and introduce an asymptotic consistent formulation for the KGS equations. In Section 3, we present a finite difference method and state our main results. Section 4 is devoted to the details of the error estimates. Numerical results are reported in Section 5 to confirm our error bounds. Finally some conclusions are drawn in Section 6. Throughout the paper, we adopt the standard Sobolev spaces as well as the corresponding norms and denote A B to represent that there exists a generic constant C > 0 independent of ε, τ , h, such that |A| ≤ C B.
2. Singular limit of the KGS equations in the nonrelativistic and massless limit regimes. In this section, we recall the limit behavior of the KGS system (1.10) when ε → 0 and give an asymptotic consistent formulation for (1.10).

2.1.
Convergence of the KGS system to the Schrödinger-Yukawa equations. Formally setting ε → 0 in the KGS equations (1.10), one can get the following nonlinear Schrödinger-Yukawa (SY) system [2,29]: (2.12) It can be derived from (2.12) that φ 0 (x, t) satisfies where I is the identity operator. Setting t = 0 in (2.13), we get (2.14) Multiplying the first equation in (2.12) by ψ 0 (x, t) and subtracting from its conjugate, we obtain where f denotes the complex conjugate of f and c 0 is the so-called current [2] for the SY equations (2.12) with its definition as Differentiating (2.13) with respect to t, setting t = 0 and noticing (2.15), we get Based on the above results, the initial data (ψ 0 , φ ε 0 , φ ε 1 ) can be decomposed as where α ≥ 0 and β ≥ −1 are parameters describing the incompatibility of the initial data of the KGS equations (1.10) with respect to that of the SY equations (2.12) in the nonrelativistic and massless limit regimes such that the Hamiltonian (1.5) is bounded, ω 0 (x) and ω 1 (x) are two given real functions independent of ε. Due to the perturbation of the wave operator 'ε 2 ∂ tt ' or the inconsistency of the initial data, the solution of the KGS equations would display high oscillation in time at O(ε)wavelength with amplitude at O(ε min{2,α,1+β} ), and propagate rapid outspreading initial layers at speed O(1/ε) in space. To illustrate the temporal oscillation and rapid outgoing wave phenomena, Figure 2.2 shows the solutions φ ε (x, 1), φ ε (1, t) of the KGS system (1.10) for d = 1 and the initial data 20) and χ being the characteristic function, α = β = 0 in (2.18) for different ε, which was obtained numerically by the exponential-wave-integrator and time-splitting sine pseudospectral method on a bounded interval [−200, 200] with homogenous Dirichlet boundary condition [5].  2.
2. An asymptotic consistent formulation. Inspired by the analysis concerning on the convergence between the Zakharov system (ZS) and the limiting cubically Schrödinger equation [28] and the uniform method for solving the ZS in the subsonic limit regime [3], we introduce where ϕ ε (x, t) = (−∆ + I) −1 |ψ ε | 2 (x, t), and ω ε (x, t) represents the initial layer caused by the incompatibility of the initial data (2.18), which is the solution of the linear wave-type equation Substituting (2.21) into the KGS equations (1.10), we can reformulate it into an asymptotic consistent formulation as (2.23) The advantage of this formulation is that the main oscillatory wave with amplitude at O(1) in φ ε , which is caused by the inconsistency of the initial data, is now removed by the initial layer ω ε in (2.22), which is easy to solve separately. In the nonrelativistic and massless limit regimes, i.e.,

(2.24)
Similar to the convergence of the Zakharov system to the nonlinear Schrödinger equation in the subsonic limit [28], we can obtain the following result concerning on the convergence from the KGS system (2.23) to the SY-OP (2.24) where 0 < T < T * with T * > 0 being the maximum common existence time of the solutions of the KGS system (2.23) and the SY-OP (2.24) and C T is a positive constant independent of ε. To illustrate this, Figure 2.3 depicts the convergence behavior between the solutions of the KGS equations (2.23) and the SY-

3.
A finite difference method and main results. In this section, we present a finite difference scheme for the reformulated KGS equations (2.23) and give its uniform error bounds.

3.1.
A uniformly accurate finite difference method. For simplicity of notation, we only present the numerical method for the KGS system on one space dimension, and extensions to higher dimensions are straightforward. Practically, similar to most works for computation of the Zakharov-type equations [3,30], (2.23) is truncated on a bounded domain Ω = (a, b) with homogeneous Dirichlet boundary condition: where ω ε (x, t) is defined as the solution of (2.22) with homogeneous Dirichlet boundary condition for d = 1, is the solution of the SY-OP equations with homogeneous boundary condition Choose a mesh size h := ∆x = b−a M with M being a positive integer and a time step τ := ∆t > 0. Denote the grid points and time steps as equipped with inner products and norms defined as The finite difference operators are the standard notations as: To simplify notations, for a function E(x, t), and a grid function In this paper, we consider the finite difference discretization of (3.26) as following where we apply an average of the oscillatory potential ω ε over the interval [t k−1 , t k+1 ] Meanwhile, the initial condition is discretized as Choice of the first step value. By Taylor expansion, we get ψ ε,1 j as where by (3.26), Noticing (2.18), the above approximation for ψ ε,1 In such case, in order to make sure ψ ε,1 is uniformly bounded for ε ∈ (0, 1], τ has to be taken as τ ε −β/2 , which is too restrictive. To rescue this, we replace ψ 2 (x) above by a modified version [3] which yields the first step value with second order accuracy as In practical computation, µ ε,k j in (3.31) can be obtained by solving the linear wavetype equation (3.27) via the sine pseudospectral discretization in space followed by integrating in time in phase space exactly [3] as , and 3.2. Main results. Let T * > 0 be the maximum common existence time for the solutions of the KGS system (3.26) and the SY-OP equations (3.28). Then for any fixed 0 < T < T * , according to the known results in [29,32], it is natural to assume that the solution (ψ ε , ϕ ε , χ ε ) of the KGS (3.26) and the solution ( ψ ε , ϕ ε ) of the SY-OP (3.28) are smooth enough over Ω T := Ω × [0, T ] and satisfy 1, where α * = min{1, α, 1 + β} ∈ [0, 1]. Moreover, we assume the initial data satisfies Then one can obtain Define the error functions e ε,k ψ , e ε,k ϕ and e ε,k Thus by taking the minimum among the two error bounds for ε ∈ (0, 1], we obtain a uniform error estimate for α ≥ 0 and β ≥ −1, For s ≥ 0, y 1 , y 2 ∈ C, define ρ B (s) = s ρ s B , with B = (M 0 + 1) 2 , and g(y 1 , y 2 ) = y 1 + y 2 2 Then ρ B (s) is globally Lipschitz and Set ψ ε,k = ψ ε,k , χ ε,k = χ ε,k , k = 0, 1, and define ψ ε,k , χ ε,k ∈ X M as following Here ( ψ ε,k , ϕ ε,k , χ ε,k ) can be viewed as another approximation of (ψ ε , ϕ ε , χ ε )| t=t k . Define the error function e ε,k ψ,j , e ε,k ϕ,j , e ε,k χ,j ∈ X M for k ≥ 0 as (4.43) For the local truncation, we have the following error bounds.
Proof. By Taylor expansion, we have By (3.26) and using Taylor expansion, we get for j ∈ T M and 1 ≤ k ≤ T τ − 1, Note that by (3.31), we have Accordingly, by the assumption (A) and (3.36), we can conclude that Applying δ + x to ξ ε,k and using the same approach, we can get |δ + x ξ ε,k j | h 2 + τ 2 ε . Similarly, we obtain Applying δ c t to η ε,k j , we have Finally, it can be easily deduced that Thus the proof is completed.
For the initial step, we have the following estimates.
Remark 4.1. The error bounds in Theorem 3.1 are still valid in high dimensions, e.g., d = 2, 3, provided that an additional condition on the time step τ is added The reason is due to the discrete Sobolev inequality [8,35]: where ψ h is a mesh function over Ω with homogeneous Dirichlet boundary condition.

5.
Numerical results. In this section, we present numerical results for the KGS equations (1.10) by our proposed finite difference method. Furthermore, we apply the method to numerically study convergence rates of the KGS equations to its limiting models (2.12) and (2.24) in the nonrelativistic and massless limit regimes. In order to do so, we take d = 1 in (1.10) and the initial condition is set as (2.19).

Accuracy test.
We mainly consider two types of initial data: Case I. well-prepared initial data, i.e., α = 1 and β = 0; Case II. ill-prepared initial data, i.e., α = 0 and β = −1.  Practically, the problem is truncated on an interval Ω ε = −30 − 1 ε , 30 + 1 ε , which is large enough such that the truncation error of (3.26) to the original whole space problem (2.23) can be ignorable due to the homogeneous Dirichlet boundary condition. Due to the rapid outspreading waves with wave speed O 1 ε (cf. Figure 2.2(b)) and the homogeneous Dirichlet boundary condition truncated at the boundary, the computational domain Ω ε has to be chosen as ε-dependent. The computational ε-dependent domain can be fixed as ε-independent if one applies other where e ε,k ψ,j = ψ ε (x j , t k ) − ψ ε,k j and φ ε,k j = ϕ ε,k j + χ ε,k j + ω ε (x j , t k ) for j ∈ T M . The "exact" solution is obtained by the phase space analytical solver & time splitting spectral method [5] with very small mesh size h = 1/32 and time step τ = 10 −6 . The errors are displayed at t = 1. For spatial error analysis, we set the time step τ = 10 −5 such that the temporal error can be neglected; for temporal error analysis, the mesh size h is set as h = 2.5 × 10 −4 such that the spatial error can be ignored. Figure 5.4(a) depicts the spatial errors for Case II initial data with different mesh size h and 0 < ε ≤ 1. It clearly demonstrates that our proposed finite difference method is second order accurate in space, which is uniformly for all ε ∈ (0, 1]. The results for other initial data are analogous, e.g., different α ≥ 0 and β ≥ −1 and thus are omitted for brevity. Figure 5.4(b), Tables 5.1 and 5.2 present the temporal errors for Case I and II initial data for different time step τ and 0 < ε ≤ 1, respectively. Figure 5.4(b) shows the temporal errors of ψ ε for Case I initial data, which suggests that the method is uniformly second order accurate for the nucleon field ψ ε with well-prepared initial data. While for the messon field φ ε , the upper and lower triangle parts of Table  5.1 suggest the error at O(τ 2 /ε) and O(τ 2 + ε 2 ), respectively. There is a resonance
6. Conclusion. We presented a uniformly accurate finite difference method for the Klein-Gordon-Schrödinger (KGS) equations in the nonrelativistic and massless limit regimes -parameterized by a dimensionless parameter 0 < ε ≤ 1 -which is inversely proportional to the speed of light. When 0 < ε ≪ 1, the solution of KGS equations propagates highly oscillatory waves in time and rapid outspreading waves in space. Our method was designed by reformulating KGS system into an asymptotic consistent formulation and applying an integral approximation for the oscillating term. By using the energy method and the limiting model, we established two independent error bounds, which depend explicitly on the mesh size h, time step τ and the parameter 0 < ε ≤ 1. From the two error bounds, a uniform error estimate was obtained, which is uniformly accurate at second order in space and at least first order in time. Numerical experiments suggest that the error bounds are sharp. By adopting our numerical method, we observed that the Schrödinger-Yukawa system with an oscillatory potential approximates the KGS system quadratically in the nonrelativistic and massless limit regimes.