{'fraction_non_alphanumeric': 0.056113332766102274, 'fraction_numerical': 0.022972857993703737, 'mean_word_length': 4.20988475177305, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 10, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Arbitrary high-order structure-preserving methods for the quantum Zakharov system

In this paper, we present a new methodology to develop arbitrary high-order structure-preserving methods for solving the quantum Zakharov system. The key ingredients of our method are: (i) the original Hamiltonian energy is reformulated into a quadratic form by introducing a new quadratic auxiliary variable; (ii) the original system is then rewritten into a new equivalent system which inherits the quadratic energy based on the energy variational principle; (iii) the resulting system is discretized by symplectic Runge-Kutta method in time combining with the Fourier pseudo-spectral method in space. The proposed method achieves arbitrary high-order accurate in time and can preserve the discrete mass and Hamiltonian energy exactly. Moreover, an efficient iterative solver is presented to solve the resulting discrete nonlinear equations. Finally, extensive numerical examples are provided to demonstrate the theoretical claims and illustrate the efficiency of our method.


Introduction
The quantum Zakharov system (QZS) [18,26] is widely used to describe the nonlinear interaction between the Langmuir waves and the ion-acoustic waves and can usually be formulated, as follows: i∂ t E(x, t) + ∆E(x, t) − ε 2 ∆ 2 E(x, t) = N (x, t)E(x, t), ∂ tt N (x, t) − ∆N (x, t) + ε 2 ∆ 2 N (x, t) = ∆|E(x, t)| 2 , x ∈ Ω, t > 0, E(x, 0) = E 0 (x), N (x, 0) = N 0 (x), N t (x, 0) = N 1 (x), x ∈ Ω. (1.1) The complex-valued function E := E(x, t) denotes the slowly varying envelope of the rapidly oscillatory electric field, and the real-valued function N := N (x, t) represents the low-frequency variation of the density of the ions. The quantum effect ε > 0 is the ratio of the ion plasma and the temperature of electrons, E 0 (x), N 0 (x) and N 1 (x) are given initial functions, and the periodic boundary condition is imposed on ∂Ω. In the special case ε = 0, it reduces to the classical Zakharov system (ZS), which have been widely applied to various physical problems, such as plasma [54,24], hydrodynamics [13], the theory of molecular chains [12] and so on. When either the electrons temperature is low or the ion-plasma frequency is high, the quantum effect can be characterized by the fourth-order perturbation with a quantum parameter ε. For more physical background, we refer the reader to [25,37,40]. The QZS (1.1) also conserves the mass and energy where △v = N t . Extensive mathematical and numerical studies have been carried out for the above QZS (1.1) in the literature. Along the mathematical front, Fang et al [16] showed that the QZS is locally well-posed in L 2 (R d ) data for dimension up to eight and globally well-posed for 1 ≤ d ≤ 5, which is different from the classical ZS where the global and local well-posedness of the Cauchy problem is known only for 1 ≤ d ≤ 3 [15,19]. The blow-up in finite time of the solution for high-dimensional classical ZS was investigated in [38]. Along the numerical front, the numerical studies of the classical or generalized ZS are very rich, such as time splitting method [3,2,31,32], scalar auxiliary variable approach [48], finite difference method [1,8,20,52], discontinuous Galerkin method [51], etc. Recently, there has been growing interest in developing accurate and efficient numerical methods for the QZS (1.1). Baumstark and Schratz [5] developed a new class of asymptotic preserving trigonometric integrators for the QZS (1.1), and the scheme converges to the classical ZS in the limit ε → 0 uniformly. Zhang [56] proposed an explicit time-splitting combined with exponential wave integrator Fourier pseudo-spectral method for the QZS (1.1). However, these mentioned schemes cannot conserve the mass (1.2) and the energy (1.3) of the QZS (1.1). It has been shown that in the numerical simulation of the collision of solitons, the solution of mass-and energy-conserving schemes cannot produce nonlinear blow-up [44,55]. Zhang and Su [57] proposed a linearly-implicit conservative compact finite difference method for the QZS (1.1), however it is shown rigorously in mathematics that the scheme is only second-order accurate in time. In [11,23,29], it is clear to see that high-order accurate structure-preserving scheme will provide much smaller numerical error and more robust than the second-order accurate one as the large time step is chose. Thus, it is desirable to propose high-order accurate mass-and energy-preserving methods for solving the QZS (1.1).
As a matter of the fact, in the past few decades, how to devise high-order accurate energy-preserving schemes for conservative systems have attracted plenty of attention. The excellent ones include Hamiltonian Boundary Value Methods (HBVMs) [6,7], energy-preserving variant of collocation methods [10,27], highorder averaged vector field (AVF) methods [34,42,39], functionally fitted energy-preserving methods [41,36,50] and energy-preserving continuous stage Runge-Kutta (RK) methods [41,49], etc. These methods can be easily extended to propose high-order accurate energy-preserving scheme for the QZS (1.1), which however cannot preserve the mass exactly [4,35]. Based on the basic principle of the structure-preserving method where the numerical method should preserve the intrinsic properties of the original problems as much as possible, it is valuable to expect that the high-order mass and energy-preserving discretizations for the QZS (1.1) will produce richer information on the continuous system. Very recently, inspired by the ideas of the invariant energy quadratization (IEQ) approach [53], a new class of high-order accurate energypreserving methods are proposed in [21,33]. Especially, the term "quadratic auxiliary variable (QAV) approach" was coined by Gong et al in [21]. The QAV approach is a special subclass of the IEQ approach and the key differences between them are: (i) the auxiliary variable introduced by the QAV approach shall be quadratic; (ii) when a high-order accurate method in time which can preserve the quadratic invariant is applied to the equivalent system, the resulting method can preserve both the quadratic invariant and the original Hamiltonian energy [21,22] instead of a modified energy [17,30]. Thus, the QAV approach will be an efficient strategy to develop high-order accurate mass-and energy-preserving scheme for the QZS (1.1), however to our knowledge, there has been no reference considering this issue.
In this paper, motivated by the QAV approach, we first introduce a quadratic auxiliary variable to transform the Hamiltonian energy into a quadratic from, and the original system is then reformulated into a new equivalent system. Finally, a fully-discrete method is obtained by using the RK method in time and the Fourier pseudo-spectral method in space for the reformulated system. We show that when the symplectic RK method is selected, the proposed method can conserve the discrete mass and original Hamiltonian energy exactly. In addition, to solve the discrete nonlinear equations of our method efficiently, a fast iterative solver is presented.
The rest of the paper is organized as follows. In section 2, we first reformulate the original QZS (1.1) into an equivalent form, and both the energy and the mass conservation laws of the reformulated system are investigated. In section 3, a new class of high-order structure-preserving schemes are proposed based on the symplectic RK method in time and the Fourier pseudo-spectral method in space, respectively. In Section 4, an efficient implementention for the proposed scheme is presented. In Section 5, extensive numerical experiments for the QZS are carried out to illustrate the capability and accuracy of the method, and show some complex dynamical behaviors. Finally, some conclusions are drawn in Section 6.

Model reformulation
Denote u = (E, N, v) T , the QSZ (1.1) can be written into the following energy-conserving system where u is the complex conjugate of u, is the skew-adjoint operator, and the Hamiltonian energy functional Next, based on the QAV approach, we first introduce a quadratic auxiliary variable, as follows: where (f, g) = Ω f gdx, for any f, g ∈ L 2 (Ω) and the L 2 -norm is defined as f = (f, f ). According to the energy variational principle, we rewritten the QZS (2.1) into a new equivalent system with the consistent initial conditions where Re(•) represents the real part of •.
At present, we focus on investigating the structure-preserving properties of the reformulated system (2.5)-(2.6).
Theorem 2.1. Under the periodic boundary condition, the system (2.5) possesses the following conservation laws where M(t) and E(t) are defined by (1.2) and (2.4), respectively.
Proof. It is clear to see that which implies (2.7). Combining the initial conditions (2.6) with fourth equation of the system (2.5), we can deduce which directly yields (2.8).
According to (2.4) and (2.5), a routine computation gives rise to the modified energy conservation law, as follows: This immediately gives (2.9) and the proof is completed.

High-order accurate mass-and energy-preserving method
In this section, a family of symplectic RK methods is first used to discretize the temporal derivatives in the system (2.5) to establish semi-discrete high-order structure-preserving schemes, and rigorously proved to maintain the original energy and mass conservation laws of the QZS (1.1). Furthermore, the Fourier pseudo-spectral method is then used to discretize the spatial derivatives for extending to the fully-discrete case.

Time semi-discretization
Let the time step τ = T J with T is the final time, and denote t n = nτ for 0 ≤ n ≤ J. Let w n and w ni be the numerical approximations of the function w(x, t) at t n and t n + c i τ , respectively. Applying an s-stage RK method to discrete the system (2.5) in time and the following semi-discrete scheme is presented. Scheme 3.1. Consider an s-stages RK method (A, b, c) can be described by the abscissa c = (c 1 , · · · , c s ) T ∈ R s , the weight vector b = (b 1 , · · · , b s ) T ∈ R s , and the coefficients matrix A = (a ij ) s i,j=1 ∈ R s×s , where c i = s j=1 a ij . For given (E n , N n , v n , q n ), the following intermediate values are calculated by (3.1) Theorem 3.1. If the coefficients of the RK method (3.1) satisfy with the periodic boundary condition, the proposed method preserves the following semi-discrete conservation laws where M n = |E n | 2 , 1 , Proof. According to Eq. (3.2), we have a ij k 1 j into the right hand side of (3.10) and using (3.6), we have Furthermore, we can deduce Then, combining the above equation with (3.11), we obtain the discrete mass conservation law (3.7).

Based on (3.2) and (3.6) together with
Noticing that It follows from (3.12) and (3.13) that (3.14) With the help of the initial condition q 0 = |E 0 | 2 , we can obtain (3.8). Similar to Eq. (3.11), we can obtain from Scheme 3.1 that Thus, together with (3.1), we can derive which immediately gives the modified energy conservation law (3.9). This proof is finished.
Under the consistent initial condition (2.6) and the periodic condition as well as the symplectic condition (3.6), the symplectic RK method conserves the original energy conservation law at the discrete-intime level, that is, Proof. By using the property (3.8) and noting that the consistent initial condition, we have q n = |E n | 2 , ∀ n ≥ 0, (3.16) which leads to This together with the modified energy conservation law (3.9), one can obtain the original energy conservation law (3.15). This completes the proof.
Thanks to q n = |E n | 2 , the original energy and mass conservation law of the symplectic RK methods can be preserved. For convenience, this new class of structure-preserving schemes is called SPRK methods. It is easy to find that the SPRK method is equivalent to the following scheme.
Remark 3.1. We should note that the Gaussian collocation method where the RK coefficients c 1 , c 2 , · · · , c s are chosen as the Gaussian quadrature nodes, i.e., the zeros of the s-th shifted Legendre polynomial d s dx s (x s (x− 1) s ) satisfies the condition (3.6) [28]. And especially, the RK coefficients for 2-stage and 3-stage Gauss collocation methods can be given explicitly [28,43] (see Tab. 1).  Remark 3.2. It is well-known that the scalar auxiliary variable (SAV) approach [46,47] is also an efficient method for developing high-order accurate structure-preserving methods for the conservative systems [11,29]. However, we should note that it is challenging for introducing a special scalar auxiliary variable to construct high-order accurate methods in time which can preserve the original energy of the system.

Full discretization
For Scheme 3.2, the Fourier pseudo-spectral method is used to discretize the spatial derivatives for extending to the fully discrete SPRK schemes. Let the domain Ω = (−L, L) × (−L, L) be uniformly partition with sizes h = h x = h y =: 2L N , where N is a positive even integer. Then, the spatial grid points ∈ Ω h be the space of grid functions defined on Ω h . Note that an element of V h can be regarded as a vector, and its basic rules of operation are the same as the vector, unless otherwise stated. For any two grid functions u, v ∈ V h , define the discrete inner product and norms as Let X j (x) and Y k (y) be the interpolation basis functions given by Then, we define S N as the interpolation space and the interpolation operator I N : C(Ω) → S N is given by Taking the second-order derivative with respect to θ (θ = x or y), and it follows from the resulting expression at the collocation points (x j , y k ) that Lemma 3.1. [45] For the matrix D θ 2 (θ = x or y), there exists the following relation 22) where F N denotes the discrete Fourier transform (DFT), and satisfies F T N = F −1 N , Here we only introduce the operator ∆ h = 1 2 I y ⊗ D 2 x + D 2 y ⊗ I x to denote the second-order Fourier differential matrix and omit the details due to save space. Interested readers are referred to [9] for details. Applying the Fourier pseudo-spectral method to Scheme 3.2, we obtain the following fully discrete scheme. (3.24) For the fully discrete scheme, we can obtain the following theorems can be carried out similarly as Theorem 3.1 and Theorem 3.2, we skip the proof and only give the results here. where Remark 3.3. On the one hand, with noting (3.26), we can deduce that the proposed method can preserve the discrete Hamiltonian energy exactly at every time step, as follows: On the other hand, we need to obtain v n at every time step for computing the discrete Hamiltonian energy (3.28). Thus, we require the zeroth Fourier coefficient on v is zero, so that the Poisson equation N t = ∆v with periodic boundary conditions only has a solution.

