Loading [MathJax]/jax/element/mml/optable/BasicLatin.js
Research article Special Issues

Wronskian-type determinant solutions of the nonlocal derivative nonlinear Schrödinger equation

  • A nonlocal derivative nonlinear Schrödinger (DNLS) equation is analytically studied in this paper. By constructing Darboux transformations (DTs) of arbitrary order, new determinant solutions of the nonlocal DNLS equation in the form of Wronskian-type are derived from both zero and nonzero seed solutions. Periodic solitons are obtained with different parameter choices. When one eigenvalue tends to another one, generalized DTs are constructed, leading to rogue waves. Due to complex parametric constraints, the derived solutions may have singularities. Despite this, the work presented in this paper can still provide a valuable reference for the study of nonlocal integrable systems.

    Citation: Dexin Meng. Wronskian-type determinant solutions of the nonlocal derivative nonlinear Schrödinger equation[J]. AIMS Mathematics, 2025, 10(2): 2652-2667. doi: 10.3934/math.2025124

    Related Papers:

    [1] Hongying Jiao, Shuhai Zhu, Jinguo Zhang . Existence of infinitely many solutions for critical sub-elliptic systems via genus theory. Communications in Analysis and Mechanics, 2024, 16(2): 237-261. doi: 10.3934/cam.2024011
    [2] Jizheng Huang, Shuangshuang Ying . Hardy-Sobolev spaces of higher order associated to Hermite operator. Communications in Analysis and Mechanics, 2024, 16(4): 858-871. doi: 10.3934/cam.2024037
    [3] Xiulan Wu, Yaxin Zhao, Xiaoxin Yang . On a singular parabolic p-Laplacian equation with logarithmic nonlinearity. Communications in Analysis and Mechanics, 2024, 16(3): 528-553. doi: 10.3934/cam.2024025
    [4] Erlend Grong, Irina Markina . Harmonic maps into sub-Riemannian Lie groups. Communications in Analysis and Mechanics, 2023, 15(3): 515-532. doi: 10.3934/cam.2023025
    [5] Leandro Tavares . Solutions for a class of problems driven by an anisotropic (p,q)-Laplacian type operator. Communications in Analysis and Mechanics, 2023, 15(3): 533-550. doi: 10.3934/cam.2023026
    [6] Velimir Jurdjevic . Time optimal problems on Lie groups and applications to quantum control. Communications in Analysis and Mechanics, 2024, 16(2): 345-387. doi: 10.3934/cam.2024017
    [7] Zhiyong Wang, Kai Zhao, Pengtao Li, Yu Liu . Boundedness of square functions related with fractional Schrödinger semigroups on stratified Lie groups. Communications in Analysis and Mechanics, 2023, 15(3): 410-435. doi: 10.3934/cam.2023020
    [8] Ming Liu, Binhua Feng . Grand weighted variable Herz-Morrey spaces estimate for some operators. Communications in Analysis and Mechanics, 2025, 17(1): 290-316. doi: 10.3934/cam.2025012
    [9] Xiaotian Hao, Lingzhong Zeng . Eigenvalues of the bi-Xin-Laplacian on complete Riemannian manifolds. Communications in Analysis and Mechanics, 2023, 15(2): 162-176. doi: 10.3934/cam.2023009
    [10] Eleonora Amoroso, Angela Sciammetta, Patrick Winkert . Anisotropic (p,q)-Laplacian problems with superlinear nonlinearities. Communications in Analysis and Mechanics, 2024, 16(1): 1-23. doi: 10.3934/cam.2024001
  • A nonlocal derivative nonlinear Schrödinger (DNLS) equation is analytically studied in this paper. By constructing Darboux transformations (DTs) of arbitrary order, new determinant solutions of the nonlocal DNLS equation in the form of Wronskian-type are derived from both zero and nonzero seed solutions. Periodic solitons are obtained with different parameter choices. When one eigenvalue tends to another one, generalized DTs are constructed, leading to rogue waves. Due to complex parametric constraints, the derived solutions may have singularities. Despite this, the work presented in this paper can still provide a valuable reference for the study of nonlocal integrable systems.



    In this paper, we focus on the following hydrodynamically consistent Cahn-Hilliard phase field model for incompressible two-phase flows with variable density

    ρt+(ρu)=0, (1.1a)
    ϕt+uϕγΔw=0, (1.1b)
    w=Δϕ+1ε2(ϕ3ϕ), (1.1c)
    ρ(ut+uu)ηΔu+pλwϕ=0, (1.1d)
    u=0, (1.1e)

    where ϕ, w, u and p are the phase field function, the chemical potential, the velocity of flow and pressure respectively. Moreover, ρ=12(ρ1+ρ2)+ϕ2(ρ1ρ2) is the density, where ρ1 and ρ2 are the densities of the two fluids. The parameters γ, ε, η and λ are the mobility parameter related to the relaxation time scale, the interface thickness [1], the viscosity of the field, and the mixing energy density respectively. The system (1.1) is supplemented with appropriate boundary and initial conditions which are given as

    ϕ(x,0)=ϕ0(x), ρ(x,0)=ρ0(x), u(x,0)=u0(x),ϕ|Ω=0,         u|Ω=0.

    This model is also called the Cahn-Hilliard-Navier-Stokes model when the density is constant, which has many practical applications in physical and engineering, such as wetting, coating, and painting. By considering the influence of variable density, this model has broad applicability, which include highly stratified flows, interfaces between fluids of different densities and some problems of inertial confinement.

    The phase field method which has a wide range of applications is one of the main methods to deal with the fluid interface in two-phase flow modeling; see [2,3,4,5] and the references therein. It was initially developed to simulate solid-liquid phase transitions, where the interface is treated as a thin, smooth transition layer[6,7,8] to remove the singularities between two phases. The basic framework is to use phase field variables to represent the volume fraction of fluid components and then adopt a variational form to derive the model. Recently, it has increasingly attracted researchers' interest, mainly because the phase field method is superior to other available methods in some aspects of two-phase flow [9]. Various material properties or complex interface behaviors can be simulated directly by introducing suitable energy functions. Numerical solutions play a crucial role in their study and applications because the analytic solutions are usually not available.

    A few works have been devoted to the design, analysis and implementation of numerical schemes for the phase field model, although this model is very classical and canonical. We briefly review the available methods used in the phase field model. It is worth noting that there are projection/gauge/penalty methods [10,11,12], scalar auxiliary variables [13,14], linear stability [15,16], convex splitting [17,18,19], invariable energy quadratization(IEQ) [20,21], nonlinear quadratic [22], exponential time differencing [23], etc. In practical applications, we often couple the flow-field equation with the phase field equation. Typically two-phase incompressible flow models are coupled to phase-field models. The Cahn-Hilliard model is taken into account in this paper, because it is effective in the following two aspects: (ⅰ) Cahn-Hilliard models can accurately conserve the volume and dynamics; (ⅱ) the equation is one of the most important models in mathematical physics. Because of these reasons, here we use the Cahn-Hilliard equation developed in [24] to couple the incompressible flows with variable density.

    There also have been a lot of works on numerical approximates for the Cahn-Hilliard phase field model for incompressible two-phase flows with variable density. Hohenberg and Halperin proposed the model in [25] to simulate two incompressible viscous fluids with constant density. In [26], Gurtin et al. obtained the equal model by using the framework of rational continuum mechanics. A fully adaptive energy stabilization scheme is proposed in [27]. An efficient Picard iteration procedure was designed in [28] to further decouple the model. In the last years, many authors have been concerned with designing incompressible two-phase flow models with variable densities. Several efficient and energy-stable time discretization schemes for the coupled nonlinear Cahn-Hilliard phase field system with variable density are constructed; see [29]. Yang and Dong presented an energy-stable scheme in [30] for the numerical approximation of the two-phase governing equations with variable density and viscosity for the two fluids by introducing a scalar-valued variable related to the total of the kinetic energy and the potential free energy. A second-order accurate, coupled, energy-stable schemeis proposed in [31], where the Crank-Nicolson method and the IEQ method were used. In [32], the conservation scheme of the first-order energy law was established, in which the Cahn-Hilliard solver was used to decouple from the two-phase incompressible flows solver through the use of the fractional step method. Ye et al.[33] have designed a fully-decoupled type scheme to solve the Cahn-Hilliard phase field model for a two-phase incompressible fluid flow system with constant density. They only give a detailed practical implementation method and also prove the solvability. None of the various existing schemes have been subjected to error analysis, where the main difficulty lies in the delicate treatment of a several of nonlinear terms. Rigorous error estimates of models with variable densities, using an optimal order error bound, may seem to be a difficult prospect, but is a very interesting direction for future research.

    In this paper, we finally arrive at an unconditionally stable in energy, first-order time-accurate scheme for the incompressible Cahn-Hilliard two phase flows with variable density by coming up with a fractional step method. This method has an advantage over the projection method in that the original boundary conditions of the problem can be implemented in all substeps of the scheme. The popular approach to discretizing the Cahn-Hilliard phase field model (1.1b) and (1.1c) in time is based on the convex-splitting of the free energy functional, i.e., an idea that can be traced back to [32]. In the convex-splitting framework, one treats the contribution from the convex part implicitly and the contribution from the concave part explicitly. This treatment promotes the energy stability of the scheme and this property is unconditional in terms of time steps. We also give a rigorous proof of the convergence results and error estimates in the theoretical analysis. The main contribution of this paper is a rigorous error analysis, particularly under the condition that energy stability is available. To the best of the authors' knowledge, the proof developed in this article is the first to have the description of error estimation. The accuracy and stability are also demonstrated through the simulation of various numerical examples, where the challenge is in creating an efficient and easy to implement numerical scheme that preserves the energy dissipation law.

    The rest of this article is organized as follows. In Section 2, we construct an efficient time discrete scheme for variable density and derive unconditional energy stability. In Section 3, the error analysis of the semi-discrete scheme in time is provided. Some numerical experimentations are given in Section 4. Finally, conclusions are drawn in Section 5.

    For the sake of simplicity, some notations are needed for the following content. We assume that the domain ΩR2 is open, sufficiently smooth and bounded. For any two functions ϕ(x) and ψ(x), their L2 inner product on Ω is denoted by (ϕ,ψ)=Ωϕ(x)ψ(x)dx, and the L2 norm of ϕ(x) is denoted by ϕ=(ϕ,ϕ)12. Let τ>0 be the time step size and set tn=nτ for 0nN with T=Nτ. Moreover, we introduce the following spaces,

    H={uL2(Ω)|divu in Ω, un=0 on Ω},V=H10, V0={uV|divu=0 in Ω},M=L20(Ω)={qL2(Ω)|Ωqdx=0},

    where n is the outward normal vector of Ω. Next, we reformulate the system (1.1) as follows:

    ρn+1ρnτ+ρn+1un=0, (2.1a)
    ϕn+1ϕnτ+˜un+1ϕnγΔwn+1=0, (2.1b)
    wn+1=1ε2((ϕn+1)3ϕn)Δϕn+1, (2.1c)
    ρn˜un+1unτηΔ˜un+1+ρn+1(un)˜un+1λwn+1ϕn+14ρn+1(un)˜un+1=0, (2.1d)
    ρn+1un+1˜un+1τη(Δun+1Δ˜un+1)+pn+1=0, (2.1e)
    un+1=0 (2.1f)

    Remark 1. In (2.1d), the term 14ρn+1(un)˜un+1 is added to obtain the unconditional stability, as it is 0 if un=0.

    Theorem 1. (Stabilityofρ) For any τ>0 and any sequence {un}n=0,,N satisfying the boundary condition unn=0 on Ω, the solution {ρn}n=1,,N to (2.1a) satisfies

    ρN2+N1n=1ρn+1ρN2=ρ02. (2.2)

    Proof. Testing (2.1a) with 2τρn+1 gives

    ρn+12ρn2+ρn+1ρn2+2τΩ(unρn+1+12ρn+1un)ρn+1dx=0. (2.3)

    Owing to the boundary conditions on un, we note that

    Ω(unρn+1+12ρn+1un)ρn+1dx=12Ω(|ρn+1|2un)dx=12Ω|ρn+1|2unndx=0. (2.4)

    Thus, we get

    ρn+12+ρn+1ρn2=ρn2.

    Summing all indices n ranging 0 to N1, the proof is completed.

    Next, the stability of (2.1b)–(2.1e) will be proved in the theorem below. Moreover, since the kinetic energy of the fluid is 12ρnun2, it is more suitable to establish a bound based on ρnun2 than on the velocity itself. For simplicity, let us say that σn=ρn for all 1nN and σ0=ρ0.

    Theorem 2. (Stabilityofenergy) For any τ>0, (2.1b)–(2.1e) satisfy the the conditions of following energy estimates:

    σNuN2+λϕN2λε2ϕN2+λ2ε2ϕN4L4+2τγλN1n=0wn+12+λN1n=0((ϕn+1ϕn)2+12ε2ϕn+1ϕn2)+λ2ε2N1n=0(ϕn+12ϕn)2+N1n=0(σn(˜un+1un)2+σn+1(un+1˜un+1)2)+ητN1n=0(un+12+˜un+12+(un+1˜un+1)2)+λε2N1n=0(ϕn+1(ϕn+1ϕn)2)=σ0u02+λϕ02λε2ϕ02+λ2ε2ϕ04L4.

    Proof. Multiplying (2.1b) by 2τλwn+1 and integrating over Ω, we get

    2λ(ϕn+1ϕn,wn+1)+2τλγwn+12+2τλΩ˜un+1ϕnwn+1dx=0. (2.5)

    Multiplying (2.1c) by 2λ(ϕn+1ϕn) yields

    2λ(ϕn+1ϕn,wn+1)=λ2ε2(ϕn+14L4ϕn4L4)+λ2ε2(ϕn+12ϕn2)2+λε2ϕn+1(ϕn+1ϕn)2+λε2(ϕn2ϕn+12+ϕn+1ϕn2)+λ(ϕn+12ϕn2+(ϕn+1ϕn)2), (2.6)

    where we use the following identity:

    2a(ab)=a2b2+(ab)2,a3(ab)=14(a4b4+(a2b2)2+2a2(ab)2). (2.7)

    Testing (2.1a) with τ|˜un+1|2 leads to

    σn+1˜un+12σn˜un+12+τΩρn+1un|˜un+1|2dx+τ2Ωρn+1(un)|˜un+1|2dx=0. (2.8)

    By taking the inner product of (2.1d) with 2τ˜un+1, we have

    σn˜un+12σnun2+σn(˜un+1un)22τλΩwn+1ϕn˜un+1dx+2τη˜un+12+τΩρn+1(un)|˜un+1|2dx+τ2Ωρn+1(un)|˜un+1|2dx=0. (2.9)

    Testing (2.1e) with 2τun+1 yields

    σn+1un+12σn+1˜un+12+σn+1(un+1˜un+1)2+ητun+12+ητ˜un+12+ητ(un+1˜un+1)2=0. (2.10)

    Taking into account the boundary condition on un and using integration by parts, we have

    Ωρn+1un|˜un+1|2dx+Ωρn+1(un)|˜un+1|2dx+Ωρn+1un|˜un+1|2dx=Ω(ρn+1un|˜un+1|2)dx=Ωρn+1un|˜un+1|2ndx=0. (2.11)

    Summing the above inequality, we arrive at

    2τλγwn+12+λ2ε2(ϕn+14L4ϕn4L4)+λ2ε2(ϕn+12ϕn2)2+λε2ϕn+1(ϕn+1ϕn)2+λε2(ϕn2ϕn+12+ϕn+1ϕn2)+λ(ϕn+12ϕn2+(ϕn+1ϕn)2)+σn+1un+12σnun2+σn(˜un+1un)2+σn+1(un+1˜un+1)2+ητun+12+ητ˜un+12+ητ(un+1˜un+1)2=0. (2.12)

    Adding up the above inequality from n=0 to N1, we obtain Theorem 2.3.

    The bound on the pressure p is proved in the following theorem.

    Theorem 3. (Stabilityofp) For any τ>0, the solution pn+1 to (2.1e) satisfies the following inequality:

    τ2N1n=0pn+12C(σ0u02+λϕ02λε2ϕ02+λ2ε2ϕ04L4)(ρ0+1). (2.13)

    Proof. Under the inf-sup condition, there exists a positive constant β such that

    βpn+1sup (2.14)

    Testing (2.14) with all v \in V leads to

    \begin{equation} \begin{aligned} \left(\nabla \cdot v, p^{n+1}\right) = & \eta\left(\nabla\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right), \nabla v\right)+\tau^{-1}\left(\rho^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right), v\right) \\ \leq & \eta\left(\left\|\nabla \boldsymbol{u}^{n+1}\right\|+\left\|\nabla \tilde{\boldsymbol{u}}^{n+1}\right\|\right)\|\nabla v\| \\ &+\tau^{-1}\left\|\sigma^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right)\right\|_{L^{2}}\left\|\sigma^{n+1}\right\|_{L^{3}}\|v\|_{L^{6}}. \end{aligned} \end{equation} (2.15)

    Given that \left\|\rho^{n}\right\| \leq\left\|\rho_{0}\right\| for all 1 \leq n \leq N , and by using Hölder's inequality, we have

    \begin{equation} \|\sigma^{n+1}\|_{L^{3}} = (\int{\|p^{n+1}\|^{\frac{3}{2}}})^{\frac{1}{3}} \leq C \| \rho^{n+1} \|^{\frac{1}{2}} \leq C \| \rho^{0} \|^{\frac{1}{2}}. \end{equation} (2.16)

    Then, by the Sobolev embedding inequality \|v\|_{L^{6}} \leq C\|\nabla v\|_{L^{2}} for any v \in V , we get

    \begin{equation} \begin{aligned} \left(\nabla \cdot v, p^{n+1}\right) \leq & \eta\left(\left\|\nabla \boldsymbol{u}^{n+1}\right\|+\left\|\nabla \tilde{\boldsymbol{u}}^{n+1}\right\|\right)\|\nabla v\| \\ &+C \tau^{-1}\left\|\sigma^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right)\right\|\left\|\rho_{0}\right\|^{\frac{1}{2}}\|\nabla v\|. \end{aligned} \end{equation} (2.17)

    Substituting the above inequalities into Eq (2.14), we obtain

    \begin{equation} \begin{aligned} \beta\left\|p^{n+1}\right\| \leq & \eta\left(\left\|\nabla \boldsymbol{u}^{n+1}\right\|+\left\|\nabla \tilde{\boldsymbol{u}}^{n+1}\right\|\right) \\ &+C \tau^{-1}\left\|\sigma^{n+1}\left(\boldsymbol{u}^{n+1}-\tilde{\boldsymbol{u}}^{n+1}\right)\right\|\left\|\rho_{0}\right\|^{\frac{1}{2}}. \end{aligned} \end{equation} (2.18)

    And by using Theorem 2.2, we get the desired result.

    In this section, we will give the time error estimates and show that the scheme has a first-order convergence rate. Although we verified that the scheme (2.1) is unconditionally stable in the previous chapter, we need to make the following assumptions[34,35] when conducting temporal error analysis:

    \begin{equation} \{\rho^n\}_{n = 0, \cdots, N} \ is\ uniformly\ bounded\ in \ L^\infty, \end{equation} (3.1)
    \begin{equation} for \ all\ n = 0, \cdots, N, \ {it\ holds\ that}\ \rho^n\geq\chi \ a.e. \ in \ \Omega , \end{equation} (3.2)

    where \chi is a number in (0, \rho_0^{min}] .

    We assume that the exact solution (\rho, \boldsymbol{u}, \phi, w, p) is sufficiently smooth. To be more precise,

    \begin{align} \rho &\in H^{2}\left(0, T ; L^{2}(\Omega)\right) \cap L^{2}\left(0, T ; W^{1, \infty}(\Omega)\right), \\ \phi &\in L^{\infty}\left(0, T ; H^{3}(\Omega)\right) \cap W^{1, \infty}\left(0, T ; H^{2}(\Omega)\right)\cap W^{2, \infty}\left(0, T ; H^{1}(\Omega)\right) \cap W^{3, \infty}\left(0, T ; L^{2}(\Omega)\right), \\ w& \in L^{\infty}\left(0, T ; H^{2}(\Omega)\right) \cap W^{1, \infty}\left(0, T ; L^{2}(\Omega)\right), \\ \boldsymbol{u} &\in \boldsymbol{H}^{2}\left(0, T ; \boldsymbol{L}^{2}(\Omega)\right) \cap \boldsymbol{L}^{\infty}\left(0, T ; \boldsymbol{V} \cap \boldsymbol{H}^{2}(\Omega)\right), \\ p& \in W^{2, \infty}\left(0, T ; H^{1}(\Omega)\right). \end{align} (3.3)

    We denote

    \begin{equation*} \begin{aligned} &e_{\phi}^{n} = \phi\left(t_{n}\right)-\phi^{n}, e_{w}^{n} = w\left(t_{n}\right)-w^{n}, e_{\rho}^{n} = \rho\left(t_{n}\right)-\rho^{n}, \\ &e_{u}^{n} = \boldsymbol{u} \left(t_{n}\right)-\boldsymbol{u} ^{n}, \tilde{e}_{u}^{n} = \boldsymbol{u} \left(t_{n}\right)-\tilde{\boldsymbol{u} }^{n}, e_{p}^{n} = p\left(t_{n}\right)-p^{n}.\\ \end{aligned} \end{equation*}

    In (1.1), taking t = t_{n+1} and subtracting from (2.1), we get the following error equations:

    \begin{align} &\frac{e_{\rho}^{n+1}-e_{\rho}^{n}}{\tau} +\nabla e_{\rho}^{n+1}\cdot\boldsymbol{u}(t_{n+1}) + \nabla\rho^{n+1}\cdot \left(\boldsymbol{u}(t_{n+1})-\boldsymbol{u}(t_{n})\right) +\nabla\rho^{n+1}\cdot e_{u}^{n} = R_\rho^{n+1}, \end{align} (3.4a)
    \begin{align} &\frac{e^{n+1}_\phi-e^{n}_\phi}{\tau} -\gamma\triangle e^{n+1}_w+\boldsymbol{u}(t_{n+1})\nabla\phi(t_{n+1})-\tilde{\boldsymbol{u}}^{n+1}\nabla\phi^{n} = R_\phi^{n+1}, \end{align} (3.4b)
    \begin{align} &e^{n+1}_w+\triangle e^{n+1}_\phi = \frac{1}{\varepsilon^2}(\phi(t_{n+1})^3-(\phi^{n+1})^3-\phi(t_{n+1})+\phi^n), \end{align} (3.4c)
    \begin{align} &\rho^{n} \frac{\tilde{e}_{u}^{n+1}-e_{u}^{n}}{\tau}-\eta \triangle \tilde{e}_{u}^{n+1}+\nabla p\left(t_{n+1}\right) +\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right) \\ &\qquad \quad -\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1} -\lambda w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)+ \lambda w^{n+1} \nabla \phi^{n} = R_{u}^{n+1}, \end{align} (3.4d)
    \begin{align} &\rho^{n+1} \frac{e_{u}^{n+1}-\tilde{e}_{u}^{n+1}}{\tau}-\eta\triangle \left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)-\nabla p^{n+1} = 0, \end{align} (3.4e)
    \begin{align} &\nabla \cdot e_{u}^{n+1} = 0, \end{align} (3.4f)

    where

    \begin{equation*} \begin{split} &R_{\phi}^{n+1} = \frac{\phi\left(t_{n+1}\right)-\phi\left(t_{n}\right)}{\tau}-\phi_{t}\left(t_{n+1}\right), \\ &R_{\rho}^{n+1} = \frac{\rho\left(t_{n+1}\right)-\rho\left(t_{n}\right)}{\tau}-\rho_{t}\left(t_{n+1}\right), \\ &R_{u}^{n+1} = \rho^{n} \frac{\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)}{\tau}-\rho\left(t_{n+1}\right) \boldsymbol{u}_{t}\left(t_{n+1}\right). \\ \end{split} \end{equation*}

    If the exact solution is sufficiently smooth, it is easy to establish the following estimate of the truncation error.

    Lemma 1. Under the regularity assumptions given by (3.3), the truncation errors satisfy:

    \begin{align*} \left\|R_{\rho}^{n+1}\right\|^{2}& \leq C \tau \int_{t_{n}}^{t_{n+1}}\left\|\rho_{t t}(t)\right\|^{2} d t \leq C \tau^{2}, \\ \left\|R_{\phi}^{n+1}\right\|^{2} &\leq C \tau \int_{t_{n}}^{t_{n+1}}\left\|\phi_{t t}(t)\right\|^{2} d t \leq C \tau^{2}, \\ \left\|R_{u}^{n+1}\right\|^{2} &\leq C \tau \int_{t_{n}}^{t_{n+1}}\left(\left\|\boldsymbol{u}_{t t}(t)\right\|^{2}+\left\|\rho_{t}(t)\right\|^{2}\right) d t+C\left\|e_{\rho}^{n}\right\|^{2} \\ &\leq C \tau^{2}+C\left\|e_{\rho}^{n}\right\|^{2}, \end{align*}

    for all 0 \leq n \leq N-1 .

    Proof. By using the integral residual of the Taylor formula, we have

    \begin{equation} \begin{split} R_\rho^{n+1} = \frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_n)\rho_{tt}(t)dt. \end{split} \end{equation} (3.5)

    By Hölder's inequality, we can derive

    \begin{equation} \begin{split} \Vert R_\rho^{n+1}\Vert^2 = & \int_\Omega(\frac{1}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_n)\rho_{tt}(t)dt)^2dx\\ \leq&\frac{1}{\tau^2}\int_\Omega\int_{t_{n}}^{t_{n+1}}(t-t_n)^2dt^{\frac{1}{2}\cdot2}\int_{t_{n}}^{t_{n+1}}\rho_{tt}(t)^2dt^{\frac{1}{2}\cdot2}dx\\ \leq&\frac{1}{\tau^2}(\frac{\tau^3}{3}\int_{t_{n}}^{t_{n+1}}\int_\Omega\rho_{tt}(t)^2dxdt)\\ \leq&C\tau\int_{t_{n}}^{t_{n+1}}\Vert \rho_{tt}(t)\Vert^2 dt\\ \leq&C\tau^2.\\ \end{split} \end{equation} (3.6)

    Similarly, we can prove the inequality of R_\phi^{n+1} . For R_u^{n+1} , we can rewrite

    \begin{equation} \begin{split} R_u^{n+1}& = \rho^n(\frac{\boldsymbol{u}(t_{n+1})-\boldsymbol{u}(t_{n})}{\tau}-\boldsymbol{u}_t(t_{n+1})) -(e_\rho^n+\int_{t_{n}}^{t_{n+1}}\rho_{t}(t)dt)\boldsymbol{u}_t(t_{n+1})\\ & = \frac{\rho^n}{\tau}\int_{t_{n}}^{t_{n+1}}(t-t_n)\boldsymbol{u}_{tt}(t)dt-(e_\rho^n+\int_{t_{n}}^{t_{n+1}}\rho_{t}(t)dt)\boldsymbol{u}_t(t_{n+1}).\\ \end{split} \end{equation} (3.7)

    Using R_\rho^{n+1} estimation and Hölder's inequality can yield the result for R_u^{n+1} .

    We introduce the following Gronwall's inequality, which will frequently be used in error estimates.

    Lemma 2. Let a_{k} , b_{k} , c_{k} and \gamma_{k} , for integers k\geq0 , be the nonnegative numbers such that

    \begin{equation*} a_{n}+\tau\sum^{n}\limits_{k = 0}b_{k}\leq\tau\sum^{n}\limits_{k = 0}\gamma_{k}a_{k}+\tau\sum^{n}\limits_{k = 0}c_{k}+B\ \ for\ \ \geq0. \end{equation*}

    Suppose that \tau\gamma_{k} < 1 , for all k, and set \sigma_{k} = (1-\tau\gamma_{k})^{-1} . Then,

    \begin{equation*} a_{n}+\tau\sum^{n}\limits_{k = 0}b_{k}\leq \exp\bigg(\tau\sum^{n}\limits_{k = 0}\gamma_{k}\sigma_{k}\bigg)\bigg(\tau\sum^{n}\limits_{k = 0}c_{k}+B\bigg)\ \ for\ \ \geq0. \end{equation*}

    We verify Lemma 5 by the following lemma.

    Lemma 3. Define

    \begin{equation} G_{c}^{n+1} = \left( \phi(t_{n+1})\right)^3-(\phi^{n+1})^3 = 3\left( \phi(t_{n+1})\right)^2 e_{\phi}^{n+1}-3\phi(t_{n+1})\left( e_{\phi}^{n+1}\right)^2+\left( e_{\phi}^{n+1}\right)^3. \end{equation} (3.8)

    Then, for n < \frac{T}{\tau}-1 , we have

    \begin{equation} \begin{split} &\left\|G_{c}^{n+1}\right\| \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{1}}, \\ &\left\|e_{\phi}^{n+1}\right\|_{H^{2}}\leq C\left(\tau+\left\| e_w^{n+1} \right\|+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}+\left\| e_\phi^{n} \right\| \right).\\ \end{split} \end{equation} (3.9)

    Proof. For \left\|G_{c}^{n+1}\right\| , we use (3.8) to conclude that

    \begin{equation} \begin{split} \left\|G_{c}^{n+1}\right\| & \leq C\left(\left\|e_{\phi}^{n+1}\right\|\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}^{2}+\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{2}\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}+\left\|e_{\phi}^{n+1}\right\|_{L^{6}}^{3}\right) \\ & \leq C\left(\left\|e_{\phi}^{n+1}\right\|\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}^{2}+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}^{2}\left\|\phi\left(t^{n+1}\right)\right\|_{L^{\infty}}+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}^{3}\right) \\ & \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{1}}, \end{split} \end{equation} (3.10)

    where we have used the a priori bound \left\|\phi^{n}\right\|_{H^{1}} \leq C implied by the stability result given by Theorem 2. Using the H^{2} regularity results for elliptic equations, we conclude that

    \left\|e_{\phi}^{n+1}\right\|_{H^{2}} \leq c\left(\left\|e_{\phi}^{n+1}\right\|_{L^{2}}+\left\|\triangle e_{\phi}^{n+1}\right\|_{L^{2}}\right){;}

    from (3.4c), we know that

    \begin{equation*} \begin{split} \triangle e_{\phi}^{n+1} = -e_w^{n+1}+\frac{1}{\varepsilon^2}G_c^{n+1}-\frac{1}{\varepsilon^2}e_{\phi}^{n} -\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt{;} \end{split} \end{equation*}

    thus,

    \begin{equation*} \begin{split} \left\|e_{\phi}^{n+1}\right\|_{H^{2}}\leq& C\left( \left\|e_{\phi}^{n+1}\right\|+\left\| e_w^{n+1} \right\| +\left\|G_c^{n+1} \right\|+\left\| e_\phi^{n} \right\|+\tau\right)\\ \leq& C\left(\tau+\left\| e_w^{n+1} \right\|+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}+\left\| e_\phi^{n} \right\| \right).\\ \end{split} \end{equation*}

    The error estimate for the discrete density \rho^{n+1} is derived in the following lemma.

    Lemma 4. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) hold. Then, we have

    \begin{equation} \left\|e_{\rho}^{n+1}\right\|^{2}+2 \sum\limits_{m = 0}^{n}\left\|e_{\rho}^{m+1}-e_{\rho}^{m}\right\|^{2} \leq C\left(\tau^{2}+\tau \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}\right) \end{equation} (3.11)

    for all 0 \leq n \leq N-1 .

    Proof. Multiplying (3.4a) by 2 \tau e_{\rho}^{n+1} and integrating over \Omega , we have

    \begin{equation} \begin{split} &\left\|e_{\rho}^{n+1}\right\|^{2}-\left\|e_{\rho}^{n}\right\|^{2}+\left\|e_{\rho}^{n+1}-e_{\rho}^{n}\right\|^{2}\\ = &-2 \tau\left(\nabla e_{\rho}^{n+1} \cdot \boldsymbol{u}\left(t_{n+1}\right), e_{\rho}^{n+1}\right)+2 \tau\left(R_{\rho}^{n+1}, e_{\rho}^{n+1}\right)\\ &-2 \tau\left(\nabla \rho^{n+1} \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right)-2 \tau\left(\nabla \rho^{n+1} \cdot e_{u}^{n}, e_{\rho}^{n+1}\right)\\ = &\sum\limits_{i = 1}^{4} K_{i}.\\ \end{split} \end{equation} (3.12)

    Using \nabla \cdot \boldsymbol{u}\left(t_{n+1}\right) = 0 in \Omega and \boldsymbol{u}\left(t_{n+1}\right) = 0 on \partial \Omega , we have

    \begin{equation} K_{1} = \frac{1}{2} \int_{\partial \Omega}\vert e_{\rho}^{n+1}\vert^{2} \boldsymbol{u}\left(t_{n+1}\right) \cdot \boldsymbol{n} d s = 0 . \end{equation} (3.13)

    By using Young's inequality, the Cauchy-Schwarz inequality and Lemma 1, we have

    K_{2} \leq 2 \tau\left\|R_{\rho}^{n+1}\right\|\left\|e_{\rho}^{n+1}\right\| \leq C \tau\left\|R_{\rho}^{n+1}\right\|^{2}+\varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2} \leq C \tau^{3}+\varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2},

    where \varepsilon > 0 is a sufficiently small constant.

    Then, by the Sobolev inequality and Young inequality, we get

    \begin{equation} \begin{aligned} K_{3} & = -2 \tau\left(\nabla \rho\left(t_{n+1}\right) \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right)+2 \tau\left(\nabla e_{\rho}^{n+1} \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right) \\ & = -2 \tau\left(\nabla \rho\left(t_{n+1}\right) \cdot\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right), e_{\rho}^{n+1}\right) \\ & \leq 2 \tau\left\|\nabla \rho\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}} \boldsymbol{u}_{t}(t) d t\right\|\left\|e_{\rho}^{n+1}\right\| \\ & \leq \varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2}+C \tau^{2} \int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{t}(t)\right\|^{2} d t \\ & \leq C \tau^{3}+\varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2}. \end{aligned} \end{equation} (3.14)

    Similarly, the last term can be estimated as follows:

    \begin{equation} \begin{aligned} K_{4} & \leq 2 \tau\left\|\nabla \rho\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\frac{1}{\sigma^{n}}\right\|_{L^{\infty}}\|\| \sigma^{n} e_{u}^{n}\|\| e_{\rho}^{n+1} \| \\ & \leq \varepsilon \tau\left\|e_{\rho}^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}. \end{aligned} \end{equation} (3.15)

    If \varepsilon is sufficiently small such that \varepsilon \tau \leq \frac{1}{6} , substituting the estimates of K_{i}\ (1 \leq i \leq 4) into (3.12), we have

    \begin{equation} \left\|e_{\rho}^{n+1}\right\|^{2}-\left\|e_{\rho}^{n}\right\|^{2}+2\left\|e_{\rho}^{n+1}-e_{\rho}^{n}\right\|^{2} \leq C \tau^{3}+C \tau\left\|e_{\rho}^{n}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}. \end{equation} (3.16)

    Using the discrete Gronwall inequality, we obtain the desired result.

    Lemma 5. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) are valid. For sufficiently small \tau , there are the following error estimates:

    \begin{equation} \begin{split} &\tau\gamma\sum\limits_{n = 0}^{N-1}\left(\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+\left\| e_{w}^{n+1} \right\|^2 \right)+\left\| e_{\phi}^{N}\right\|^2+\lambda\left\| \nabla e_{\phi}^{N}\right\|^2+\frac{\lambda}{2 \varepsilon^{2}}\left\|e_{\phi}^{N}\right\|_{L^{4}}^{4}\\ &+\left\|\sigma^{N} e_{u}^{N}\right\|^{2}+ \sum\limits_{n = 0}^{N-1}\left(\frac{1}{2}\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\right)\\ &+\tau\eta\sum\limits_{n = 0}^{N-1}\left(\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \frac{1}{2}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \right)\\ \leq& C\tau.\\ \end{split} \end{equation} (3.17)

    Proof. Let us multiply (3.4a) by \tau\vert\tilde{e}_{u}^{n+1}\vert^{2} , (3.4b) by 2 \tau e_{\phi}^{n+1} and 2 \tau \lambda e_{w}^{n+1} , (3.4c) by 2\tau\gamma e_{w}^{n+1} and -2 \lambda(e_{\phi}^{n+1}-e_{\phi}^{n}) , (3.4d) by 2 \tau \tilde{e}_{u}^{n+1} and (3.4e) by 2 \tau e_{u}^{n+1} . Summing up all of the above equations, we have

    \begin{align} &2\tau\gamma\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+ \left\| e_{\phi}^{n+1}\right\|^2- \left\| e_{\phi}^{n}\right\|^2+ \left\| e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &+2\tau\gamma\left\| e_{w}^{n+1} \right\|^2+\lambda\left\| \nabla e_{\phi}^{n+1}\right\|^2- \lambda\left\| \nabla e_{\phi}^{n}\right\|^2+ \lambda\left\| \nabla e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\\ &+ \tau\eta\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \tau\eta\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \tau\eta\left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \\ = &-2 \tau\left(\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) +2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1}, \tilde{e}_{u}^{n+1}\right)\\ &-\tau\left(\nabla \rho^{n+1} \cdot {\bf u}^{n}, \vert e_{u}^{n+1}\vert^{2}\right)-2 \tau\left(\nabla p\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right)+2 \tau\left(R_{u}^{n+1}, \tilde{e}_{u}^{n+1}\right) \\ &+2\tau\lambda \left( w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- w^{n+1} \nabla \phi^{n}, \tilde{e}_{u}^{n+1}\right) \\ &-2\tau\lambda \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)+ 2 \tau\lambda\left(R_{\phi}^{n+1}, e_{w}^{n+1}\right)\\ &-2\tau \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right) + 2 \tau\left(R_{\phi}^{n+1}, e_{\phi}^{n+1}\right)\\ &+\frac{2\tau\gamma}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)^3-\left(\phi^{n+1}\right)^3, e_{w}^{n+1}\right) -\frac{2\tau\gamma}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)-\phi^n, e_{w}^{n+1}\right)\\ &-\frac{2\lambda}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)^3-\left(\phi^{n+1}\right)^3, e_{\phi}^{n+1}-e_{\phi}^{n}\right) +\frac{2\lambda}{\varepsilon^2} \left(\phi\left(t_{n+1}\right)-\phi^n, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ = &\sum\limits_{i = 1}^{i = 14}A_i.\\ \end{align} (3.18)

    Thanks to \nabla \cdot {\bf u}^{n} = 0 in \Omega , we have

    \begin{equation} \begin{aligned} A_{2}+A_{3} & = 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1}, \tilde{e}_{u}^{n+1}\right)-\tau\left(\nabla \rho^{n+1} \cdot \boldsymbol{u}^{n}, \vert\tilde{e}_{u}^{n+1}\vert^{2}\right) \\ & = 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \widetilde{e}_{u}^{n+1}\right)-\tau\left(\nabla \cdot\left(\boldsymbol{u}^{n} \rho^{n+1}\vert\tilde{e}_{u}^{n+1}\vert^{2}\right), 1\right) \\ & = 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right). \end{aligned} \end{equation} (3.19)

    Therefore, we get

    \begin{align} &A_{1}+A_{2}+A_{3} \end{align} (3.20)
    \begin{align} = &2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) \\ = &-2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right)-2 \tau\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(\rho\left(t_{n+1}\right)\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), \tilde{e}_{u}^{n+1}\right) \\ \leq& C \tau\left\|\rho^{n+1}\right\|_{L^{\infty}}\left\|\frac{1}{\sigma^{n}}\right\|_{L^{\infty}}\left\|\sigma^{n} e_{u}^{n}\right\|\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|\tilde{e}_{u}^{n+1}\right\|_{L^{6}} \end{align} (3.21)
    \begin{align} &+C \tau\left\|e_{\rho}^{n+1}\right\|\left\|\boldsymbol{u}\left(t_{n}\right)\right\|_{L^{\infty}} \left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|\tilde{e}_{u}^{n+1}\right\|_{L^{6}}\\ &+C \tau\left\|\rho\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}} \boldsymbol{u}_{t}(t) d t\right\|\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|\tilde{e}_{u}^{n+1}\right\|_{L^{6}} \\ \leq& C \tau\left(\left\|\sigma^{n} e_{u}^{n}\right\|+\left\|e_{\rho}^{n+1}\right\|+\tau \int_{t_{n}}^{t_{n+1}}\left\|\boldsymbol{u}_{t}(t)\right\| d t\right)\left\|\nabla \tilde{e}_{u}^{n+1}\right\| \\ \leq& \frac{\eta \tau}{12}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+C \tau\left(\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|e_{\rho}^{n+1}\right\|^{2}+\tau^{2}\right) \\ \leq& C \tau^{3}+\frac{\eta \tau}{12}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2} , \\ \end{align} (3.22)

    where the Sobolev embedding inequality \|v\|_{L^{6}} \leq C\|\nabla v\|_{L^{2}} for any v \in V is used.

    For A_4 , we have

    \begin{equation} \begin{split} A_4&\leq2\tau \left\|\nabla p\left(t_{n+1}\right)\right\|\left\|\frac{1}{\sigma^n}\right\|_{L^\infty} \left\|\sigma^n\left(\tilde e_u^{n+1}-e_u^{n} \right)\right\|\\ &\leq C\tau^2+\frac{1}{2} \left\| \sigma^n\left(\tilde e_u^{n+1}-e_u^{n} \right) \right\|^2 .\\ \end{split} \end{equation} (3.23)

    Using the Poincar\acute{e} inequality and Cauchy-Schwarz inequality, we can deal with A_6 , A_7 and A_9 :

    \begin{align} A_6 = &2\tau\lambda \left( w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- w^{n+1} \nabla \phi^{n}, \tilde{e}_{u}^{n+1}\right)\\ = &2\tau\lambda \left( w\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right), \tilde{e}_{u}^{n+1}\right) +2\tau\lambda \left( w\left(t_{n+1}\right)\nabla e_\phi^n, \tilde{e}_{u}^{n+1}\right)\\ &+2\tau\lambda \left( \nabla\phi^n e_w^{n+1}, \tilde{e}_{u}^{n+1}\right)\\ \leq& 2\tau\lambda \left\|w\left(t_{n+1}\right)\right\|_{L^\infty} \left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_tdt\right\|\left\|\tilde{e}_{u}^{n+1}\right\| +2\tau\lambda \left\|w\left(t_{n+1}\right)\right\|_{L^\infty} \left\| \nabla e_\phi^n \right\| \left\|\tilde{e}_{u}^{n+1} \right\|\\ &+2\tau\lambda \left( \nabla\phi^n e_w^{n+1}, \tilde{e}_{u}^{n+1}\right)\\ \leq& C\tau^3+\frac{\eta\tau}{12}\left\|\nabla\tilde{e}_{u}^{n+1}\right\|^2+C\tau\left\| \nabla e_\phi^n \right\|^2+2\tau\lambda \left( \nabla\phi^n e_w^{n+1}, \tilde{e}_{u}^{n+1}\right), \\ A_7 = &-2\tau\lambda \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)\\ = &-2\tau\lambda \left( \boldsymbol{u}\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right) +\boldsymbol{u}\left(t_{n+1}\right)\nabla e_\phi^n+\nabla\phi^n\tilde{e}_{u}^{n+1}, e_{w}^{n+1}\right)\\ \leq& 2\tau\lambda\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_tdt\right\| \left\|e_{w}^{n+1}\right\| +2\tau\lambda \left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty} \left\| \nabla e_\phi^n \right\| \left\|e_{w}^{n+1} \right\|\\ & -2\tau\lambda \left(\nabla\phi^n\tilde{e}_{u}^{n+1}, e_{w}^{n+1}\right)\\ \leq& C\tau^3+\frac{\tau\gamma}{3} \left\|e_w^{n+1}\right\|^2+C\tau \left\|\nabla e_\phi^{n}\right\|^2 -2\tau\lambda \left(\nabla\phi^n\tilde{e}_{u}^{n+1}, e_{w}^{n+1}\right), \end{align} (3.24)

    and

    \begin{equation} \begin{split} A_9 = &-2\tau \left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)- \tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right)\\ \leq& 2\tau \left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_tdt\right\| \left\|e_{\phi}^{n+1}\right\| +2\tau \left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^\infty} \left\| \nabla e_\phi^n \right\| \left\|e_{\phi}^{n+1} \right\|\\ & +2\tau\left\|\nabla\phi^n \right\|\left\| \tilde{e}_{u}^{n+1} \right\|_{L_6} \left\| e_{\phi}^{n+1} \right\|_{L_6}\\ \leq& C\tau^3+C\tau \left\|\nabla e_\phi^{n+1}\right\|^2+C\tau \left\|\nabla e_\phi^{n}\right\|^2 +\frac{\eta\tau}{12} \left\|\nabla\tilde{e}_{u}^{n+1}\right\|^2.\\ \end{split} \end{equation} (3.25)

    From Lemma 1, we have

    \begin{equation} \begin{split} A_5+A_8+A_{10} \leq& 2\tau \left\|R_u^{n+1} \right\| \left\| \tilde{e}_{u}^{n+1} \right\| +2\tau\lambda\left\|R_\phi^{n+1} \right\| \left\| e_w^{n+1} \right\| +2\tau \left\| R_\phi^{n+1} \right\| \left\| e_\phi^{n+1} \right\|\\ \leq& C\tau^3+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2} +\frac{\eta\tau}{12} \left\|\nabla \tilde{e}_{u}^{n+1} \right\|^2+\frac{\tau\gamma}{3} \left\| e_w^{n+1} \right\|^2 +\frac{1}{2} \left\| e_\phi^{n+1} \right\|^2.\\ \end{split} \end{equation} (3.26)

    From Lemma 3, we obtain

    \begin{equation} \begin{aligned} A_{11}+A_{12} & = \frac{2\tau\gamma}{\varepsilon^2}\left(G_{c}^{n+1}-e_{\phi}^n-\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{w}^{n+1}\right) \\ & \leq \frac{2\tau\gamma}{\varepsilon^2}\left(\left\|e_{\phi}^n\right\|+\left\|G_{c}^{n+1}\right\|+ \left\|\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt\right\|\right)\left\|e_{w}^{n+1}\right\| \\ & \leq C\tau^3+C \tau\left(\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}\right) +\frac{\tau\gamma}{3}\left\|e_{w}^{n+1}\right\|^{2}, \\ \end{aligned} \end{equation} (3.27)

    and

    \begin{align} A_{13}+A_{14} = &- \frac{2\lambda}{\varepsilon^2}\left(G_{c}^{n+1}-e_{\phi}^n-\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{\phi}^{n+1}-e_{\phi}^{n}\right) \\ = &-\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|\left(e_{\phi}^{n+1}\right)^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e_{\phi}^{n}\right)\right\|^{2}\right)-\frac{2\lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2} \\ &-\frac{2\lambda}{\varepsilon^{2}}\left(3\left(\phi\left(t_{n+1}\right)\right)^{2} e_{\phi}^{n+1}-3 \phi\left(t_{n+1}\right)\left(e_{\phi}^{n+1}\right)^{2}-e_{\phi}^{n+1}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ &+\frac{2\lambda}{\varepsilon^{2}} \left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{\phi}^{n+1}-e_{\phi}^{n}\right), \end{align} (3.28)

    where we use this identity

    \begin{equation} \begin{split} a^{3}(a-b) = \frac{1}{4}\left(a^{4}-b^{4}+\left(a^{2}-b^{2}\right)^{2}+2 a^{2}(a-b)^{2}\right). \end{split} \end{equation} (3.29)

    We denote

    \begin{equation} \tilde{G}^{n+1} = -\frac{2\lambda}{\varepsilon^{2}}\left(3\left(\phi\left(t_{n+1}\right)\right)^{2} e_{\phi}^{n+1}-3 \phi\left(t_{n+1}\right)\left(e_{\phi}^{n+1}\right)^{2}-e_{\phi}^{n+1}\right) \text {. } \end{equation} (3.30)

    Similar to the method estimated from Lemma 3, we can obtain

    \begin{equation} \left\|\tilde{G}^{n+1}\right\| \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{1}} . \end{equation} (3.31)

    Taking the gradient of \tilde{G}^{n+1} , we get

    \begin{equation} \begin{aligned} \nabla \tilde{G}^{n+1} = &-\frac{2\lambda}{\varepsilon^{2}}\left[\left(3\left(\phi\left(t_{n+1}\right)\right)^{2}-1\right) \nabla e_{\phi}^{n+1}+6 \phi\left(t_{n+1}\right) e_{\phi}^{n+1} \nabla \phi\left(t_{n+1}\right)\right.\\ &\left.-3\left(e_{\phi}^{n+1}\right)^{2} \nabla \phi\left(t_{n+1}\right)-6 \phi\left(t_{n+1}\right) e_{\phi}^{n+1} \nabla e_{\phi}^{n+1}\right]. \end{aligned} \end{equation} (3.32)

    Since H^{2}(\Omega) \subset L^{\infty}(\Omega) and by utilizing the bound of \left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}} implied by Theorem 2, we conclude that

    \begin{equation} \left\|e_{\phi}^{n+1} \nabla e_{\phi}^{n+1}\right\|_{L^{2}} \leq\left\|e_{\phi}^{n+1}\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}} \leq C\left\|e_{\phi}^{n+1}\right\|_{H^{2}} . \end{equation} (3.33)

    In view of (3.10) and the bound \left\|\phi^{n}\right\|_{H^{1}} < C , we have

    \begin{equation} \begin{array}{rl} \left\|\nabla \tilde{G}^{n+1}\right\| \leq & C\left[\left(\left\|\phi\left(t_{n+1}\right)\right\|_{L^{\infty}}^{2}+1\right)\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}}\right.\\ &+\left\|\phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{3}}\left\|e_{\phi}^{n+1}\right\|_{L^{6}} \\ & \left.+\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{6}}\left\|e_{\phi}^{n+1}\right\|_{L^{6}}^{2}+\left\|\phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|e_{\phi}^{n+1} \nabla e_{\phi}^{n+1}\right\|_{L^{2}}\right] \\ \leq &C\left(\left\|\nabla e_{\phi}^{n+1}\right\|_{L^{2}}+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}\right. \left.+\left\|e_{\phi}^{n+1}\right\|_{H^{1}}^{2}+\left\| e_{\phi}^{n+1} \nabla e_{\phi}^{n+1} \right\|_{L^{2}}\right) \\ \leq &C\left(\tau+\left\|e_{w}^{n+1}\right\|+\left\|e_{\phi}^{n+1}\right\|+\left\|\nabla e_{\phi}^{n+1}\right\| +\left\|e_{\phi}^{n}\right\| \right) . \end{array} \end{equation} (3.34)

    Dealing with the penultimate term of (3.28) and in view of (3.31) and (3.33), we get

    \begin{equation} \begin{split} &\left(\tilde{G}^{n+1}, e_{\phi}^{n+1}-e_{\phi}^{n}\right) \\ = &\tau\left(\tilde{G}^{n+1}, \gamma \Delta e_{w}^{n+1}-\boldsymbol{u}\left(t_{n+1}\right) \nabla e_{\phi}^{n}-\tilde{e}_{u}^{n+1} \nabla \phi^{n}\right)\\ &+\tau\left(\tilde{G}^{n+1}, -\boldsymbol{u}\left(t_{n+1}\right)\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt+R_{\phi}^{n+1}\right)\\ \leq& \gamma\tau\left\|\nabla e_{w}^{n+1}\right\|\left\|\nabla \tilde{G}^{n+1}\right\| +\tau\left\| \nabla\phi^n\right\|\left\| \tilde{e}_u^{n+1}\right\|_{H_1}\left\| \tilde{G}^{n+1} \right\|_{H_1}\\ &+\tau\left\|\tilde{G}^{n+1}\right\|\left(\left\|R_{\phi}^{n+1}\right\|+\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\|+\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt\right\| \right) \\ \leq& C\tau^3+\frac{\gamma\lambda\tau}{2}\left\|\nabla e_{w}^{n+1}\right\|^2+\frac{\eta\tau}{12}\left\| \nabla\tilde{e}_u^{n+1}\right\|^2\\ &+C\tau\left(\left\|e_{w}^{n+1}\right\|^2+\left\|e_{\phi}^{n+1}\right\|^2+\left\|\nabla e_{\phi}^{n+1}\right\|^2 +\left\|e_{\phi}^{n}\right\| ^2+\left\|\nabla e_{\phi}^{n}\right\| ^2 \right).\\ \end{split} \end{equation} (3.35)

    Then, we estimate the last term of (3.28),

    \begin{align} & \frac{2\lambda}{\varepsilon^2}\left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, e_{\phi}^{n+1}-e_{\phi}^{n}\right) \\ = &\frac{2\lambda\tau}{\varepsilon^2}\left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, \gamma \Delta e_{w}^{n+1}-\boldsymbol{u}\left(t_{n+1}\right) \nabla e_{\phi}^{n}-\tilde{e}_{u}^{n+1} \nabla \phi^{n}\right)\\ &+\frac{2\lambda\tau}{\varepsilon^2}\left(\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt, -\boldsymbol{u}\left(t_{n+1}\right)\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt+R_{\phi}^{n+1}\right)\\ \leq& \frac{2\lambda\gamma\tau}{\varepsilon^2}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt\right\|\left\|\nabla e_w^{n+1}\right\| +\frac{2\lambda\tau}{\varepsilon^2}\left\| \nabla\phi^n\right\|\left\| \tilde{e}_u^{n+1}\right\|_{H_1}\left\| \int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt \right\|_{H_1}\\ &+\frac{2\lambda\tau}{\varepsilon^2}\left\|\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt\right\|\left(\left\|R_{\phi}^{n+1}\right\|+\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\| \right) \\ &+\frac{2\lambda\tau}{\varepsilon^2}\left\|\int_{t_{n}}^{t_{n+1}}\phi_t\left(t\right)dt\right\|\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}}\nabla\phi_t\left(t\right)dt\right\| \\ \leq& C\tau^3+\frac{\gamma\lambda\tau}{2}\left\|\nabla e_{w}^{n+1}\right\|^2+\frac{\eta\tau}{12}\left\|\nabla \tilde{e}_{u}^{n+1}\right\| ^2 +C\tau\left\|\nabla e_{\phi}^{n}\right\| ^2 .\\ \end{align} (3.36)

    Combining (3.35) and (3.36), we can obtain

    \begin{equation} \begin{split} A_{13}+A_{14} \leq&-\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|\left(e_{\phi}^{n+1}\right)^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e_{\phi}^{n}\right)\right\|^{2}\right)-\frac{2\lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2} \\ &+ C\tau^3+\gamma\lambda\tau\left\|\nabla e_{w}^{n+1}\right\|^2+\frac{\eta\tau}{6}\left\|\nabla \tilde{e}_{u}^{n+1}\right\| ^2\\ &+C\tau \left( \left\|\nabla e_{\phi}^{n}\right\| ^2+\left\|e_{w}^{n+1}\right\|^2+\left\|e_{\phi}^{n+1}\right\|^2+\left\|\nabla e_{\phi}^{n+1}\right\|^2 +\left\|e_{\phi}^{n}\right\| ^2+\left\|\nabla e_{\phi}^{n}\right\| ^2 \right) .\\ \end{split} \end{equation} (3.37)

    Substituting the above estimates into (3.18), we have

    \begin{equation} \begin{split} &\tau\gamma\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+ \frac{1}{2}\left\| e_{\phi}^{n+1}\right\|^2- \left\| e_{\phi}^{n}\right\|^2+ (1+\frac{2\lambda}{\varepsilon^{2}})\left\| e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &+\tau\gamma\left\| e_{w}^{n+1} \right\|^2+\lambda\left\| \nabla e_{\phi}^{n+1}\right\|^2- \lambda\left\| \nabla e_{\phi}^{n}\right\|^2+ \lambda\left\| \nabla e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^2\\ &+\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\frac{1}{2}\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\\ &+ \tau\eta\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \frac{\tau\eta}{2}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \tau\eta\left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \\ &+\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|\left(e_{\phi}^{n+1}\right)^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e_{\phi}^{n}\right)\right\|^{2}\right)\\ \leq& C \tau^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}\\ &+C\tau \left( \left\|\nabla e_{\phi}^{n}\right\| ^2+\left\|e_{w}^{n+1}\right\|^2+\left\|e_{\phi}^{n+1}\right\|^2+\left\|\nabla e_{\phi}^{n+1}\right\|^2 +\left\|e_{\phi}^{n}\right\| ^2+\left\|\nabla e_{\phi}^{n}\right\| ^2 \right).\\ \end{split} \end{equation} (3.38)

    Adding up from 0 to N-1 , and applying Gronwall's inequality, we infer that

    \begin{align*} \label{eq12} &\tau\gamma\sum\limits_{n = 0}^{N-1}\left(\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+\left\| e_{w}^{n+1} \right\|^2 \right)+\left\| e_{\phi}^{N}\right\|^2+\lambda\left\| \nabla e_{\phi}^{N}\right\|^2+\frac{\lambda}{2 \varepsilon^{2}}\left\|e_{\phi}^{N}\right\|_{L^{4}}^{4}\\ &+\left\|\sigma^{N} e_{u}^{N}\right\|^{2}+ \sum\limits_{n = 0}^{N-1}\left(\frac{1}{2}\left\|\sigma^{n}\left(\tilde{e}_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2}\right)\\ &+\tau\eta\sum\limits_{n = 0}^{N-1}\left(\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \frac{1}{2}\left\|\nabla \tilde{e}_{u}^{n+1}\right\|^{2}+ \left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|^{2} \right)\\ \leq& C\tau. \end{align*}

    Lemma 5 shows that the fractional scheme converges at a rate of O(\tau^{\frac{1}{2}}) . However, since we use the first-order backward Euler method of time discretization, it is not optimal from the perspective of theoretical analysis. Next, we will improve the convergence speed to the first order.

    Lemma 6. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) are valid. For sufficiently small \tau , there are the following error estimates:

    \begin{equation} \begin{split} &\tau\gamma\sum\limits_{n = 0}^{N-1}\left(\lambda\left\| \nabla e_{w}^{n+1} \right\|^2+\left\| e_{w}^{n+1} \right\|^2 \right)+\left\| e_{\phi}^{N}\right\|^2+\lambda\left\| \nabla e_{\phi}^{N}\right\|^2\\ &+\frac{\lambda}{2 \varepsilon^{2}}\left\|e_{\phi}^{N}\right\|_{L^{4}}^{4} +\left\|\sigma^{N} e_{u}^{N}\right\|^{2}+\tau\eta\sum\limits_{n = 0}^{N-1}\left\|\nabla e_{u}^{n+1}\right\|^{2}\\ \leq& C\tau^2.\\ \end{split} \end{equation} (3.39)

    Proof. Taking the sum of (3.4d) and (3.4e), we get

    \begin{equation} \begin{split} &\rho^{n} \frac{e_u^{n+1}-e_u^ n}{\tau}+\frac{\rho^{n+1}-\rho^{n}}{\tau}\left(e_{u}^{n+1}-\tilde{e}_{n}^{n+1}\right)-\eta \triangle e_{u}^{n+1}+\rho\left(t_{n+1}\right)\left(u\left(t_{n+1}\right) \cdot \nabla\right) u(t_{n+1})\\ &-\rho^{n+1}\left(u^{n} \cdot \nabla\right) \tilde{u}^{n+1}+\Delta e_{p}^{n+1}-\lambda w\left(t _{n+1}\right) \nabla \phi(t_{n+1})+\lambda w^{n+1}\nabla\phi^{n} = R_{u}^{n+1}.\\ \end{split} \end{equation} (3.40)

    Let us multiply (3.4a) by \tau\tilde{e}_{u}^{n+1} , (3.4b) by 2 \tau e_{\phi}^{n+1} and 2 \tau \lambda e_{w}^{n+1} , (3.4c) by 2\tau\gamma e_{w}^{n+1} and -2 \lambda(e_{\phi}^{n+1}-e_{\phi}^{n}) and (3.40) by 2 \tau e_{u}^{n+1} . Summing up all of the above equations, we obtain

    \begin{align} &\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|\sigma^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}+ 2 \eta\tau\left\|\nabla e_{u}^{n+1}\right\|^{2}+2 \tau \gamma \lambda\left\|\nabla e_{w}^{n+1}\right\|^{2}\\ &+\left\|e_{\phi}^{n+1}\right\|^{2}-\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2}+2 \tau \gamma\left\|e_{w}^{n+1}\right\|^{2}+\lambda\left\|\nabla e_{\phi}^{n+1}\right\|^{2}-\lambda\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+\lambda\left\|\nabla e _{\phi}^{n+1}-\nabla e_{\phi}^{n} \right\|^{2}\\ = &-2\tau\left(\nabla\rho^{n+1}\cdot u^n, \tilde{e}_{u}^{n+1}\cdot e_{u}^{n+1} \right)+\tau\left( \nabla\rho^{n+1}\cdot u^n, \vert e_{u}^{n+1}\vert^2 \right)\\ &-2 \tau\left(\rho\left(t_{n+1}\right)\left(u\left(t_{n+1}\right)\cdot\nabla\right) u\left(t_{n+1}\right)-\rho^{n+1}\left(u^{n}\cdot\nabla\right) \tilde{u}^{n+1}, e_u^{n+1}\right)\\ &+2 \tau \lambda\left(w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-w^{n+1} \nabla \phi^{n}, e_{u}^{n+1}\right)+2 \tau\left(R_{u}^{n+1}, e_u^{n+1}\right)\\ &-2 \tau \lambda\left(u\left(t_{n+1}\right) \nabla \phi\left(t _{n+1}\right)-\tilde{u}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)+2 \tau \lambda \left(R_{\phi}^{n+1}, e_w^{n+1}\right)\\ &-2 \tau\left(u\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-\tilde{u}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right)+2 \tau\left(R_{\phi}^{n+1}, e_{\phi}^{n+1}\right)\\ &+\frac{2 \tau r}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)^{3}-\left(\phi^{n+1}\right)^{3}, e_w^{n+1}\right)-\frac{2 \tau r}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)-\phi^{n}, e_{w}^{n+1}\right)\\ &-\frac{2 \lambda}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)^{3}-\left(\phi^{n+1}\right)^{3}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)+\frac{2 \lambda}{\varepsilon^{2}}\left(\phi\left(t_{n+1}\right)-\phi^{n}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ = &\sum\limits_{i = 1}^{13} L_{i}. \end{align} (3.41)

    Using integration by parts, we have

    \begin{equation} \begin{split} L_{1}+L_{2} = & 2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{e}_{u}^{n+1}, e_{ u}^{n+1}\right)+2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}\left(t_{n}\right) \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right){;} \end{split} \end{equation} (3.42)

    we rewrite the L_{3} as

    \begin{equation} \begin{split} L_{3} = &-2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde e_{u}^{n+1}, e_{u}^{n+1}\right)-2 \tau\left(\rho^{n+1}\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right)-2 \tau\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right){;} \end{split} \end{equation} (3.43)

    then,

    \begin{equation} \begin{split} L_{1}+L_{2}+L_{3} = & 2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(\boldsymbol{u}\left(t_{n}\right) \cdot \nabla\right) e_{u}^{n+1}, e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ &-2 \tau\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), e_{u}^{n+1}\right) \\ = & \sum\limits_{i = 1}^{5} J_{i}. \\ \end{split} \end{equation} (3.44)

    Lemma 5 shows that

    \begin{equation} \begin{split} \left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\| \leq C.\quad\quad \end{split} \end{equation} (3.45)

    According to Young's inequality and the Cauchy-Schwarz inequality, we can deduce that

    \begin{align} J_{1} & \leq \frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{u}^{n}\right\|^{2}\left\|\nabla\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\|\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\| \\ & \leq \frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau^{\frac{3}{2}}\left\|\nabla e_{u}^{n}\right\|^{2}, \\ J_{2} & \leq 2 \tau\left\|e_{\rho}^{n+1}\right\|\left\|\boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{3}}\left\|e_{u}^{n+1}\right\|_{L^{6}} \\ & \leq \frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|e_{\rho}^{n+1}\right\|^{2} \\ & \leq C \tau^{3}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}, \\ J_{3} & \leq C \tau\left\|\nabla e_{u}^{n+1}\right\|\left\|\boldsymbol{u}\left(t_{n}\right)\right\|_{{L}^{\infty}}\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\tilde{e}_{u}^{n+1}\right)\right\| \\ & \leq C \tau^{3}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2, \\ J_{4} & \leq 2 \tau\left\|\rho^{n+1}\right\|_{{L}^{\infty}}\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\int_{t_{n}}^{t_{n+1}} \boldsymbol{u}_{t} d t\right\|\left\|e_{u}^{n+1}\right\| \\ & \leq C \tau^{3}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}, \\ J_{5} & \leq 2 \tau\left\|\rho^{n+1}\right\|_{L^{\infty}}\left\|\nabla \boldsymbol{u}\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\frac{1}{\sigma^{n}}\right\|_{L^{\infty}}\left\|\sigma^{n} e_{u}^{n}\right\|\left\|e_{u}^{n+1}\right\| \\ & \leq C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\frac{\eta \tau}{10}\left\|\nabla e_{u}^{n+1}\right\|^{2}. \end{align} (3.46)

    Therefore, we have

    \begin{equation} \begin{split} M_{1}+M_{2}+M_{3} \leq& C \tau^{3}+C \tau^{\frac{3}{2}}\left\|\nabla e_{u}^{n}\right\|^{2}+\frac{\eta \tau}{2}\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n} e_{u}^{n}\right\|^{2}\\ &+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2. \end{split} \end{equation} (3.47)

    Using the embedding inequality and Cauchy-Schwarz inequality, we can deal with M_4 M_9

    \begin{equation} \begin{split} M_{4} = &2 \tau \lambda\left(w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-w^{n+1} \nabla \phi^{n}, e_{u}^{n+1}\right)\\ = &2 \tau \lambda\left(w\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right)+w\left(t_{n+1}\right) \nabla e_{\phi}^{n}+\nabla \phi^{n} e_{w}^{n+1}, e_{u}^{n+1}\right)\\ \leq& 2 \tau \lambda \left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\| \int_{t_{n}}^{t_{n+1}} \nabla \phi_{t}(t) d t\right\|\| e_{u}^{n+1}\|+2 \tau \lambda\left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\|\left\|e_{u}^{n+1}\right\|\\ &+2 \tau \lambda\left\|\nabla \phi^{n}\right\|_{L^2}\left\|e_w^{n+1}\right\|_{L^{3}}\left\|e_u^{n+1}\right\|_{L^{6}}\\ \leq& C \tau^{3}+\tau \eta\left\|\nabla e_{u}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}+C \tau\left\|\nabla e_{w}^{n+1}\right\|^{2}, \\ \end{split} \end{equation} (3.48)
    \begin{equation} \begin{split} M_{5}+M_{7}+M_{9}& = 2 \tau\left(R_{u}^{n+1}, e_{u}^{n+1}\right)+2\tau \lambda\left(R_{\phi}^{n+1}, e_{w}^{n+1}\right)+2\tau\left(R_{\phi}^{n+1}, e_{\phi}^{n+1}\right)\\ &\leq2\tau\left\|R_{u}^{n+1}\right\|\left\|e_{u}^{n+1}\right\|+2 \tau \lambda\left\|R_{\phi}^{n+1}\right\|\left\|e_{w}^{n+1}\right\|+2\tau\left\|R_{\phi}^{n+1}\right\|\left\|e_{\phi}^{n+1}\right\|\\ &\leq C \tau^{3}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}+\frac{\tau\lambda\gamma}{4}\left\|\nabla e_w^{n+1}\right\|^{2}+\frac{1}{4}\| e_\phi^{n+1} \|^{2}, \end{split} \end{equation} (3.49)
    \begin{equation} \begin{split} M_{6} = &-2 \tau \lambda\left( \boldsymbol{u}\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-\tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{w}^{n+1}\right)\\ = &-2 \tau \lambda\left(\left(\tilde{e}_{u}^{n+1}-e_{u}^{n+1}\right) \nabla \phi\left(t_{n+1}\right), e_w^{n+1}\right)-2 \tau\lambda\left(e_{u}^{n+1} \nabla \phi\left(t_{n+1}\right), e_{w}^{n+1}\right)\\ &-2 \tau \lambda\left(\tilde{\boldsymbol{u}}^{n+1} \int_{t_n}^{t_{n+1}} \nabla \phi_t(t) d t, e_{w}^{n+1}\right)-2 \tau \lambda \left(\tilde{\boldsymbol{u}}^{n+1} \nabla e_{\phi}^{n}, e_{w}^{n+1}\right)\\ \leq& 2 \tau \lambda\left\|\frac{1}{\sigma^{n+1}}\right\|_{L^{\infty}}\left\|\sigma^{n+1}\left(\tilde{e}_{u}^{{n+1}}-e_{u}^{n+1}\right)\right\|\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|e_{w}^{n+1}\right\|\\ &+2\tau \lambda\left\|\frac{1}{\sigma^{n+1}}\right\|_{L^{\infty}}\left\|\sigma^{n+1} e_{u}^{n+1}\right\|\left\|\nabla \phi\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|e_w^{n+1}\right\|\\ &+2 \tau \lambda\left\|\int_{t_{n}}^{t_{n+1}} \nabla \phi_{t}(t) d t\right\|_{L^{2}}\left\|\tilde{\boldsymbol{u}}^{n+1}\right\|_{L^{3}}\left\|e_{w}^{n+1}\right\|_{L^{6}}\\ &+2 \tau \lambda\left\|\nabla e_{\phi}^{n}\right\|\left\|\tilde{\boldsymbol{u}}^{n+1}\right\| _{L^{3}}\left\|e_{w}^{n+1}\right\|_{L^{6}}\\ \leq& C \tau^{3}+\frac{\tau\lambda}{2}\left\|e_w^{n+1}\right\|^{2}+C \tau\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}+\frac{\tau\lambda\gamma}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2, \\ \end{split} \end{equation} (3.50)

    and

    \begin{equation} \begin{split} M_8 & = -2 \tau\left(u\left(t_{n+1}\right) \nabla \phi(t_{n+1} )-\tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}, e_{\phi}^{n+1}\right) \\ & \leq C\tau^{3}+\frac{1}{4}\left\|e_{\phi}^{n+1}\right\|^{2}+C\tau\left\|\sigma^{n+1} e_u^{n+1}\right\|^{2}+\frac{\lambda}{2}\left\|\nabla e_{\phi}^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}.\\ \end{split} \end{equation} (3.51)

    Next, we estimate M_{10} + M_{11} and M_{12} + M_{13} as follows:

    \begin{equation} \begin{split} M_{10}+M_{11} \leq C\tau^{3}+C \tau\left(\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}\right)+\frac{\tau\gamma}{2} \left\|e_{w}^{n+1}\right\|^{2}. \end{split} \end{equation} (3.52)
    \begin{align} &\left(\tilde{G}^{n+1}, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\\ = &\tau\left(\tilde{G}^{n+1}, \gamma \triangle e_w^{n+1}-u(t_{n+1} ) \nabla \phi(t _{n+1})+\tilde{\boldsymbol{u}}^{n+1} \nabla \phi^{n}+R_{\phi}^{n+1}\right)\\ = &\tau\left(\tilde{G}^{n+1}, \gamma \triangle e_w^{n+1}-\left(\tilde{e}_{u}^{n+1}-e_{u}^{n+1}\right) \nabla \phi\left(t_{n+1}\right)-e_{u}^{n+1} \nabla \phi(t _{n+1})-\tilde{\boldsymbol{u}}^{n+1} \int_{t _n}^{t _{n+1}} \nabla \phi_t(t) d t\right)\\ &-\tau\left(\tilde{G}^{n+1}, \tilde{\boldsymbol{u}} \nabla e_{\phi}^{n}-R_{\phi}^{n+1}\right) \\ \leq &C \tau^{3}+C\tau\left\|\tilde{G}^{n+1}\right\|^{2}+\frac{\tau\gamma \lambda}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+C\tau\left\|\nabla \tilde{G}^{n+1}\right\|^{2}+\frac{1}{4}\left\|\sigma^{n+1} e_u^{n+1}\right\|^{2}+C \tau\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2\\ \leq& C\tau^{3}+\frac{\tau \gamma \lambda}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+\frac{1}{4}\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2\\ &+C\tau\left(\left\|e_{w}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{\phi}^{n}\right\|^{2}\right), \end{align} (3.53)

    and

    \begin{align} \left(\int_{t_{n}}^{t_{n+1}} \phi_{t}(t) d t, e_{\phi}^{n+1}-e_{\phi}^{n}\right)\leq C \tau^{3}+\frac{\tau \gamma \lambda}{4}\left\|\nabla e_{w}^{n+1}\right\|^{2}+\frac{1}{4}\left\|\sigma^{n+1} e^{n+1}\right\|^{2}+C \tau\|\nabla e _{\phi}^{n}\|^{2}. \end{align}

    From (3.29), combining the two inequalities above, we deduce that

    \begin{equation} \begin{split} M_{12}+M_{13} \leq&-\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|(e_{\phi}^{n+1})^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e _{\phi}^{n}\right)\right\|^{2}\right)-\frac{2 \lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e _{\phi}^{n}\right\|^{2}\\ &+C \tau^{3}+\frac{\tau\gamma \lambda}{2}\left\|\nabla e_{w}^{n+1}\right\|^{2}+\frac{1}{2}\left\|\sigma^{n+1} e_u^{n+1}\right\|^{2} \\ &+C\tau\left(\left\|e_{w}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{\phi}^{n}\right\|^{2}\right).\\ \end{split} \end{equation} (3.54)

    Plugging the above inequality into (3.41) gives

    \begin{align} &\frac{1}{2}\left(\left\|\sigma^{n+1} e_{u}^{n+1}\right\|^{2}-\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|\sigma^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right)\right\|^{2}\right)+ \eta\tau\left\|\nabla e_{u}^{n+1}\right\|^{2}+ \tau \gamma \lambda\left\|\nabla e_{w}^{n+1}\right\|^{2} \\ &+\frac{1}{2}\left(\left\|e_{\phi}^{n+1}\right\|^{2}-\left\|e_{\phi}^{n}\right\|^{2}+\left\|e_{\phi}^{n+1}-e_{\phi}^{n}\right\|^{2}\right)+ \tau \gamma\left\|e_{w}^{n+1}\right\|^{2}+\frac{\lambda}{2}\left\|\nabla e_{\phi}^{n+1}\right\|^{2}-\lambda\left\|\nabla e_{\phi}^{n}\right\|^{2}\\ &+\lambda\left\|\nabla e _{\phi}^{n+1}-\nabla e_{\phi}^{n} \right\|^{2}+\frac{\lambda}{2 \varepsilon^{2}}\left(\left\|e_{\phi}^{n+1}\right\|_{L^{4}}^{4}-\left\|e_{\phi}^{n}\right\|_{L^{4}}^{4}+\left\|(e_{\phi}^{n+1})^{2}-\left(e_{\phi}^{n}\right)^{2}\right\|^{2}\right.\\ &\left.+2\left\|e_{\phi}^{n+1}\left(e_{\phi}^{n+1}-e _{\phi}^{n}\right)\right\|^{2}\right)+\frac{2 \lambda}{\varepsilon^{2}}\left\|e_{\phi}^{n+1}-e _{\phi}^{n}\right\|^{2}\\ \leq& C \tau^{3}+C \tau^{2} \sum\limits_{m = 0}^{n}\left\|\sigma^{m} e_{u}^{m}\right\|^{2}+C\tau\left\|\sigma^{n+1}\left(e_{u}^{n+1}- \tilde{e}_{u}^{n+1} \right)\right\|^2\\ &+C\tau\left(\left\|\sigma^{n} e_{u}^{n}\right\|^{2}+\left\|e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{\phi}^{n}\right\|^{2}+\left\|\nabla e_{w}^{n+1}\right\|^{2} +\left\|e_{w}^{n+1}\right\|^{2}+\left\|e_{\phi}^{n+1}\right\|^{2}+\left\|\nabla e_{\phi}^{n+1}\right\|^{2}\right){;} \end{align} (3.55)

    for sufficiently small \tau , taking the sum of (3.55) from 0 to N -1 and using the discrete Gronwall inequality, we can obtain Lemma 6.

    Theorem 4. Suppose that the solution to (1.1) satisfies the regularity assumptions given by (3.3), and suppose that (3.1)–(3.2) are valid. For sufficiently small \tau , there are the following error estimates:

    \begin{equation} \begin{aligned} &\left\|\sigma\left(t_{n}\right)\boldsymbol{u}\left(t_{n}\right)-\sigma^{n}\boldsymbol{u}^{n}\right\|^{2}+\left\|\rho\left(t_{n}\right)-\rho^{n}\right\|^{2}+\left\|\nabla\left(w\left(t_{n}\right)-w^{n}\right)\right\|^{2} \\ &+\left\|\nabla\left(\phi\left(t_{n}\right)-\phi^{n}\right)\right\|^{2}+\eta \tau \sum\limits_{m = 1}^{n}\left\|\nabla\left(\boldsymbol{u}\left(t_{m}\right)-\boldsymbol{u}^{m}\right)\right\|^{2} \leq C \tau^{2} \end{aligned} \end{equation} (3.56)

    for all 1 \leq n \leq N .

    Proof. From Lemmas 5 and 6, we obtain

    \begin{equation} \left\|\rho\left(t_{n}\right)-\rho^{n}\right\|^{2}+\left\|\nabla\left(w\left(t_{n}\right)-w^{n}\right)\right\|^{2}+\left\|\nabla\left(\phi\left(t_{n}\right)-\phi^{n}\right)\right\|^{2} +\eta \tau \sum\limits_{m = 1}^{n}\left\|\nabla\left(\boldsymbol{u}\left(t_{m}\right)-\boldsymbol{u}^{m}\right)\right\|^{2} \leq C \tau^{2}. \end{equation} (3.57)

    Thus, we only prove that

    \begin{equation} \left\|\sigma\left(t_{n}\right)\boldsymbol{u}\left(t_{n}\right)-\sigma^{n}\boldsymbol{u}^{n}\right\|^{2} \leq C \tau^{2} . \end{equation} (3.58)

    In fact, we have

    \begin{equation} \begin{aligned} \left\|\sigma\left(t_{n}\right)\boldsymbol{u}\left(t_{n}\right)-\sigma^{n}\boldsymbol{u}^{n}\right\|^{2} & \leq C\left\|\left(\sigma\left(t_{n}\right)-\sigma^{n}\right)\boldsymbol{u}\left(t_{n}\right)\right\|^{2}+\left\|\sigma^{n} e_{u}^{n}\right\|^{2} \\ & \leq C\left(\left\|e_{\rho}^{n}\right\|^{2}+\left\|\sigma^{n} e_{u}^{n}\right\|^{2}\right) \leq C \tau^{2}. \end{aligned} \end{equation} (3.59)

    Theorem 3.6 states that both \sigma^{n}\boldsymbol{u}^{n}, \; \rho^{n}, \; \boldsymbol{u}^{n} , \phi^{n} and w^{n} are order 1 approximations to \sigma\boldsymbol{u}, \; \rho, \; \boldsymbol{u} , \phi and w in l^{\infty}\left(L^{2}(\Omega)\right) , l^{\infty}\left(L^{2}(\Omega)\right) , l^{2}\left(H_{0}^{1}(\Omega)\right) , l^{\infty}\left(H_0^{1}(\Omega)\right) and l^{\infty}\left(H_0^{1}(\Omega)\right) , respectively. Finally, we can obtain order \frac{1}{2} error estimates for p approximation in l^{\infty}\left(L^{2}(\Omega)\right) .

    Theorem 5. Under the assumptions in Theorem 4, the following holds true:

    \begin{equation} \tau\sum\limits_{m = 1}^{N}\left\|p\left(t_{m}\right)-p^{m}\right\| \leq C \tau. \end{equation} (3.60)

    Proof. Let us rewrite (3.40) as

    \begin{equation} \begin{split} -\nabla e_{p}^{n+1} = &\rho^{n} \frac{e_{u}^{n+1}-e_{u}^{n}}{\tau}+\frac{\rho^{n+1}-\rho^{n}}{\tau}\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)-\eta \triangle e_{u}^{n+1}-R_{u}^{n+1}\\ &+\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{\boldsymbol{u}}^{n+1}\\ &- \lambda w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)+ \lambda w^{n+1} \nabla \phi^{n}.\\ \end{split} \end{equation} (3.61)

    To prove the theorem, we introduce the discrete inf-sup condition, i.e.,

    \begin{equation} \beta\left\|e_{p}^{n+1}\right\| \leq \frac{\left(\nabla e_{p}^{n+1}, v\right)}{\|\nabla v\|}. \end{equation} (3.62)

    Then, we can constrain the products of the right-hand side of (3.61) with an arbitrary v \in V as follows:

    \begin{equation*} \begin{aligned} \frac{1}{\tau}\left(\rho^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right), v\right) &\leq C \tau^{-1}\left\|\sigma^{n}\left(e_{u}^{n+1}-e_{u}^{n}\right)\right\| \| \nabla v \|, \\ -\eta\left(\Delta e_{u}^{n+1}, v\right)-\left(R_{u}^{n+1}, v\right) &\leq C\left(\tau^{2}+\left\|\nabla e_{u}^{n+1}\right\|\right) \| \nabla v \|, \end{aligned} \end{equation*}

    and

    \begin{equation} \begin{split} &- \lambda\left(w\left(t_{n+1}\right) \nabla \phi\left(t_{n+1}\right)-w^{n+1} \nabla \phi^{n}, v\right)\\ = & -\lambda\left(w\left(t_{n+1}\right)\left(\nabla \phi\left(t_{n+1}\right)-\nabla \phi\left(t_{n}\right)\right)+w\left(t_{n+1}\right) \nabla e_{\phi}^{n}+\nabla \phi^{n} e_{w}^{n+1}, v\right)\\ \leq& \lambda \left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\| \int_{t_{n}}^{t_{n+1}} \nabla \phi_{t}(t) d t\right\|\| v\|+ \lambda\left\| w\left(t_{n+1}\right)\right\|_{L^{\infty}}\left\|\nabla e_{\phi}^{n}\right\|\left\|v\right\|\\ &+ \lambda\left\|\nabla \phi^{n}\right\|_{L^2}\left\|e_w^{n+1}\right\|_{L^{3}}\left\|v\right\|_{L^{6}}\\ \leq& C\left( \tau+\left\|\nabla e_{\phi}^{n}\right\|+\left\|\nabla e_{w}^{n+1}\right\|\right)\left\| \nabla v \right\|.\\ \end{split} \end{equation} (3.63)

    For the second term on the right-hand side, using \|v\|_{L^{3}} \leq\|v\|^{\frac{1}{2}}\|\nabla v\|^{\frac{1}{2}} , we have

    \begin{equation} \begin{aligned} &\left(\frac{\rho^{n+1}-\rho^{n}}{\tau}\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right), v\right) \\ &\leq C \tau^{-1}\left\|\rho^{n+1}-\rho^{n}\right\|\left\|e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right\|_{L^{3}}\|v\|_{L^{6}} \\ &\leq C \tau^{-1}\left\|\rho^{n+1}-\rho^{n}\right\|\left\|e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right\|^{\frac{1}{2}}\left\|\nabla\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)\right\|^{\frac{1}{2}}\|\nabla v\| \\ &\leq C\left\|\sigma^{n+1}\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)\right\|^{\frac{1}{2}}\left\|\nabla\left(e_{u}^{n+1}-\widetilde{e}_{u}^{n+1}\right)\right\|^{\frac{1}{2}}\|\nabla v\|. \end{aligned} \end{equation} (3.64)

    By the split method, we arrive at

    \begin{equation} \begin{aligned} &\left(\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \widetilde{\boldsymbol{u}}^{n+1}, v\right) \\ & = \left(\rho^{n+1}\left(\left(\boldsymbol{u}\left(t_{n+1}\right)-\boldsymbol{u}\left(t_{n}\right)\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), v\right)+\left(\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \tilde{e}_{u}^{n+1}, v\right) \\ &\quad+\left(\rho^{n+1}\left(e_{u}^{n} \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), v\right)+\left(e_{\rho}^{n+1}\left(\boldsymbol{u}\left(t_{n+1}\right) \nabla\right) \boldsymbol{u}\left(t_{n+1}\right), v\right). \end{aligned} \end{equation} (3.65)

    Utilizing H \ddot{o} lder's inequality and Young's inequality, we derive

    \begin{equation} \begin{split} &\left(\rho\left(t_{n+1}\right)\left(\boldsymbol{u}\left(t_{n+1}\right) \cdot \nabla\right) \boldsymbol{u}\left(t_{n+1}\right)-\rho^{n+1}\left(\boldsymbol{u}^{n} \cdot \nabla\right) \widetilde{\boldsymbol{u}}^{n+1}, v\right) \\ \leq& C\left(\tau+\left\| e_{\rho}^{n+1}\right\|+\left\|\nabla e_{u}^{n}\right\|+\left\|\nabla \tilde{e}_{u}^{n+1}\right\|\right) \| \nabla v\|.\\ \end{split} \end{equation} (3.66)

    By adding up the inequalities and incorporating (3.62), we get the desired result.

    In this section, we will present some numerical experiments to prove the validity and accuracy of our method. In the following simulation, for phase field \phi , chemical potential w and pressure p , we take the P1 finite element space (continuous piecewise linear), and for fluid velocity \boldsymbol{u} , we take the P2 finite element space. All experiments were conducted in Freefem++. We fixed \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1, T = 0.1, \rho_1 = 1 and \rho_2 = 3. The computational domain and the initial conditions were taken as

    \begin{equation*} \begin{split} &\Omega = \{(x, y)\in R^2:x^2+y^2 < 1\}, \\ &\phi_0 = cos(\pi x)cos(\pi y), \\ &\boldsymbol{u}_0 = ( \pi cos(\pi x)sin(\pi y), -\pi sin(\pi x)cos(\pi y)).\\ \end{split} \end{equation*}

    Tables 1 and 2 verify that (3.4) is the first-order convergence rate O(\tau) of (\phi, \sigma\boldsymbol{u}, \boldsymbol{u}, \rho) , which is consistent with the conclusion obtained from theoretical analysis. It is only proved in Theorem 3.2 that the pressure is the half-order convergence rate O(\tau^{\frac{1}{2}}) because of technical reasons. However, the numerical results on p in Table 1 still reach the first-order optimal convergence rate O(\tau) . Tables 3 and 4 show the convergence rate with another set of parameters.

    Table 1.  The order of temporal convergence with \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1 .
    \tau \Vert \phi\Vert_{H^1} Rate \Vert p\Vert_{L^2} Rate
    0.007812 0.102393 0.469438
    0.003906 0.0653769 0.647262 0.282333 0.733534
    0.001953 0.0367805 0.829842 0.143976 0.971574
    0.000976 0.0196276 0.906053 0.0725853 0.988076

     | Show Table
    DownLoad: CSV
    Table 2.  The order of temporal convergence with \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1 .
    \tau \Vert \sigma \boldsymbol{u}\Vert_{L^2} Rate \Vert \boldsymbol{u}\Vert_{H^1} Rate \Vert \rho\Vert_{L^2} Rate
    0.007812 0.114448 0.554732 0.0779992
    0.003906 0.0610184 0.90737 0.2956 0.90814 0.0463154 0.751968
    0.001953 0.02533 1.2684 0.132451 1.15819 0.022386 1.04889
    0.000976 0.0127883 0.98602 0.071371 0.89204 0.0107557 1.05749

     | Show Table
    DownLoad: CSV
    Table 3.  The order of temporal convergence with \eta = 0.1, \lambda = 0.2, \gamma = 0.0003, \varepsilon = 0.08 .
    \tau \Vert \phi\Vert_{H^1} Rate \Vert p\Vert_{L^2} Rate
    0.007812 0.041787 0.0794722
    0.003906 0.0215011 0.958642 0.0435045 0.869285
    0.001953 0.0106192 1.01774 0.0210772 1.04548
    0.000976 0.0054229 0.96953 0.0106356 0.986787

     | Show Table
    DownLoad: CSV
    Table 4.  The order of temporal convergence with \eta = 0.8, \lambda = 0.7, \gamma = 0.0006, \varepsilon = 0.1 .
    \tau \Vert \sigma \boldsymbol{u}\Vert_{L^2} Rate \Vert \boldsymbol{u}\Vert_{H^1} Rate \Vert \rho\Vert_{L^2} Rate
    0.007812 0.0308723 0.387768 0.0585448
    0.003906 0.0161585 0.934018 0.201806 0.942225 0.0300705 0.961192
    0.001953 0.00587585 1.45943 0.0765909 1.39772 0.00932566 1.68907
    0.000976 0.00302534 0.957701 0.0398452 0.942768 0.00468607 0.992826

     | Show Table
    DownLoad: CSV

    Figure 1 shows the evolution of the total energy at \tau = 0.02 . The downward trend of the energy curve confirms that our scheme is unconditionally energy stable. We also see a downward trend in energy when using different parameters. The energy curves for different time steps are shown in Figure 2 as a result of keeping the other parameters unchanged. It can been found that the curves are very similar, which means that the scheme is robust against different time steps.

    Figure 1.  Energy evolution with different parameters.
    Figure 2.  Energy evolution with different time step sizes.

    To solve the Cahn-Hilliard phase field model for two-phase incompressible flows with variable density, we have designed a novel time marching scheme, which can significantly improve the calculation efficiency. The method is efficient because we decoupled the pressure from the velocity and phase field. We have also proved the unconditional energy stability, presented the error analysis and provided various numerical examples to demonstrate the stability and accuracy of the scheme. In addition, the decoupling method developed in this paper is universally applicable, and this method is always applicable for the generation of an effective fully decoupling scheme.

    The authors declare they have not used Artificial Intelligence tools (AI) in the creation of this article.

    This work was supported by the Research Project Supported by Shanxi Scholarship Council of China (No. 2021-029), Shanxi Provincial International Cooperation Base and Platform Project (202104041101019), Shanxi Province Natural Science Research (202203021211129), Shanxi Province Natural Science Research (No. 202203021212249) and Special/Youth Foundation of Taiyuan University of Technology (No. 2022QN101).

    All authors declare that they have no competing interests in this paper.



    [1] M. Ablowitz, H. Segur, Solitons and inverse scattering transform, Philadelphia: SIAM, 1981. https://doi.org/10.1137/1.9781611970883
    [2] Z. Li, Z. Zhao, N-periodic wave solutions of the N = 2 supersymmetric KdV equation, Appl. Math. Lett., 160 (2025), 109313. https://doi.org/10.1016/j.aml.2024.109313 doi: 10.1016/j.aml.2024.109313
    [3] P. Xin, Z. Zhao, Y. Wang, Quasi-periodic breathers and their dynamics to the Fokas system in nonlinear optics, Wave Motion, 133 (2025), 103449. https://doi.org/10.1016/j.wavemoti.2024.103449 doi: 10.1016/j.wavemoti.2024.103449
    [4] Y. Liu, Z. Zhao, Periodic line wave, rogue waves and the interaction solutions of the (2+1)-dimensional integrable Kadomtsev-Petviashvili-based system, Chaos Soliton. Fract., 183 (2024), 114883. https://doi.org/10.1016/j.chaos.2024.114883 doi: 10.1016/j.chaos.2024.114883
    [5] M. Li, T. Xu, Dark and antidark soliton interactions in the nonlocal nonlinear Schrödinger equation with the self-induced parity-time-symmetric potential, Phys. Rev. E, 91 (2015), 033202. https://doi.org/10.1103/PhysRevE.91.033202 doi: 10.1103/PhysRevE.91.033202
    [6] M. Li, T. Xu, D. Meng, Rational solitons in the parity-time-symmetric nonlocal nonlinear Schrödinger model, J. Phys. Soc. Jpn., 85 (2016), 124001. https://doi.org/10.7566/JPSJ.85.124001 doi: 10.7566/JPSJ.85.124001
    [7] B. Guo, L. Ling, Q. Liu, High-order solutions and generalized Darboux transformations of derivative nonlinear Schrödinger equations, Stud. Appl. Math., 130 (2013), 317–344. https://doi.org/10.1111/j.1467-9590.2012.00568.x doi: 10.1111/j.1467-9590.2012.00568.x
    [8] Z. Wu, J. He, New hierarchies of derivative nonlinear Schrödinger-type equation, Rom. Rep. Phys., 68 (2016), 1425–1446. https://doi.org/10.48550/arXiv.1704.02070 doi: 10.48550/arXiv.1704.02070
    [9] Z. Zhou, Darboux transformations and global solutions for a nonlocal derivative nonlinear Schrödinger equation, Commun. Nonlinear Sci., 62 (2018), 480–488. https://doi.org/10.1016/j.cnsns.2018.01.008 doi: 10.1016/j.cnsns.2018.01.008
    [10] M. Ablowitz, Z. Musslimani, Integrable nonlocal nonlinear Schrödinger equation, Phys. Rev. Lett., 110 (2013), 064105. https://doi.org/10.1103/PhysRevLett.110.064105 doi: 10.1103/PhysRevLett.110.064105
    [11] C. Bender, S. Boettcher, Real spectra in non-Hermitian Hamiltonians having PT symmetry, Phys. Rev. Lett., 80 (1998), 5243. https://doi.org/10.1103/PhysRevLett.80.5243 doi: 10.1103/PhysRevLett.80.5243
    [12] T. Xu, H. Li, H. Zhang, M. Li, S. Lan, Darboux transformation and analytic solutions of the discrete PT-symmetric nonlocal nonlinear Schrödinger equation, Appl. Math. Lett., 63 (2017), 88–94. https://doi.org/10.1016/j.aml.2016.07.024 doi: 10.1016/j.aml.2016.07.024
    [13] D. Meng, K. Li, Darboux transformation of the second-type nonlocal derivative nonlinear Schrödinger equation, Mod. Phys. Lett. B, 33 (2019), 1950123. https://doi.org/10.1142/S0217984919501239 doi: 10.1142/S0217984919501239
    [14] A. Fokas, Integrable multidimensional versions of the nonlocal nonlinear Schrödinger equation, Nonlinearity, 29 (2016), 319. https://doi.org/10.1088/0951-7715/29/2/319 doi: 10.1088/0951-7715/29/2/319
    [15] X. Huang, L. Ling, Soliton solutions for the nonlocal nonlinear Schrödinger equation, Eur. Phys. J. Plus, 131 (2016), 148. https://doi.org/10.1140/epjp/i2016-16148-9 doi: 10.1140/epjp/i2016-16148-9
    [16] M. Ablowitz, Z. Musslimani, Inverse scattering transform for the integrable nonlocal nonlinear Schrödinger equation, Nonlinearity, 29 (2016), 915. https://doi.org/10.1088/0951-7715/29/3/915 doi: 10.1088/0951-7715/29/3/915
    [17] A. Khare, A. Saxena, New solutions of coupled nonlocal NLS and coupled nonlocal mKdV equations, Ann. Phys., 463 (2024), 169627. https://doi.org/10.1016/j.aop.2024.169627 doi: 10.1016/j.aop.2024.169627
    [18] H. Chen, H. Zhu, Higher-dimensional vector two-component solitons of a nonautonomous partially nonlocal coupled NLS model in a linear and harmonic potential, Nonlinear Dyn., 111 (2023), 581–590. https://doi.org/10.1007/s11071-022-07629-7 doi: 10.1007/s11071-022-07629-7
    [19] J. Chen, Y. Lu, B. Wang, Global Cauchy problems for the nonlocal (derivative) NLS in E^{s}_\sigma, J. Differ. Equations, 344 (2023), 767–806. https://doi.org/10.1016/j.jde.2022.11.002 doi: 10.1016/j.jde.2022.11.002
    [20] A. Regensburger, M. Miri, C. Bersch, J. Näger, G. Onishchukov, D. Christodoulides, et al., Observation of defect states in PT-symmetric optical lattices, Phys. Rev. Lett., 110 (2013), 223902. https://doi.org/10.1103/PhysRevLett.110.223902 doi: 10.1103/PhysRevLett.110.223902
    [21] B. Peng, Ş. Özdemir, F. Lei, F. Monifi, M. Gianfreda, G. Long, et al., Nonreciprocal light transmission in parity-time-symmetric whispering-gallery microcavities, Nature Phys., 10 (2014), 394–398. https://doi.org/10.1038/nphys2927 doi: 10.1038/nphys2927
    [22] L. Feng, Y. Xu, W. Fegadolli, M. Lu, J. Oliveira, V. Almeida, et al., Experimental demonstration of a unidirectional reflectionless parity-time metamaterial at optical frequencies, Nature Mater., 12 (2013), 108–113. https://doi.org/10.1038/nmat3495 doi: 10.1038/nmat3495
    [23] K. Chen, X. Deng, S. Lou, D. Zhang, Solutions of nonlocal equations reduced from the AKNS hierarchy, Stud. Appl. Math., 141 (2018), 113–141. https://doi.org/10.1111/sapm.12215 doi: 10.1111/sapm.12215
    [24] W. Liu, D. Qiu, Z. Wu, J. He, Dynamical behavior of solution in integrable nonlocal Lakshmanan-Porsezian-Daniel equation, Commun. Theor. Phys., 65 (2016), 671. https://doi.org/10.1088/0253-6102/65/6/671 doi: 10.1088/0253-6102/65/6/671
    [25] J. Rao, Y. Cheng, J. He, Rational and semirational solutions of the nonlocal Davey-Stewartson equations, Stud. Appl. Math., 139 (2017), 568–598. https://doi.org/10.1111/sapm.12178 doi: 10.1111/sapm.12178
    [26] B. Yang, J. Yang, Transformations between nonlocal and local integrable equations, Stud. Appl. Math., 140 (2018), 178–201. https://doi.org/10.1111/sapm.12195 doi: 10.1111/sapm.12195
  • This article has been cited by:

    1. Yuxia Guo, Yichen Hu, Infinitely many solutions for Hamiltonian system with critical growth, 2024, 13, 2191-950X, 10.1515/anona-2023-0134
    2. Xingyue He, Chenghua Gao, Jingjing Wang, k-convex solutions for multiparameter Dirichlet systems with k-Hessian operator and Lane-Emden type nonlinearities, 2024, 13, 2191-950X, 10.1515/anona-2023-0136
    3. Hongying Jiao, Shuhai Zhu, Jinguo Zhang, Existence of infinitely many solutions for critical sub-elliptic systems via genus theory, 2024, 16, 2836-3310, 237, 10.3934/cam.2024011
    4. Jinli Yang, Jiajing Miao, Algebraic Schouten solitons of Lorentzian Lie groups with Yano connections, 2023, 15, 2836-3310, 763, 10.3934/cam.2023037
  • Reader Comments
  • © 2025 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Metrics

Article views(633) PDF downloads(62) Cited by(0)

Other Articles By Authors

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return

Catalog