An efficient implementention for the proposed scheme
In this section, we will present a fast solver to solve the resulting nonlinear equations of the proposed schemes efficiently, which is mainly based on the matrix diagonalization method (see Ref. [45] and references therein) and the discrete Fourier transform (FFT) algorithm. For simplicity, we only consider the 2-stage Gauss method (i.e., s = 2) and the RK coefficients a ij , b i (i, j = 1, 2, · · · , s) are given in Table 1. For given E n , N n and ∆ h V n , 2-stage Gauss method is equivalent to

Numerical results
The purpose of this section is to test validity and efficiency of the newly proposed schemes for solving the QZS (1.1) with the periodic boundary condition, where the SPRK scheme with s-stage Gaussian collocation coefficients will be referred to as SPRK-s. Numerical examples including accuracy tests, convergence of the QZS (1.1) to the classical ZS in the semi-classical limit, dynamic simulation of QZS (1.1) in two-dimension, soliton collisions and pattern dynamics of one-dimensional QZS (1.1).  Scheme

Accuracy test
Let E n h,τ and N n h,τ be the numerical solutions obtained by the SPRK-s with mesh size h and time step τ for approximating the exact solutions E(·, t n ) and N (·, t n ). To quantify the numerical methods, we define the error functions as e ε (t n ) = E(·, t n ) − E n h,τ h,∞ , n ε (t n ) = N (·, t n ) − N n h,τ h,∞ .
Meanwhile, we also define the relative residuals on the mass and energy as where B is a constant, x 0 and v represent the initial displacement and the propagation velocity of the soliton, respectively. In this test, we choose B = 1, v = 1 2 , x 0 = 0, and the initial conditions are taken as      Table 2 reports the spatial errors of the SPRK-s at T = 1 with a very small time step τ = 1/2 12 such that the discretization error in time is negligible. We observe that the spatial errors converge exponentially. Table 3 lists the temporal errors of the SPRK-s with a fixed mesh size h = 1/2 4 , for the classical ZS at T = 1 under different τ . The results clearly show that fourth (s = 2) and sixth (s = 3) order accuracy are achieved. In addition, Figures 1 and 2 show the relative residuals on the mass and energy calculated by the SPRK-s with h = 1/2 4 , τ = 1/20 until T = 200. It is found that the SPRK-s preserve the mass and energy conservations of the classical ZS very well.   To demonstrate the advantages of our proposed high-order SPRK-s with the TS-EWI-FP method [56], we summarized the evolution of numerical errors on long-time simulations using the large time step τ = 1/20, h = 1/2 4 and the final time T = 200 in Figure 3. Compared to the TS-EWI-FP method, the errors of the SPRK-s are much smaller at the same time steps, and the SPRK-s have a good numerical performance on long time simulations.  To verify the accuracy and efficiency, we also investigated the maximum norm errors and computational time using the SPRK-s and the TS-EWI-FP method [56] at T = 1 under ε = 1/8 and h = 1/2 4 . The results are summarized in Table 6. These computational results illustrate that the errors of the SPRK-s are much smaller at the same time steps, and TS-EWI-FP method converges at the second order in time. Furthermore, generally speaking, the SPRK-s is superior to the TS-EWI-FP scheme in accuracy for fixed h and ε.
Subsequently, we apply the SPRK-s to study the convergence rate of the QZS to its limiting model, i.e. the classical ZS (ε = 0). In order to do so, we take the same initial conditions (5.3), and choose various ε = 2 −(2j+1) , j = 2, 3, 4, 5, 6, 7. Figure 6 shows the maximum norm errors between the QZS and the corresponding classical ZS at T = 10, which are calculated by using the same SPRK-s with h x = h y = 1/2 4 and τ = 1/100. It is clear to see that the solution of the QZS with (5.3) converges to the classical ZS quadratically with respect to ε. This is consistent with the exact result presented in [14]. with the initial conditions

Dynamic simulation of two-dimensional QZS
When the quantum effect is absent (ε = 0) and 2ν − ν 2 µ > 4π, mathematically, the solution of the classical ZS will blow-up at some time [38,32], and this problem have been utilized in [3,32] to demonstrate their numerical behavior. Moreover, taking quantum effects into account, we apply the SPRK-2 to investigate the dynamical behaviors of the QZS (5.4)-(5.5), and test the following cases: (i). ε = 0, µ = ν = 20; (ii). ε = 1/2 5 , µ = ν = 20; (iii). ε = 1, µ = ν = 20. be observed that a singularity indeed starts to form as time evolves, and this numerical results are in good agreement with those presented in [3,32,56]. At the same time, one can see that no blow-up occurs with the strong quantum regime ε = 1 (cf. case (iii) in Figure 9), but when we decrease ε to 1/2 5 while keeping other parameters to be the same, we observe clearly the blow-up process in Figure 8.

Dynamic simulation of one dimensional QZS
Example 5.4. The initial data are selected as (cf. [57,56]) where x = x j (j = 1, 2) are initial locations of the two solitary waves, v j (j = 1, 2) are the propagation velocity and the sign means moving to the right or left. All computations are performed by using the SPRK-2 with The simulation of the surface plots of the interaction of two solitary waves for QZS (ε = 0) and the classical ZS (ε = 0) under cases (I)-(III) are demonstrated in Figures 10-12, respectively. We can observe that all the collisions between two solitons are not elastic. When two initially well-separated solitons with opposite propagation velocities (cf. case (I) in Figure 10) or different propagation speeds (cf. case (II) in Figure 11), they collide and fuse into a new soliton with the strengthened amplitude and the narrower width. Meanwhile, the amplitude-weakened solitons with propagation speeds changed and some small radiation are generated during the collision. Subsequently, when the initial locations are not initially well-separated (cf. case (III) in Figure 12), the dynamics is considerably more complicated, it seems that a periodic perturbation on the position of some localized pulse. From Figures 10-12, the soliton-soliton collisions of the QZS are more unstable than the corresponding classical ZS after collision, and the quantum effect makes the chaos much more obvious. The larger the quantum effect is, the more obvious the spatiotemporal chaos is. Particularly, for the strong quantum regime ε = 1, there are small outgoing waves emitting before colliding, and the chaos is much more obvious than the classical one. where E 0 = (k/ √ 2)(1 + ε 2 k 2 ) is the amplitude of the pump Langmuir wave, 0 < k < √ 2E 0 and β represents a relatively small constant to emphasize that the perturbation is small.
The QZS with the initial conditions (5.7) was advanced in time by using the SPRK-2 with time step τ = 1/20 and mesh size h = 1/10 on the interval Ω = [−100, 100]. Take parameters T = 200, k = 0.7, β = 0.001. Figure 13 shows that many solitary patterns can be generated and excited through the modulational instability of unstable harmonic modes. It can be clearly seen that the QZS is more unstable than the corresponding classical ZS. In particular, for the strong quantum regime ε = 1, numerical simulation also indicates that the motion of the trains leads to more collision among the neighboring coherent solitary patterns, and fuse into the newer pattern with strengthened amplitude. This space-time evolution reveals that the spatiotemporal chaos is more obvious as the classical one. The numerical results are in good agreement with the results given in [40]. We finally remark that similar results of the SPRK-3 are obtained for Examples 5.3, 5.4 and 5.5, so for the sake of brevity, we do not list them here.

Conclusion
In this paper, we develop a novel class of high-order accurate structure-preserving methods for solving the QZS (1.1), which is based on the idea of QAV approach, the symplectic RK method as well as the standard Fourier pseudo-spectral method. We show that the proposed scheme can preserve the discrete mass and Hamiltonian energy exactly. In addition, a fast iterative solver is presented to solve the resulting nonlinear equations of the proposed scheme efficiently. Numerical experiments for the QZS (1.1) are carried out to illustrate the capability and accuracy of the new schemes. We also use our new scheme to numerically simulate the soliton-soliton interaction and the pattern dynamics of the QZS (1.1) in 1D and the blow-up of the solution of the QZS (1.1) in 2D. We numerically observe that the numerical solution of the QZS (1.1) converges to the classical Zakharov system quadratically in the semi-classical limit.