The generalized Huygens principle for the Cauchy problem of a generic non-conservative compressible two-fluid model in \(\mathbb{R}^3\) was established. This work fills a key gap in the theory, as previous results were confined to systems with full conservation laws or “equivalent" conservative structures from specific compensatory cancellations in Green’s function. Indeed, the genuinely non-conservative model studied here falls outside these categories and presents two major analytical challenges. First, its inherent non-conservative structure blocks the direct use of techniques (e.g., variable reformulation) effective for conservative systems. Second, its Green’s function contains a -1-order Riesz operator associated with the fraction densities \(\alpha^\pm\rho^\pm\), which generates a so-called Riesz wave-IV exhibiting both slower temporal decay and poorer spatial integrability compared to the standard heat kernel, necessitating novel sharp convolution estimates with the Huygens wave. To overcome these difficulties, we develop a framework for precise nonlinear coupling, including interaction of Riesz wave-IV and Huygens wave. A pivotal step is extracting enhanced decay rates for the non-conservative pressure terms \(\alpha^\pm\nabla P(\rho^\pm)\). By reformulating these terms into a product involving the fraction densities and the specific combination \(\bar\rho^-\alpha^+\rho^++\bar\rho^+\alpha^-\rho^-\), and then proving this combination decays faster than the individual densities, we meet the minimal requirements for the crucial convolution estimates. This allows us to close the nonlinear ansatz by constructing essentially new nonlinear estimates. The success of our analysis stems from the model’s special structure, particularly the equal-pressure condition. More broadly, the sharp nonlinear estimates developed herein is applicable to a wide range of non-conservative compressible fluid models.

****Key Words: Green’s function; two-fluid model; non-conservative; wave propagation.

****MSC2010: 35A09; 35B40; 35J08; 35Q35.

1 ↩︎

1.1 Background and motivation↩︎

Multi-fluid flows, involving either immiscible phases such as air-water or miscible components forming complex mixtures like emulsions, are ubiquitous in nature and industry. A classic example is interfacial wave motion between separated fluids, traditionally described by a pair of Euler equations coupled through a moving free surface: \[\label{1461} \left\{\begin{array}{l} \partial_{t} \rho_{i}+\operatorname{div}\left(\rho_{i} \mathbf{v}_{i}\right) =0, \quad i=1,2,\\ \partial_{t}\left(\rho_{i} \mathbf{v}_{i}\right)+\operatorname{div}\left(\rho_{i} \mathbf{v}_{i} \otimes \mathbf{v}_{i}\right)+\nabla p_i =-g\rho_{i} e_3\pm F_D. \end{array}\right.\tag{1}\] In this formulation, \((\rho_1, \mathbf{v}_1)\) and \((\rho_2, \mathbf{v}_2)\) denote the density and velocity fields of the upper (e.g., air) and lower (e.g., water) fluid, respectively, \(p_i\) is the pressure, \(-g\rho_i e_3\) represents the gravitational force (\(g>0\) is the gravitational acceleration and \(e_3\) the vertical unit vector), and \(F_D\) accounts for drag. The fluids are separated by a moving interface \(z=\eta(x,y,t)\), whose evolution is determined by the kinematic condition \[\partial_t\eta=u_{1,z}-u_{1,x}\partial_x \eta-u_{1, y}\partial_y \eta\label{1463}\tag{2}\] with pressure continuity across the interface.

However, under conditions such as wave breaking, the interface becomes highly fragmented, leading to a dispersed regime where a two-fluid modeling approach becomes essential. Applying a volume-averaging procedure (see [1], [2] for details) to the microscopic equations yields a macroscopic, interface-free description—the generic compressible two-fluid model. With the inclusion of viscous and capillary effects, this model can be formulated as \[\label{1464} \left\{\begin{array}{l} \partial_{t}\left(\alpha^{\pm} \rho^{\pm}\right)+\operatorname{div}\left(\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm}\right)=0, \\ \partial_{t}\left(\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm}\right)+\operatorname{div}\left(\alpha^{\pm} \rho^{\pm} \mathbf{u}^{\pm} \otimes \mathbf{u}^{\pm}\right) +\alpha^{\pm} \nabla P=\operatorname{div}\left(\alpha^{\pm} \tau^{\pm}\right)+\sigma^\pm\alpha^{\pm} \rho^{\pm}\nabla\Delta(\alpha^{\pm} \rho^{\pm}), \\ P=P^{\pm}\left(\rho^{\pm}\right)=A^{\pm}\left(\rho^{\pm}\right)^{\bar{\gamma}^{\pm}}, \end{array}\right.\tag{3}\] where \(\alpha^\pm\), \(\rho^\pm\), and \(\mathbf{u}^\pm\) denote the volume fractions, densities, and velocities of each phase, sharing a common pressure \(P\). Here, \(\alpha^++\alpha^-=1\), the viscous stress tensors are given by \[\label{1465}\tau^{\pm} = \mu^{\pm}(\nabla \mathbf{u}^{\pm} + \nabla^{t} \mathbf{u}^{\pm}) + \lambda^{\pm} (\operatorname{div} \mathbf{u}^{\pm}) \mathrm{I},\tag{4}\] with viscosity coefficients satisfying \(\mu^{\pm} > 0\) and \(2\mu^{\pm} + 3\lambda^{\pm} \geq 0\), and \(\sigma^\pm>0\) denote the capillary coefficients. For more information about this model, we refer to [3]–[14] and references therein. However, it is well–known that as far as mathematical analysis of two–fluid model is concerned, there are many technical challenges. Some of them involve, for example:

• The corresponding linear system of the model has multiple eigenvalue, which makes linear analysis of the model become quite difficult and complicated;

• Transition to single-phase regions, i.e, regions where the mass \(\alpha^{+} \rho^{+}\) or \(\alpha^{-} \rho^{-}\) becomes zero, may occur when the volume fractions \(\alpha^{\pm}\) or the densities \(\rho^{\pm}\) become zero;

• The system is non-conservative, since the non–conservative terms \(\alpha^{\pm} \nabla P^{\pm}\) are involved in the momentum equations. This brings various mathematical difficulties for us to employ methods used for single phase models to the two–fluid model.

Bresch et al. in the seminal work [15] considered the generic two-fluid model 3 with the following special density-dependent viscosities: \[\label{1467} \mu^\pm(\rho^\pm)=\mu^\pm \rho^\pm,~~~~~\lambda^\pm(\rho^\pm)=0.\tag{5}\] They obtained the global weak solutions in periodic domain with \(1<\overline{\gamma}^{\pm}< 6\). However, as indicated by themselves, their methods rely heavily on the above special density-dependent viscosities, and particularly cannot handle the case of constant viscosity coefficients as in 4 . Later, Bresch–Huang–Li [16] established the global existence of weak solutions in one space dimension without capillary effects (i.e., \(\sigma^\pm=0\)) when \(\overline{\gamma}^{\pm}>1\) by taking advantage of the one space dimension. Recently, Cui–Wang–Yao–Zhu [17] obtained the time–decay rates of classical solutions for model 3 with the following special density-dependent viscosities with equal viscosity coefficients, and equal capillary coefficients: \[\label{1468} \mu^\pm(\rho^\pm)=\nu\rho^\pm,~~~~~\lambda^\pm(\rho^\pm)=0,~~~\sigma^+=\sigma^-=\sigma.\tag{6}\] Based on the above special choice for viscosities and capillary coefficients, they can take a linear combination of model 3 to reformulate it into two \(4\times 4\) systems whose linear parts are decoupled with each other and possess the same dissipation structure as that of the compressible Navier–Stokes-Korteweg system, and then employ the similar arguments as in [18], [19] to prove their main results. However, since this reformulation played a crucial role in their analysis, the case of constant viscosities, even if the equal constant viscosities (i.e., \(\mu^\pm(\rho^\pm)=\nu,~~\lambda^\pm(\rho^\pm)=\lambda\)), cannot be handled in their settings. Very recently, The third author of this article and the coauthors in [20] investigated the optimal \(L^2\)-decay rate for all the derivatives of the solution for this model with general constant viscosities and capillary coefficients based on the spectral analysis and energy estimates. More precisely, they verified the following: Assume that \(R_{0}^{+}-1,~ R_{0}^{-}-1\in H^{\ell+1}(\mathbb{R}^3)\) and \(\mathbf{u}_{0}^{+},~ \mathbf{u}_{0}^{-}\in H^\ell(\mathbb{R}^3)\) for an integer \(\ell\geq 3\), where \(R^\pm=\alpha^\pm\rho^\pm\). Then there exists a constant \(\delta_0\) such that if \[\label{1469} K_0:=\left\|\left(R_{0}^{+}-1, R_{0}^{-}-1 \right)\right\|_{H^{\ell+1}\cap L^1}+\left\|\left(\mathbf{u}_{0}^{+}, \mathbf{u}_{0}^{-}\right)\right\|_{H^{\ell}\cap L^1} \leq \delta_0,\tag{7}\] then the Cauchy problem of 3 –4 admits a unique solution \(\left(R^{+}, u^{+}, R^{-}, u^{-}\right)\) globally in time, and the solution satisfies the following decay rate for any \(t\geq 0,\) and \(0\leq k\leq \ell\) that \[\label{14610}\left\|\nabla^k\left(\rho^{+}-\bar{\rho}^+, \rho^{-}-\bar{\rho}^-\right)(t)\right\|_{H^{\ell-k}}\lesssim K_0(1+t)^{-\frac{3}{4}-\frac{k}{2}},\tag{8}\] \[\label{14611}\left\|\nabla^{k}\left(\mathbf{u}^{+}, \mathbf{u}^{-}\right)(t)\right\|_{H^{\ell-k}} \lesssim K_0(1+t)^{-\frac{3}{4}-\frac{k}{2}},\tag{9}\] \[\label{14612}\left\|\nabla^{k}\left[\beta^+(R^+-1)+\beta^-(R^--1)\right](t)\right\|_{H^{\ell-k}} \lesssim K_0(1+t)^{-\frac{3}{4}-\frac{k}{2}},\tag{10}\] \[\label{14613}\left\|\nabla^k\left(R^{+}-1, R^{-}-1 \right)(t)\right\|_{H^{\ell-k+1}}\lesssim K_0(1+t)^{-\frac{1}{4}-\frac{k}{2}},\tag{11}\] where \(\bar{\rho}^\pm=\rho^\pm(1,1)\) denote equilibrium states of \(\rho^\pm\) respectively, and \(\beta^\pm=\sqrt{\frac{\bar{\rho}^\mp}{{\bar{\rho}}^\pm}}\).

The main purpose of this work is to establish the long-time wave propagation of classical solutions to the generic viscous-capillary two-fluid model 3 –4 . Specifically, we aim to prove that its solutions obey the so-called generalized Huygens principle—a property previously established for the one-phase compressible Navier-Stokes equations in Liu-Wang [21]. This principle states that the time-asymptotic profile of the solution decomposes into a superposition of a stationary diffusion wave (D-wave) of the form \((1+t)^{-\frac{3}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-\frac{3}{2}}\), and a moving diffusion wave (H-wave) of the form \((1+t)^{-2}\big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\big)^{-\frac{3}{2}}\), where “c" is the base sound speed. The H-wave originates from the wave-type components in the low-frequency of Green’s function. Since the H-wave decays more slowly in \(L^p(\mathbb{R}^3)\) for \(p<2\) than the D-wave, its nonlinear coupling requires”borrowing” one derivative from the nonlinear terms—a step that typically relies on the conservative (divergence) structure of the equations, as seen in earlier works [21], [22] and recent works [23], [24]. We would like to refer to other related results for the long-time behaviors of the one-phase model in [19], [25]–[33] and the kinetic equations in [34]–[39].

Compared with classical one-phase models, the analysis of the two-fluid system faces two basical structural obstacles. First, the system is inherently non-conservative due to the pressure terms \(\alpha^\pm\nabla P(\rho^\pm)\) in the momentum equations. This prevents the direct use of cancellations and reformulation techniques that are essential in conservative settings. Second, a detailed spectral analysis of the linearized system reveals that Green’s function contains a singular low-frequency term of the type \(|\xi|^{-1}e^{-|\xi|^2t}\) in the entries corresponding to the two densities. This term gives rise to a so-called Riesz wave with pointwise behavior like \((1+t)^{-1}\big(1+\frac{|x|^2}{1+t}\big)^{-1}\) (labeled Riesz wave-IV), which exhibits slower temporal decay and worse spatial regularity than the standard D-wave. Such Riesz waves also appear in other conservative systems such as the unipolar compressible Navier–Stokes–Poisson equations in Wang and Wu [40], where they verified that the electronic field containing the -1-order operator \(\nabla(-\Delta)^{-1}\) will impede the propagation of acoustic wave and ultimately results that its solution does not obey the generalized Huygens principle. On the other hand, although there exists the H-waves and Riesz wave-IV for the bipolar compressible Navier-Stokes-Poisson equations in Wu and Wang [41], these two waves do not interact with each other when dealing with the nonlinear coupling, since one can rewrite the original system in [41] into two subsystem by using some linear combinations of unknowns such that the Riesz wave-IV is only confronted with the D-wave (the propagation speed is zero).

The key innovation of this paper lies in overcoming the simultaneous interaction between the H-wave and the Riesz wave within a non-conservative framework. Unlike previous models where these waves could be decoupled, here they are intrinsically coupled through the nonlinear terms. To handle this, we develop new pointwise estimates for the nonlinear coupling of distinct wave patterns, particularly the interaction between the Riesz wave and the H-wave. A crucial technical step is the reformulation of the non-conservative pressure term \(\alpha^\pm\nabla P\) into a product structure that exposes improved decay properties, allowing us to close the nonlinear ansatz. This approach not only verifies the generalized Huygens principle for the non-conservative two-fluid system but also provides a set of convolution tools applicable to other non-conservative fluid models.

In the following, we first formulate the linearized system and present the Fourier representation of its Green’s function. We then explain how the pointwise space-time description of Green’s function, combined with refined nonlinear convolution estimates for different wave patterns, leads to the large-time asymptotics of the nonlinear problem.

1.2 New formulation of system 3 and Main Results↩︎

In this subsection, we devote ourselves to reformulating the system 3 and stating the main results. The relations between the pressures of \(\eqref{1464}_3\) implies \[\label{jpwfsmiu} {\rm d}P = s_+^2 {\rm d}\rho^+ = s_-^2 {\rm d}\rho^-, \quad {\rm where} \quad s_\pm := \sqrt{ \frac{{\rm d}P}{{\rm d}\rho^\pm}(\rho^\pm)}.\tag{12}\] Here \(s_\pm\) represent the sound speed of each phase respectively. As in [15], we introduce the fraction densities \[\label{hlkjusdc} R^\pm = \alpha^\pm \rho^\pm,\tag{13}\] which together with the relation \(\alpha^++\alpha^-=1\) leads to \[\label{ocktyrgn} {\rm d} \rho^+ = \frac{1}{\alpha_+} ({\rm d}R^+ - \rho^+ {\rm d}\alpha^+), \quad {\rm d} \rho^- = \frac{1}{\alpha_-} ({\rm d}R^- + \rho^- {\rm d}\alpha^+).\tag{14}\] From 7 –8 , we finally get \[\notag {\rm d}\alpha^+ = \frac{\alpha^-s_+^2}{\alpha^-\rho^+s_+^2+\alpha^+\rho^-s_-^2} {\rm d}R^+ - \frac{\alpha^+s_-^2}{\alpha^-\rho^+s_+^2+\alpha^+\rho^-s_-^2} {\rm d}R^-.\] Substituting the above equality into 9 yields \[\notag {\rm d}\rho^+ = \frac{s_-^2}{\alpha^-\rho^+s_+^2+\alpha^+\rho^-s_-^2} (\rho^-{\rm d}R^++\rho^+{\rm d}R^-),\] and \[\notag {\rm d}\rho^- = \frac{s_+^2}{\alpha^-\rho^+s_+^2+\alpha^+\rho^-s_-^2} (\rho^-{\rm d}R^++\rho^+{\rm d}R^-),\] which combined with \((\ref{1469})\) imply \[\label{qicoysnh} {\rm d}P = \mathcal{C}^2(\rho^- {\rm d} R^+ +\rho^+ {\rm d}R^-),\tag{15}\] where \[\notag \mathcal{C}^2 := \frac{s_+^2s_-^2}{\alpha^-\rho^+s_+^2+\alpha^+\rho^-s_-^2}, \quad {\rm and} \quad s_\pm^2 = \frac{{\rm d}P(\rho^\pm)}{{\rm d}\rho^\pm} = \tilde{\gamma}^\pm \frac{P(\rho^\pm)}{\rho^\pm}.\] Next, by using the relation: \(\alpha^+ +\alpha^- =1\) again, we can get \[\label{jbnaghpl} \frac{R^+}{\rho^+} + \frac{R^-}{\rho^-} =1, \quad {\rm and \;therefore} \quad \rho^- = \frac{R^-\rho^+}{\rho^+-R^+}.\tag{16}\] From \(\eqref{1464}_3\), we have \[\notag \varphi(\rho^+, R^+, R^-) := P(\rho^+) - P\left(\frac{R^-\rho^+}{\rho^+-R^+}\right) = 0.\] Consequently, for any given two positive constants \(\tilde{R}^+\) and \(\tilde{R}^-\), there exists \(\tilde{\rho}^+>\tilde{R}^+\) such that \[\notag \varphi(\tilde{\rho}^+, \tilde{R}^+, \tilde{R}^-)=0.\] Differentiating \(\varphi\) with respect to \(\rho^+\), we get \[\notag \frac{\partial\varphi}{\partial\rho^+}(\rho^+, R^+, R^-) = s_+^2+ s_-^2 \frac{R^-R^+}{(\rho^+-R^+)^2},\] which implies \[\notag \frac{\partial\varphi}{\partial\rho^+}(\tilde{\rho}^+, \tilde{R}^+, \tilde{R}^-)>0.\] Thus, this together with Implicit Function Theorem and 8 , 11 implies that the unknowns \(\rho^\pm\), \(\alpha^\pm\) and \(\mathcal{C}\) can be given by \[\notag \rho^\pm = \varrho^\pm(R^+,R^-), \qquad \alpha^\pm = \alpha^\pm(R^+,R^-), \quad {\rm and \;therefore} \quad \mathcal{C}=\mathcal{C}(R^+, R^-).\] We refer to [[15], pp. 614] for the details.

Note that the non-conservation of system (3 ) is mainly from two pressure terms \(\alpha^\pm\nabla P^\pm\) in the nonlinear part, since the other nonlinear terms are of the divergence form. As mentioned above, the conservation is critical when deducing the generalized Huygens’ principle. To this end, we still need to use the divergence form of the other nonlinear terms when dealing with the nonlinear coupling. Therefore, we consider the following system with respect to the unknowns \((R^\pm,\;\mathbf{m}^\pm=R^\pm u^\pm)\): \[\label{14614} \left\{\begin{array}{l} \partial_{t} R^{\pm}+\operatorname{div}\mathbf{m}^\pm=0, \\ \partial_{t}\mathbf{m}^++\operatorname{div}\big(\frac{\mathbf{m}^+\otimes\;\!\mathbf{m}^+}{R^+}\big)+\alpha^{+} \mathcal{C}^{2}\big[\rho^{-} \nabla R^{+}+\rho^{+}\nabla R^{-}\big] \\ =\operatorname{div}\big\{\alpha^{+}\big[\mu^{+}\big(\nabla \frac{\mathbf{m}^+}{R^+}+\nabla^{T} \frac{\mathbf{m}^+}{R^+}\big) +\lambda^{+} \operatorname{div} \frac{\mathbf{m}^+}{R^+} \mathbb{I}\big]\big\}+\sigma^+R^+\nabla\Delta R^+, \\ \partial_{t}\mathbf{m}^-+\operatorname{div}\big(\frac{\mathbf{m}^-\otimes\;\!\mathbf{m}^-}{R^-}\big)+\alpha^{-} \mathcal{C}^{2}\big[\rho^{-} \nabla R^{+}+\rho^{+}\nabla R^{-}\big] \\ =\operatorname{div}\big\{\alpha^{-}\big[\mu^{-}\big(\nabla \frac{\mathbf{m}^-}{R^+}+\nabla^{T} \frac{\mathbf{m}^-}{R^-}\big) +\lambda^{-} \operatorname{div} \frac{\mathbf{m}^-}{R^-} \mathbb{I}\big]\big\}+\sigma^-R^-\nabla\Delta R^-, \end{array}\right.\tag{17}\] subject to the initial condition \[\label{14615} (R^{+}, \mathbf{m}^{+}, R^-, \mathbf{m}^{-})(x, 0)=(R_{0}^{+}, \mathbf{m}_{0}^{+}, R_{0}^+, \mathbf{m}_{0}^{-})(x)\rightarrow(\bar R^{+}, \mathbf{0}, \bar R^{-}, \mathbf{0}) \quad as\quad |x|\rightarrow\infty,\tag{18}\] where two positive constants \(\bar{R}^{+}\) and \(\bar{R}^{-}\) denote the background doping profile, and in the present paper are taken as 1 for simplicity.

Now, we are in a position to state our main result.

At last, when the two pressures are the same (the common pressure) and there are not capillary forces effects, it is a more challenging open problem due to the fact that the system is partially dissipative. More recently, the first result on the global existence and the temporal decay rate of the strong (or smooth) solutions under suitably small perturbation of initial data was established in [44] in Sobolev framework and [10] in critical Besov framework, respectively. As for this case, the zero eigenvalues of Green’s matrix result in the non-dissipative properties for densities and fraction densities. Thus, it is rather hard to deduce the similar results as in the present paper due to some nonlinear terms containing the densities (or fraction densities) only exhibiting the same decay rate as the linear term, since there are no ready-made convolution estimates available on the nonlinear interaction of a Huygens wave and a non-decay diffusion wave up to now. This will be considered in future.

Now, we outline the proof strategy for Theorem 1 and highlight the main difficulties and key techniques. Note that only the two pressure terms are not of the divergence form in the nonlinear terms, to utilize the divergence form of nonlinear terms except these two terms in nonlinear analysis, we consider the unknowns \((R^\pm,\;\mathbf{m}^\pm=R^\pm \mathbf{u}^\pm)\) instead of the unknowns \((R^\pm,\;\mathbf{u}^\pm)\) in Li \(et.\;al.\) [20]. This transformation alters the structure of the nonlinear terms. To see this, by taking \(n^{\pm}=R^{\pm}-1\), one can write the corresponding linearized system of model 17 in terms of the variables \((n^+, \mathbf{m}^+, n^-, \mathbf{m}^-)\): \[\label{14627} \left\{\begin{array}{l} \partial_{t} n^++\operatorname{div}\mathbf{m}^+=0, \\ \partial_{t}\mathbf{m}^{+}+\beta_1\nabla n^++\beta_2\nabla n^--\nu^+_1\Delta \mathbf{m}^+-\nu^+_2\nabla\operatorname{div} \mathbf{m}^+-\sigma^+\nabla\Delta n^+=F_1, \\ \partial_{t} n^-+\operatorname{div}\mathbf{m}^-=0, \\ \partial_{t}\mathbf{m}^{-}+\beta_3\nabla n^++\beta_4\nabla n^--\nu^-_1\Delta \mathbf{m}^--\nu^-_2\nabla\operatorname{div} \mathbf{m}^--\sigma^-\nabla\Delta n^-=F_2, \\ \end{array}\right.\tag{19}\] where \(\nu_{1}^{\pm}=\frac{\mu^{\pm}}{\bar{\rho}^{\pm}}\), \(\nu_{2}^{\pm}=\frac{\mu^{\pm}+\lambda^{\pm}}{\bar{\rho}^{\pm}}>0\), and the coefficients \(\beta_i\) with \(i=1,2,3,4\) related to the pressure terms satisfy that \[\label{rkfxqavt} \beta_{1}=\frac{\mathcal{C}^{2}(1,1) \bar{\rho}^{-}}{\bar{\rho}^{+}},\;\; \beta_{2}=\beta_{3}=\mathcal{C}^{2}(1,1),\;\; \beta_{4}=\frac{\mathcal{C}^{2}(1,1) \bar{\rho}^{+}}{\bar{\rho}^{-}}.\tag{20}\] We shall also rewrite the nonlinear terms \(F_1\) and \(F_2\) as \[\label{14629} F_1=-\underbrace{[\alpha^+\nabla P-(\beta_1\nabla n^++\beta_2\nabla n^-)]}_{Q_1}+\cdots,\;F_2=-\underbrace{[\alpha^-\nabla P-(\beta_3\nabla n^++\beta_4\nabla n^-)]}_{Q_2}+\cdots,\tag{21}\] where “\(\cdots\)" are the terms of the divergence form in \(F_1\) and \(F_2\), and they will not bring us any new difficulty. The detailed expressions of \(F_1\) and \(F_2\) will be given in the next section.

The main difficulties stem from three interconnected features of the two-phase model 3 : (i) the presence of a nonlocal component of the form \(\frac{1}{|\xi|}e^{-|\xi|^2t}\) in the low frequency of Green function; (ii) the existence of Huygens waves in the low frequency of its Green’s function; (iii) the inherent non-conservative structure of the system. Superficially, this bears a resemblance to the bipolar Navier-Stokes-Poisson system studied by Wu and Wang [41], which also exhibits Huygens waves and a non-conservative term arising from the nonlocal electric field \(\nabla\phi=-\nabla(-\Delta)^{-1}\rho\). While the analysis in [41] relied on rewriting the system via linear combinations of unknowns to exploit a partially conservative structure, coupled with the assumption \(\nabla\phi_0\in L^1\) to enhance decay, such techniques are not directly transferable to the two-phase model 3 . The non-conservative pressure terms \(\alpha^\pm\nabla P\) in our model are intrinsically coupled and cannot be decoupled through a linear transformation, presenting a fundamental obstacle that necessitates a novel analytical framework.

The main analytical difficulties arise from the low-frequency behavior of Green’s function for the compressible part, analyzed via Hodge decomposition. The slower pointwise decay of the fraction densities in (?? ) can be attributed to the following two key reasons.

First, the low-frequency spectrum contains two distinct real eigenvalues \(\lambda_3\) and \(\lambda_4\): \[\label{14630}\begin{cases} \lambda_3=\frac{-(\beta_1\nu^-+\beta_4\nu^+)+\sqrt{(\beta_1\nu^-+\beta_4\nu^+)^2-4(\beta_1+\beta_4)(\beta_1\sigma^-+\beta_4\sigma^+)}}{2(\beta_1+\beta_4)}|\xi|^2+\mathcal{O}(|\xi|^4),\\ \lambda_4=\frac{-(\beta_1\nu^-+\beta_4\nu^+)-\sqrt{(\beta_1\nu^-+\beta_4\nu^+)^2-4(\beta_1+\beta_4)(\beta_1\sigma^-+\beta_4\sigma^+)}}{2(\beta_1+\beta_4)}|\xi|^2+\mathcal{O}(|\xi|^4),\\ \end{cases}\tag{22}\] which directly influence the diffusion structure and lead to weakened temporal decay. Their corresponding project operators \(P^3\) and \(P^4\) are given by \[\begin{align}\label{14631}P^3(\xi)=&\begin{pmatrix} \frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{r}_4}{\mathcal{R}}&-\frac{\beta_4}{\mathcal{R}|\xi|}&\frac{\beta_2\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}} &\frac{\beta_2}{\mathcal{R}|\xi|}\\ 0&-\frac{\beta_4(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_4}{\mathcal{R}}&0& \frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{r}_4}{\mathcal{R}}\\ \frac{\beta_1\beta_2(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}&\frac{\beta_2}{\mathcal{R}|\xi|}&\frac{\beta_1\beta_4(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_4}{\mathcal{R}} &-\frac{\beta_1}{\mathcal{R}|\xi|}\\ 0&\frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}&0& -\frac{\beta_1(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_4}{\mathcal{R}} \end{pmatrix}\\ &+\mathcal{O}(|\xi|),\end{align}\tag{23}\] \[\begin{align}\label{14632}P^4(\xi)=&-\begin{pmatrix} \frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}}&-\frac{\beta_4}{\mathcal{R}|\xi|}&\frac{\beta_2\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}} &\frac{\beta_2}{\mathcal{R}|\xi|}\\ 0&-\frac{\beta_4(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}}&0& \frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}\\ \frac{\beta_1\beta_2(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}&\frac{\beta_2}{\mathcal{R}|\xi|}&\frac{\beta_1\beta_4(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_3}{\mathcal{R}} &-\frac{\beta_1}{\mathcal{R}|\xi|}\\ 0&\frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}&0& -\frac{\beta_1(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_3}{\mathcal{R}} \end{pmatrix}\\ &+\mathcal{O}(|\xi|).\end{align}\tag{24}\] It is clear that there exists the term like \(\frac{1}{|\xi|}e^{-|\xi|^2t}\) in the corresponding entries of \(\hat{G}^{l}\) to the two fraction densities \(R^\pm\), which will result that \[\label{14633} \begin{array}{rl} \|n^\pm\|_{L^2}\lesssim (1+t)^{-\frac{1}{4}}. \end{array}\tag{25}\] Besides, as for the pointwise space-time behavior, the term like \(|\xi|^{-1}e^{-|\xi|^2t}\) will give rise to the so-called Riesz wave-IV, whose pointwise decay is captured by the profile \((1+t)^{-\frac{n-1}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-\frac{n-1}{2}}\) for \(x\in\mathbb{R}^n\) with \(n\geq2\). In contrast to the treatment in [41] for the bipolar compressible Navier-Stokes-Poisson system, it is unnatural here to impose additional assumptions on the initial data \(\Lambda^{-1}n^\pm_0\) in order to improve the decay rate. Instead, we aim to derive sharp pointwise estimates for the solution—a generalized Huygens principle—under the realistic situation where the \(L^2\)-decay of certain unknowns is slower than that of the classical heat kernel. The main innovation of this work lies in handling the nonlinear interaction between waves of fundamentally different types. As seen from the representations of of \(P^3\) and \(P^4\) above, and \(P^1\) and \(P^2\) below, the nonlinear convolution necessarily involves the coupling of the slowly decaying Riesz wave-IV with the faster decaying Huygens waves. This structure is essentially different from that in [41] and raises the core difficulty: How can one obtain sharp estimates for the convolution of the Riesz wave-IV with a Huygens wave under the slower decay in the ansatz of the solution to the nonlinear problem? In fact, according to the above initial propagation, we postulate the following ansatz for the nonlinear problem: \[\begin{align} \label{14634} |\partial_x^\alpha(n^+,n^-)|\lesssim (1+t)^{-\frac{2+|\alpha|}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{2+|\alpha|}{2}} +(1+t)^{-\frac{3+|\alpha|}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}. \end{align}\tag{26}\] To meet the minimal requirements for deriving the desired results from the convolution involving the Riesz wave-IV, we must extract improved decay rates for the pressure terms \(\alpha^\pm\nabla P\). Note that the other terms, except for the two pressure terms \(\alpha^\pm\nabla P\), are in divergence form (conservative structure). Hence, we focus primarily on the nonlinear estimates for the pressure terms \(\alpha^\pm\nabla P\). At first glance, \(\alpha^\pm\nabla P \sim n^-\nabla n^+ + n^+\nabla n^- +n^+\nabla n^++ n^-\nabla n^-\). Then, according to the above ansatz and the fact that the entry \(G_{12}\) in Green’s function corresponding to \(n^+\) contains the Riesz wave-IV (which has the worst decay in both time and space), one encounters the following nonlinear convolution of the Riesz wave-IV and Huygens waves: \[\begin{align} \mathcal{N}_1= \int_0^t\!\int_{\mathbb{R}^3}(1+t-\tau)^{-1}\Big(1+\frac{|x-y|^2}{1+t-\tau}\Big)^{-1}(1+\tau)^{-\frac{7}{2}}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau.\label{9467} \end{align}\tag{27}\] When analyzing the pointwise behavior, a central task is to study the space-time convolution between the slowly decaying Riesz wave–IV and the nonlinear Huygens waves. A standard approach involves partitioning the spatial domain into the near-field \((x,t)\in D_1=\{|x|\leq \sqrt{1+t}\}\) and the wave region \(D_2=\{||x|-{\rm c}t|\leq \sqrt{1+t}\}\). However, even within the simpler region \(D_1\), the time integration must be split into the sub-intervals: \(0\leq \tau\leq\frac{t}{2}\) and \(\frac{t}{2}\leq\tau\leq t\), denoted by \(\mathcal{N}_{11}\) and \(\mathcal{N}_{12}\), respectively. A direct estimate for the latter yields \[\begin{align} \tag{28} \mathcal{N}_{12}&=&\int_0^\frac{t}{2}\!\int_{\mathbb{R}^3}(1+t-\tau)^{-1}\Big(1+\frac{|x-y|^2}{1+t-\tau}\Big)^{-1}(1+\tau)^{-\frac{7}{2}}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau\nonumber\\ &\lesssim& (1+t)^{-1}\int_0^{\frac{t}{2}}(1+\tau)^{-\frac{7}{2}}(1+\tau)^{\frac{5}{2}}d\tau\nonumber\\ &\lesssim& (1+t)^{-1}\ln(1+t)\lesssim (1+t)^{-1}\ln(1+t)\bigg(1+\frac{|x|^2}{1+t}\bigg)^{-N},\tag{29} \end{align}\] which introduces a logarithmic growth factor. This result is strictly weaker than the decay rate prescribed by the ansatz (81 ) for the unknowns \(n^\pm\), presenting a major obstacle to closing the nonlinear argument.

To overcome this difficulty, we should conduct a refined analysis of the nonlinear terms \(Q_1=(\beta_1\nabla n^++\beta_2\nabla n^-)-\alpha^+\nabla P\) and \(Q_2=(\beta_3\nabla n^++\beta_4\nabla n^-)-\alpha^-\nabla P\) in \(F_1\) and \(F_2\) arising from two pressure terms. Through careful observation and algebraic reformulation, we find that they can be recast as: \[Q_1\sim n^+( \bar\rho^-\nabla n^++\bar\rho^+\nabla n^-)+\cdots,\;\;\;Q_2\sim n^-(\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-)+\cdots,\] where “\(\cdots\)" denote cubic terms with faster decay. The crucial insight lies in the structure of the linear propagators \(P^3\) and \(P^4\). A hidden cancellation between their first and third rows endows the combined quantity \(\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-\) with an extra \((1+t)^{-\frac{1}{2}}\)-decay factor: based on the cancellation of the first row and third row in \(P^3\) and \(P^4\), that is \[\begin{align} \label{14636} |\bar\rho^-n^++\bar\rho^+n^-| \lesssim\bigg\{(1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg\}. \end{align}\tag{30}\] This cancellation is fundamental and distinguishes our analysis from prior work. It transforms the problematic convolution into a new, refined form with significantly improved integrability: \[\begin{align} \mathcal{N}_1= \int_0^t\!\int_{\mathbb{R}^3}(1+t-\tau)^{-1}\Big(1+\frac{|x-y|^2}{1+t-\tau}\Big)^{-1}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau.\label{94610} \end{align}\tag{31}\] Applying the precise convolution estimate \(K_4\) from Lemma 11 to \(\mathcal{N}_1\) finally yields the sharp, desired bound: \[\begin{align} \mathcal{N}_1\lesssim (1+t)^{-1}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-\frac{3}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1},\label{94611} \end{align}\tag{32}\] which is fully consistent with the original ansatz and devoid of any logarithmic loss.

The second difficulty comes from the non-conservative structure of the system (3 ) and the wave operators in the low frequency of Green’s function for the compressible part \(\hat{G}^{l}\). In particular, we have \[\label{gnqldkmc} \begin{cases} \lambda_1=-\frac{\beta_1\nu^++\beta_4\nu^-}{2(\beta_1+\beta_4)}|\xi|^2+\mathrm{i}\sqrt{\beta_1+\beta_4}|\xi|+\mathcal{O}(|\xi|^3),\\ \lambda_2=-\frac{\beta_1\nu^++\beta_4\nu^-}{2(\beta_1+\beta_4)}|\xi|^2-\mathrm{i}\sqrt{\beta_1+\beta_4}|\xi|+\mathcal{O}(|\xi|^3),\\ \end{cases}\tag{33}\] and their corresponding project operators \(P^1\) and \(P^2\) are given by \[\label{hwyfzask}P^1(\xi)=\begin{pmatrix} \frac{\beta_1}{2(\beta_1+\beta_4)}&\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)} &\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ -\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_1}{2(\beta_1+\beta_4)}&-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_2}{2(\beta_1+\beta_4)}\\ \frac{\beta_2}{2(\beta_1+\beta_4)}&\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_4}{2(\beta_1+\beta_4)} &\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ -\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)}&-\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_4}{2(\beta_1+\beta_4)} \end{pmatrix}+\mathcal{O}(|\xi|),\tag{34}\] \[\label{xprukzjt}P^2(\xi)=\begin{pmatrix} \frac{\beta_1}{2(\beta_1+\beta_4)}&-\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)} &-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ \frac{\beta_1}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_1}{2(\beta_1+\beta_4)}&\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_2}{2(\beta_1+\beta_4)}\\ \frac{\beta_2}{2(\beta_1+\beta_4)}&-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_4}{2(\beta_1+\beta_4)} &-\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ \frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)}&\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_4}{2(\beta_1+\beta_4)} \end{pmatrix}+\mathcal{O}(|\xi|).\tag{35}\] We see that the leading terms in each entry of \(P^1\) and \(P^2\) are constants, which is the same as the compressible Navier-Stokes equations. However, the lack of a divergence form in the nonlinear terms prevents us from directly applying classical convolution estimates for Huygens waves in [21], [22]. To see this, let us recall the classical nonlinear convolution estimates on the interaction of Huygens waves in [21], [22]: \[\begin{align} &&\int_{0}^{t}\int_{\mathbb{R}^3}\uwave{(1+t-\tau)^{-\frac{5}{2}}}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1\!+\!t\!-\!\tau}\Big)^{-N}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1\!+\!\tau}\Big)^{-3}dyd\tau\nonumber\\ &\lesssim& (1+t)^{-2}\Big(\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}+\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-\frac{3}{2}}\Big),\label{14637} \end{align}\tag{36}\] where the extra \((1+t-\tau)^{-\frac{1}{2}}\)-decay in the underline term arises precisely from the conservative (divergence) structure of the nonlinearity. Since the momenta \(\mathbf{m}^\pm\) decay at the heat‑kernel rate in \(L^2\) [20], refining estimate (36 ) in the usual way seems unattainable. To resolve this difficulty, we devise a bold strategy that preserves the temporal decay by strategically sacrificing spatial decay in sub‑region estimates, especially near the wave cone. This compensates for the derivative loss induced by the non‑conservative terms. Through a meticulous and lengthy analysis of multiple space‑time subregions, we ultimately establish the following key convolution estimate: \[\begin{align} &\int_{0}^{t}\!\int_{\mathbb{R}^3}(1+t-\tau)^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau\nonumber\\[3mm] \lesssim &\;(1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}+(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}.\label{14638} \end{align}\tag{37}\] The proof, detailed in Lemma 11, relies on a delicate partition of the integration domain and optimized pointwise bounds, enabling us to recover the sharp decay needed for the ansatz despite the system’s non‑conservation.

Discussion

From the linear analysis, it is expected that improved decay rates for the two fractional momenta \(\mathbf{m}^\pm\) in (?? ) can be achieved if two distinct linear combinations, \(\mathbf{m}^+-a_1\mathbf{m}^-\) and \(\mathbf{m}^+-a_2\mathbf{m}^-\), can be constructed to compensate for the non-conservative structure arising from the two pressure terms in the momentum equations. Ultimately, this would yield the refined pointwise estimates (?? ) for \(\mathbf{m}^\pm\): \[\label{14639} |\mathbf{m}^\pm|\lesssim (1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}+(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-\frac{3}{2}}.\tag{38}\]

The first combination is found through a direct observation of the low-frequency projection matrices \(P^1\) and \(P^2\). A crucial structural property, \(\beta_1\beta_4=\beta_2^2\), implies that the vectors \((P_{22}^1,P_{24}^1)\) and \((P_{42}^1,P_{44}^1)\) associated with \(\mathbf{m}^\pm\) are proportional. This allows us to immediately choose \[\mathbf{m}^+ - \frac{\beta_1}{\beta_2}\mathbf{m}^-\] as the first combination. In this combination, the leading low-frequency terms become \(|\xi|\), aligning it with a divergence structure and leading to the desired decay: \[\label{14640} \Big|\mathbf{m}^+-\frac{\beta_1}{\beta_2}\mathbf{m}^-\Big|\lesssim (1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}+(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-\frac{3}{2}}.\tag{39}\]

Finding a second, genuinely different linear combination is more subtle, as the proportionality used above cannot be applied again. Our approach leverages another good structure of the original system: the conservation of the total fractional momentum \(\mathbf{m}^++\mathbf{m}^-\) of the two-fluid model with equal pressures in (1 ), combined with detailed Green’s function estimates. We denote the nonlinear couplings containing the pressure terms for \(\mathbf{m}^\pm\) by \(\tilde{\mathbf{m}}^\pm\), expressed via the low-frequency Green’s functions \(G_{ij}^\ell\): \[\label{14641} \begin{align} \tilde{\mathbf{m}}^+&=\int_0^t G_{22}^\ell(\cdot,t-\tau)\ast_x Q_1(\cdot,\tau)d\tau+\int_0^t G_{24}^\ell(\cdot,t-\tau)\ast_x Q_2(\cdot,\tau)d\tau+\cdots,\\ \tilde{\mathbf{m}}^-&=\int_0^t G_{42}^\ell(\cdot,t-\tau)\ast_x Q_1(\cdot,\tau)d\tau+\int_0^t G_{44}^\ell(\cdot,t-\tau)\ast_x Q_2(\cdot,\tau)d\tau+\cdots, \end{align}\tag{40}\] where \(Q_1\) and \(Q_2\) are non-divergence terms originating from the pressures (43 ), and the rest terms don’t affect the results. We then construct a candidate combination \(\mathbf{m}^+-a_2\mathbf{m}^-\): \[\label{14642} \begin{align} \tilde{\mathbf{m}}^+-a_2\tilde{\mathbf{m}}^- = &\int_0^t (G_{22}^l-a_2G_{42}^l)(\cdot,t-\tau)\ast \underbrace{(Q_1+Q_2)}_{\text{divergence form}}d\tau \\ &+ \int_0^t \underbrace{(-G_{22}^l+G_{24}^l+a_2G_{42}^l-a_2G_{44}^l)}_{=:G^*}\ast_x Q_2 \, d\tau +\cdots, \end{align}\tag{41}\] where “\(\cdots\)” denotes middle and high-frequency parts not affecting the main Huygens wave analysis. Since \(Q_1+Q_2\) is of divergence form, the first integral is manageable. The critical step is to choose \(a_2\) such that the kernel \(G^*\) has a leading term of order \(|\xi|\) (not a constant). From the representation of \(P^1\) or \(P^2\), this requires \[\label{14643} -\beta_1+\beta_2+a_2\beta_2-a_2\beta_4=0,\tag{42}\] which yields \(a_2=\frac{\beta_1-\beta_2}{\beta_2-\beta_4}\). Using the fundamental relation \(\beta_1\beta_4=\beta_2^2\), we find \[\frac{\beta_1-\beta_2}{\beta_2-\beta_4} \color{red}{=} \frac{\beta_1}{\beta_2},\] provided \(\beta_1 \neq \beta_2\). This shows that the second candidate combination degenerates to the first one. Therefore, this direct algebraic approach fails to produce a distinct second linear combination.

Conclusion and Outlook: The discovery of a genuinely different second linear combination for \(\mathbf{m}^\pm\) remains an open and interesting problem. Its construction likely requires a deeper exploitation of the system’s structure beyond the simple algebraic cancellation attempted here. Successfully finding it is essential to finally establish the refined pointwise estimates (38 ) for the individual momenta \(\mathbf{m}^\pm\), a task left for future work.

2 ↩︎

2.1 Linearization and Reformulation↩︎

In this subsection, we first reformulate the system. Setting \[n^{\pm}=R^{\pm}-1,\] the Cauchy problem 17 –18 can be reformulated as \[\label{2461} \left\{\begin{array}{l} \partial_{t} n^++\operatorname{div}\mathbf{m}^+=0, \\ \partial_{t}\mathbf{m}^{+}+\beta_1\nabla n^++\beta_2\nabla n^--\nu^+_1\Delta \mathbf{m}^+-\nu^+_2\nabla\operatorname{div} \mathbf{m}^+-\sigma^+\nabla\Delta n^+=F_1, \\ \partial_{t} n^-+\operatorname{div}\mathbf{m}^-=0, \\ \partial_{t}\mathbf{m}^{-}+\beta_3\nabla n^++\beta_4\nabla n^--\nu^-_1\Delta \mathbf{m}^--\nu^-_2\nabla\operatorname{div} \mathbf{m}^--\sigma^-\nabla\Delta n^-=F_2, \\ \end{array}\right.\tag{43}\] where \(\nu_{1}^{\pm}=\frac{\mu^{\pm}}{\bar{\rho}^{\pm}}\), \(\nu_{2}^{\pm}=\frac{\mu^{\pm}+\lambda^{\pm}}{\bar{\rho}^{\pm}}>0\), \(\beta_{1}=\frac{\mathcal{C}^{2}(1,1) \bar{\rho}^{-}}{\bar{\rho}^{+}}\), \(\beta_{2}=\beta_{3}=\mathcal{C}^{2}(1,1)\), \(\beta_{4}=\frac{\mathcal{C}^{2}(1,1) \bar{\rho}^{+}}{\bar{\rho}^{-}}\), which imply \[\label{246140141} \beta_{1}\beta_{4}=\beta_{2}\beta_{3}=\beta_{2}^2=\beta_{3}^2.\tag{44}\] Additionally, the nonlinear terms are given by \[\label{2462} \begin{array}{rl} F_{1}=&-\operatorname{div}\big(\frac{\mathbf{m}^+\otimes \mathbf{m}^+}{n^++1}\big)+\sigma^+R^+\nabla\Delta R^+\\ &+\mu^+{\rm div}\Big\{\frac{n^+}{\rho^+}\nabla\frac{\mathbf{m}^+}{n^++1}-\frac{1}{\rho^+}\nabla\frac{n^+\mathbf{m}^+}{n^++1}+\big(\frac{1}{\rho^+}-\frac{1}{\bar{\rho}^+}\big)\nabla \mathbf{m}^+\Big\} \\ &+\mu^+{\rm div}\Big\{\frac{n^+}{\rho^+}\nabla^t\frac{\mathbf{m}^+}{n^++1}-\frac{1}{\rho^+}\nabla^t\frac{n^+\mathbf{m}^+}{n^++1}+\big(\frac{1}{\rho^+}-\frac{1}{\bar{\rho}^+}\big)\nabla^t \mathbf{m}^+\Big\} \\ &+\lambda^+{\rm div}\Big\{\frac{n^+}{\rho^+}{\rm div}\frac{\mathbf{m}^+}{n^++1}I_3-\frac{1}{\rho^+}{\rm div}\frac{n^+\mathbf{m}^+}{n^++1}I_3+\Big(\frac{1}{\rho^+}-\frac{1}{\bar{\rho}^+}\Big){\rm div} \mathbf{m}^+I_3\Big\}\\ &-\underbrace{\big\{\mathcal{C}^{2}n^+\Big(\frac{\rho^-}{\rho^+}\nabla n^++\nabla n^-\Big)+\Big(\frac{\mathcal{C}^2\rho^-}{\rho^+}-\frac{\mathcal{C}^2\rho^-}{\rho^+}(1,1)\Big)\nabla n^++(\mathcal{C}^2-\mathcal{C}^2(1,1))\nabla n^-\big\}}_{\alpha^+\nabla P-(\beta_1\nabla n^++\beta_2\nabla n^-):=Q_1}, \end{array}\tag{45}\] and \[\label{2463} \begin{array}{rl} F_{2}=&-\operatorname{div}\big(\frac{\mathbf{m}^-\otimes \mathbf{m}^-}{n^-+1}\big)+\sigma^-R^-\nabla\Delta R^-\\ &+\mu^-{\rm div}\Big\{\frac{n^-}{\rho^-}\nabla\frac{\mathbf{m}^-}{n^-+1}-\frac{1}{\rho^-}\nabla\frac{n^-\mathbf{m}^-}{n^-+1}+\big(\frac{1}{\rho^-}-\frac{1}{\bar{\rho}^-}\big)\nabla \mathbf{m}^-\Big\} \\ &+\mu^-{\rm div}\Big\{\frac{n^-}{\rho^-}\nabla^t\frac{\mathbf{m}^-}{n^-+1}-\frac{1}{\rho^-}\nabla^t\frac{n^-\mathbf{m}^-}{n^-+1}+\big(\frac{1}{\rho^-}-\frac{1}{\bar{\rho}^-}\big)\nabla^t \mathbf{m}^-\Big\} \\ &+\lambda^-{\rm div}\Big\{\frac{n^-}{\rho^-}{\rm div}\frac{\mathbf{m}^-}{n^-+1}I_3-\frac{1}{\rho^-}{\rm div}\frac{n^-\mathbf{m}^-}{n^-+1}I_3+\big(\frac{1}{\rho^-}-\frac{1}{\bar{\rho}^-}\big){\rm div} \mathbf{m}^-I_3\Big\}\\ &-\underbrace{\big\{\mathcal{C}^{2}n^-\big(\nabla n^++\frac{\rho^-}{\rho^+}\nabla n^-\big)+\big(\frac{\mathcal{C}^2\rho^+}{\rho^-}-\frac{\mathcal{C}^2\rho^+}{\rho^-}(1,1)\big)\nabla n^-+(\mathcal{C}^2-\mathcal{C}^2(1,1))\nabla n^+\big\}}_{\alpha^-\nabla P-(\beta_3\nabla n^++\beta_4\nabla n^-):=Q_2}, \end{array}\tag{46}\] where the last terms \(Q_1\) and \(Q_2\) in \(F_1\) and \(F_2\) have been given in (21 ) in another simplified form.

In terms of the semigroup theory for evolutionary equation, we will investigate the following initial value problem for the corresponding linear system of 43 on the unknown \(U:=(n^+,\mathbf{m}^+,n^-,\mathbf{m}^-)\): \[\begin{cases} U_t=\mathcal{B}U,\\ U\big|_{t=0}={U}_0, \end{cases} \label{2464}\tag{47}\] where the operator \(\mathcal{B}\) is given by \[\nonumber\mathcal{B}=\begin{pmatrix} 0&-\text{div}&0&0\\ -\beta_1\nabla+\sigma^+ \nabla\Delta&\nu_{1}^{+} \Delta+\nu_{2}^{+} \nabla \otimes \nabla&-\beta_2\nabla&0\\ 0&0&0&-\text{div}\\ -\beta_3\nabla&0&-\beta_4\nabla+\sigma^- \nabla\Delta&\nu_{1}^{-} \Delta+\nu_{2}^{-} \nabla \otimes \nabla \end{pmatrix}.\] Applying Fourier transform to the system 47 , one has \[\label{2465} \begin{cases} \widehat {{U}}_t=\mathcal{B}(\xi)\widehat {U},\\ \widehat {U}\big|_{t=0}=\widehat U_0=(\widehat {n^+_0}, \widehat{ \mathbf{m}^+_0}, \widehat {n^-_0}, \widehat {\mathbf{m}^-_0}), \end{cases}\tag{48}\] where \(\widehat {U}(\xi,t)=\mathcal{F}(U(x,t))\), \(\xi=(\xi_1,\xi_2,\xi_3)^t\) and \(\mathcal{B}(\xi)\) is defined by \[\label{2466} \mathcal{B}(\xi)=\begin{pmatrix} 0&-{\rm i}\xi^t&0&0\\ -\mathrm{i}\beta_1\xi-{\rm i}\sigma^+|\xi|^2\xi&-\nu_1^+|\xi|^2\mathbb{I}_{3\times 3}-\nu_2^+\xi\otimes\xi&-{\rm i}\beta_2\xi&0\\ 0&0&0&-{\rm i}\xi^t\\ - {\rm i}\beta_3\xi&0&- {\rm i}\beta_4-{\rm i}\sigma^-|\xi|^2\xi&-\nu_1^-|\xi|^2\mathbb{I}_{3\times 3}-\nu_2^-\xi\otimes\xi \end{pmatrix}.\tag{49}\] To facilitate narrative in later use, we also use the definition of Green’s function \(G(x,t)\) with the following standard form as our previous works: \[\label{246640141} \left\{\begin{array}{l} G_{t}=\mathcal{A} G, \\ G|_{t=0}=\delta_0(x)\mathbb{I}_{8\times8}. \end{array}\right.\tag{50}\]

To derive the linear estimates, by using a real method as in [45], one need to make a detailed analysis on the properties of Green’s function. Since the system 57 has eight equations and the matrix \(\mathcal{B}(\xi)\) may not be diagonalizable. In order to facilitate the analysis, we take Hodge decomposition to system 47 such that it can be decoupled into three systems. One has four equations, and the other two are classic heat equations.

To begin with, let \(\varphi^{\pm}=\Lambda^{-1}{\rm div}\mathbf{m}^{\pm}\) be the “compressible part" of the momenta \(\mathbf{m}^{\pm}\), and denote \(\phi^{\pm}=\Lambda^{-1}{\rm curl}\mathbf{m}^{\pm}\) (with \(({\rm curl} z)_i^j =\partial_{x_j}z^i-\partial_{x_i}z^j\)) by the”incompressible part" of the momenta \(\mathbf{m}^{\pm}\). Then, we can divide the system 47 into the compressible part \[\label{2467} \begin{cases} \partial_t{n^+}+\Lambda{\varphi^+}=0,\\ \partial_t{\varphi^+}-\beta_1\Lambda{n^+}-\beta_2\Lambda{n^-}+\nu^+\Lambda^2{\varphi^+}-\sigma^+\Lambda^3{n^+}=0,\\ \partial_t{n^-}+\Lambda{\varphi^-}=0,\\ \partial_t{\varphi^-}-\beta_3\Lambda{n^+}-\beta_4\Lambda{n^-}+\nu^-\Lambda^2{\varphi^-}-\sigma^-\Lambda^3{n^-}=0,\\ (n^+, \varphi^+, n^-, \varphi^-)\big|_{t=0}=({n}^+_0, \Lambda^{-1}{\rm div}{\mathbf{m}}^{+}_0, {n}^-_0, \Lambda^{-1}{\rm div}{\mathbf{m}}^{-}_0)(x),\\ \end{cases}\tag{51}\] with its Green’s function \(\tilde{G}(x,t)\), and the incompressible part \[\label{2468} \begin{cases} \partial_t\phi^++\nu^+_1\Lambda^2\phi^+=0,\\ \partial_t\phi^-+\nu^-_1\Lambda^2\phi^-=0,\\ (\phi^+,\phi^-)\big|_{t=0}=(\Lambda^{-1}{\rm curl}{\mathbf{m}}^{+}_0, \Lambda^{-1}{\rm curl}{\mathbf{m}}^{-}_0)(x) \end{cases}\tag{52}\] where \(\nu^{\pm}=\nu^{\pm}_1+\nu^{\pm}_2\). It is obvious that the incompressible part are two heat kernels: \[\begin{align} \label{246840141} \hat{\phi}^+=e^{-\nu_1^+|\xi|^2t},\;\hat{\phi}^-=e^{-\nu_1^-|\xi|^2t}, \end{align}\tag{53}\] and hence, we mainly focus on the compressible part of Green’s function.

2.2 Compressible part↩︎

In view of the semigroup theory, we may represent the IVP 47 for \(\mathcal{U}=(n^+, \varphi^+, n^-, \varphi^-)^t\) as \[\label{2469} \begin{cases} \mathcal{U}_t=\mathcal{A}\mathcal{U},\\ \mathcal{U}\big|_{t=0}=\mathcal{U}_0, \end{cases}\tag{54}\] where the operator \(\mathcal{A}\) is defined by \[\nonumber\mathcal{A}=\begin{pmatrix} 0&-\Lambda&0&0\\ \beta_1\Lambda+\sigma^+\Lambda^3&-\nu^+\Lambda^2&\beta_2\Lambda&0\\ 0&0&0&-\Lambda\\ \beta_3\Lambda&0&\beta_4\Lambda+\sigma^-\Lambda^3&-\nu^-\Lambda^2 \end{pmatrix}.\] Taking Fourier transform to system 54 , we obtain \[\label{24610} \begin{cases} \widehat {\mathcal{U}}_t=\mathcal{A}(\xi)\widehat {\mathcal{U}},\\ \widehat {\mathcal{U}}\;\big|_{t=0}=\widehat {\mathcal{U}}_0, \end{cases}\tag{55}\] where \(\widehat {\mathcal{U}}(\xi,t)=\mathfrak{F}({\mathcal{U}}(x,t))\) and \(\mathcal{A}(\xi)\) is given by \[\label{24611} \mathcal{A}(\xi)=\begin{pmatrix} 0&-|\xi|&0&0\\ \beta_1|\xi|+\sigma^+|\xi|^3&-\nu^+|\xi|^2&\beta_2|\xi|&0\\ 0&0&0&-|\xi|\\ \beta_3|\xi|&0&\beta_4|\xi|+\sigma^-|\xi|^3&-\nu^-|\xi|^2 \end{pmatrix}.\tag{56}\] Its eigenvalues satisfy \[\label{24612} \begin{align}&{\rm det}(\lambda{\rm I}-\mathcal{A}(\xi))\\ &=\lambda^4+(\nu^+|\xi|^2+\nu^-|\xi|^2)\lambda^3+(\beta_1|\xi|^2+\beta_4|\xi|^2+\sigma^+|\xi|^4+\sigma^-|\xi|^4 +\nu^+\nu^-|\xi|^4)\lambda^2\\&\quad+(\beta_1\nu^-|\xi|^4+\beta_4\nu^+|\xi|^4+\nu^+\sigma^-|\xi|^6+\nu^-\sigma^+|\xi|^6)\lambda+\beta_1\sigma^-|\xi|^6+\beta_4\sigma^+|\xi|^6+\sigma^+\sigma^-|\xi|^8\\ &=0. \end{align}\tag{57}\]

Low frequency part.

From discriminant of the quartic equation in one unknown, we know that when \(|\xi|\ll1\) there exist two real roots and a pair of conjugate imaginary roots. Then, by using the method of undetermined coefficients, we have the following for the spectral in the low frequency part:

We divide the issue into two cases. Since it is basically similar to [20], we just give the sketch.

Case 1: \(\lambda_3\neq \lambda_4\). In this case, we know the four eigenvalues are different, then one can immediately get the semigroup \(\mathrm{e}^{t \mathcal{A}}\) (denoted by Green’s function \(\tilde{G}(x,t)\)) can be decomposed into \[\label{nwbogcux} \tilde{G}(x,t):=\mathrm{e}^{t \mathcal{A}(\xi)}=\sum_{i=1}^{4} \mathrm{e}^{\lambda_{i} t} P^{i}(\xi),\;\;|\xi|\ll1,\tag{58}\] and the projector \(P^{i}(\xi)=\prod_{j \neq i} \frac{\mathcal{A}(\xi)-\lambda_{j} \mathbb{I}}{\lambda_{i}-\lambda_{j}}, \quad i, j=1,2,3,4\).

A direct and tedious computation gives that \[\begin{align} \label{24614} P^1(\xi)=\bar{P}^2(\xi)=\begin{pmatrix} \frac{\beta_1}{2(\beta_1+\beta_4)}&\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)} &\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ -\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_1}{2(\beta_1+\beta_4)}&-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_2}{2(\beta_1+\beta_4)}\\ \frac{\beta_2}{2(\beta_1+\beta_4)}&\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_4}{2(\beta_1+\beta_4)} &\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ -\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)}&-\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_4}{2(\beta_1+\beta_4)} \end{pmatrix}+\mathcal{O}(|\xi|), \end{align}\tag{59}\] \[\begin{align} \label{24615} P^3(\xi)=&\begin{pmatrix}\! \frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_4}{\mathcal{R}}&-\frac{\beta_4}{\mathcal{R}|\xi|}&\frac{\beta_2\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}} &\frac{\beta_2}{\mathcal{R}|\xi|}\\ 0&-\frac{\beta_4(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_4}{\mathcal{R}}&0& \frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}\\ \frac{\beta_1\beta_2(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}&\frac{\beta_2}{\mathcal{R}|\xi|}&\frac{\beta_1\beta_4(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_4}{\mathcal{R}} &-\frac{\beta_1}{\mathcal{R}|\xi|}\\ 0&\frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}&0& -\frac{\beta_1(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_4}{\mathcal{R}} \!\end{pmatrix}\nonumber\\ &+\mathcal{O}(|\xi|), \end{align}\tag{60}\] \[\begin{align} \label{24616} P^4(\xi)=&\!-\!\!\begin{pmatrix}\!\! \frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}}\! &\! -\frac{\beta_4}{\mathcal{R}|\xi|} & \frac{\beta_2\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}} & \frac{\beta_2}{\mathcal{R}|\xi|}\\ 0 \!&\! -\frac{\beta_4(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}} & 0 & \frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}\\ \frac{\beta_1\beta_2(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}} \!&\! \frac{\beta_2}{\mathcal{R}|\xi|}&\frac{\beta_1\beta_4(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_3}{\mathcal{R}} &-\frac{\beta_1}{\mathcal{R}|\xi|}\\ 0 \!&\! \frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}&0& -\frac{\beta_1(\beta_4\nu^++\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_3}{\mathcal{R}} \!\!\end{pmatrix}\nonumber\\ &+\mathcal{O}(|\xi|). \end{align}\tag{61}\] Here \(\mathcal{R}\!=\!\sqrt{(\beta_1\nu^-\!+\!\beta_4\nu^+)^2\!-\!4(\beta_1\!+\!\beta_4)(\beta_1\sigma^-\!+\!\beta_4\sigma^+)}\), \(\tilde{\lambda}_4=\frac{-(\beta_1\nu^-+\beta_4\nu^+)+\mathcal{R}}{2(\beta_1+\beta_4)}\) and \(\tilde{\lambda}_4=\frac{-(\beta_1\nu^-+\beta_4\nu^+)-\mathcal{R}}{2(\beta_1+\beta_4)}\).

Case 2: \(\lambda_3=\lambda_4\). In this case, we have the following decomposition: \[\begin{align}\label{24630}\widehat {\mathcal{U}}(\xi, t)=&P\begin{pmatrix} \text{e}^{\lambda_1t}&0&0&0\\ 0&\text{e}^{\lambda_2t}&0&0\\ 0&0&\text{e}^{\lambda_3t}&t\text{e}^{\lambda_3t}\\ 0&0&0&\text{e}^{\lambda_3t} \end{pmatrix}P^{-1}\widehat {\mathcal{U}}_0(\xi)\\=&\left(\text{e}^{\lambda_1t}P^1(\xi)+\text{e}^{\lambda_2t}P^2(\xi)+\text{e}^{\lambda_3t}P^3(\xi)+t\text{e}^{\lambda_3t} P^4(\xi)\right)\widehat {\mathcal{U}}_0(\xi). \end{align}\tag{62}\] Then, after a direct and complicated computation, we can get the explicit expressions of \(P^i(\xi)\): \[\begin{align}\label{24631}P^1(\xi)=\bar{P}^2(\xi) =&\begin{pmatrix} \frac{\beta_1}{2(\beta_1+\beta_4)}&\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)} &\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ -\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_1}{2(\beta_1+\beta_4)}&-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_2}{2(\beta_1+\beta_4)}\\ \frac{\beta_2}{2(\beta_1+\beta_4)}&\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}&\frac{\beta_4}{2(\beta_1+\beta_4)} &\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{3}{2}}\mathrm{i}\\ -\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}&\frac{\beta_2}{2(\beta_1+\beta_4)}&-\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{1}{2}}\mathrm{i}& \frac{\beta_4}{2(\beta_1+\beta_4)} \end{pmatrix}+\mathcal{O}(|\xi|), \end{align}\tag{63}\] \[\begin{align} \label{24634}P^3(\xi)=&\begin{pmatrix} \frac{\beta_4}{\beta_1+\beta_4}&0&-\frac{\beta_2}{\beta_1+\beta_4} &0\\ 0&\frac{\beta_4}{\beta_1+\beta_4}&0& -\frac{\beta_2}{\beta_1+\beta_4}\\ -\frac{\beta_2}{\beta_1+\beta_4}&0 &\frac{\beta_1}{\beta_1+\beta_4} &0\\ 0&-\frac{\beta_2}{\beta_1+\beta_4}&0& \frac{\beta_1}{\beta_1+\beta_4} \end{pmatrix}+\mathcal{O}(|\xi|), \end{align}\tag{64}\] \[\begin{align} \label{24633} P^4(\xi)=&\begin{pmatrix} \frac{\beta_1\beta_4(\nu^+-\nu^-)}{2(\beta_1+\beta_4)^2}+\frac{\beta_4\nu^+}{2(\beta_1+\beta_4)}&-\frac{\beta_4}{(\beta_1+\beta_4)|\xi|}&\frac{\beta_2\beta_4(\nu^+-\nu^-)}{2(\beta_1+\beta_4)^2} -\frac{\beta_2\nu^-}{2(\beta_1+\beta_4)} &\frac{\beta_2}{(\beta_1+\beta_4)|\xi|}\\ 0&-\frac{\beta_4(\beta_4\nu^++\beta_1\nu^-)}{2(\beta_1+\beta_4)^2}&0& \frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{2(\beta_1+\beta_4)^2}\\ \frac{\beta_1\beta_2(\nu^--\nu^+)}{2(\beta_1+\beta_4)^2}-\frac{\beta_2\nu^+}{2(\beta_1+\beta_4)}&\frac{\beta_2}{(\beta_1+\beta_4)|\xi|} &\frac{\beta_1\beta_4(\nu^--\nu^+)}{2(\beta_1+\beta_4)^2}+\frac{\beta_1\nu^-}{2(\beta_1+\beta_4)} &-\frac{\beta_1}{(\beta_1+\beta_4)|\xi|}\\ 0&\frac{\beta_2(\beta_4\nu^++\beta_1\nu^-)}{2(\beta_1+\beta_4)^2}&0& -\frac{\beta_1(\beta_4\nu^++\beta_1\nu^-)}{2(\beta_1+\beta_4)^2} \end{pmatrix}|\xi|^2\nonumber\\ &\quad+\mathcal{O}(|\xi|^3). \end{align}\tag{65}\]

Then one can see that the two cases above are almost the same in the sense of deriving the space-time behavior of Green’s function \(\tilde{G}(x,t)\), and hence, we only focus on the first case in the following without loss of generality.

High frequency part.

From the characteristic equation (57 ), we have the following expansion of spectra in the high frequency.

After a direct and tedious computation, we can get the following expansion of the projection \(P^i\) with \(i=1,2,3,4\) in the high frequency according to Lemma 2: \[\begin{align} \label{24636} P^{i}(\xi)=&\begin{pmatrix} \mathcal{O}(1)&\mathcal{O}(|\xi|^{-1})&\mathcal{O}(|\xi|^{-2})&\mathcal{O}(|\xi|^{-3})\\ \mathcal{O}(|\xi|)&\mathcal{O}(1)&\mathcal{O}(|\xi|^{-1})&\mathcal{O}(|\xi|^{-2})\\ \mathcal{O}(|\xi|^{-2})&\mathcal{O}(|\xi|^{-3})&\mathcal{O}(1)&\mathcal{O}(|\xi|^{-1})\\ \mathcal{O}(|\xi|^{-1})&\mathcal{O}(|\xi|^{-2})&\mathcal{O}(|\xi|)&\mathcal{O}(1) \end{pmatrix}+\cdots, \end{align}\tag{66}\] which together with (56 ) and Lemma 2) gives that

Middle frequency part.

For the middle frequency part, we only need to prove the following two lemmas. The first one is about the positivity of the real part of spectra based on the Routh-Hurwitz theorem.

The analyticity for the middle frequency is partially based on the idea in Li [45].

Lemma 4 and Lemma 5 directly yield the space-time estimates for middle frequency \(\tilde{G}^m\):

3 Pointwise space-time behavior of Green’s function↩︎

We shall bring Green’s function in Fourier space back to the physical space, and combine the compressible part and the incompressible part to obtain the pointwise description of Green’s function of the original system (45 ). The difficulties mainly include resolving the singularity from the Riesz operator in low frequency of Green’s matrix by using suitable combinations, and giving the description of the singular part in high frequency arising from the definition of Green’s function.

3.1 Pointwise space-time behavior of low frequency↩︎

The Hodge decomposition, Lemma 1 together with (53 ) and (59 )-(61 ) directly give the following: \[\begin{align} \label{3461} \hat{n}^{+,l}=&\Big[\frac{\beta_1(e^{\lambda_1t}+e^{\lambda_2t})}{2(\beta_1+\beta_4)}+\Big(\frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}+\Big(\frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big]\hat{n}_0^+\nonumber\\ \!\!&\!+\!\Big[\mathrm{i}\frac{\beta_1(e^{\lambda_1t}\!-\!e^{\lambda_2t})}{2(\beta_1+\beta_4)^\frac{3}{2}}-\frac{\beta_4}{\mathcal{R}|\xi|}\Big(e^{\lambda_3t}+e^{\lambda_4t}\Big)\Big]\frac{\mathrm{i}\xi^T}{|\xi|}\hat{\mathbf{m}}_{0}^{+}\nonumber\\ &+\Big[\frac{\beta_2(e^{\lambda_1t}\!+\!e^{\lambda_2t})}{2(\beta_1+\beta_4)}+\Big(\frac{\beta_2\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}+\Big(\frac{\beta_2\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big]\hat{n}_0^{-} \nonumber\\[-0.5mm] &+\Big[\mathrm{i}\frac{\beta_2(e^{\lambda_1t}-e^{\lambda_2t})}{2(\beta_1+\beta_4)^\frac{3}{2}}+\frac{\beta_2}{\mathcal{R}|\xi|}(e^{\lambda_3t}+e^{\lambda_4t})\Big]\frac{\mathrm{i}\xi^T}{|\xi|}\hat{\mathbf{m}}_{0}^{-}+\cdots\nonumber\\[-0.5mm] :=&(G_{11}^l,G_{12}^l,G_{13}^l,G_{14}^l)\cdot(\hat{n}_{0}^+,\hat{\mathbf{m}}_0^+,\hat{n}_{0}^-,\hat{\mathbf{m}}_{0}^-), \end{align}\tag{67}\] \[\begin{align} \label{3462} \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!&\!\!\!\!\!\!\hat{\mathbf{m}}^{+,l} =-\widehat{\Lambda^{-1}\nabla \varphi^+}-\widehat{\Lambda ^{-1}{\rm div} \phi^+}\nonumber\\ =&\Big[-\frac{\beta_1}{2(\beta_1+\beta_4)^\frac{1}{2}}(e^{\lambda_1t}-e^{\lambda_2t})\Big]\frac{\mathrm{i}\xi\hat{n}_{0}^{+}}{|\xi|}\nonumber\\ &+\bigg[\Big(\frac{\beta_1(e^{\lambda_1t}\!+\!e^{\lambda_2t})}{2(\beta_1+\beta_4)}-\Big(\frac{\beta_4(\beta_4\nu^+\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_4\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}-\Big(\frac{\beta_4(\beta_4\nu^+\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big)\frac{\xi\xi^T}{|\xi|^{2}}\nonumber\\ &\;\;\;\;\;\;\;\;\;\;+e^{-\nu^+|\xi|^2t}\Big(I-\frac{\xi\xi^T}{|\xi|^{2}}\Big)\bigg]\hat{\mathbf{m}}_{0}^+\nonumber\\ &+\Big[-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}(e^{\lambda_1t}+e^{\lambda_2t})\mathrm{i}\xi\Big]\hat{n}_{0}^{-}\nonumber\\[-0.5mm] &\!+\!\bigg[\frac{\beta_2(e^{\lambda_1t}\!+\!e^{\lambda_2t})}{2(\beta_1+\beta_4)}\!+\!\Big(\frac{\beta_2(\beta_4\nu^+\!\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}\!+\!\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}\!+\!\Big(\frac{\beta_2(\beta_4\nu^+\!\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}\!+\!\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\bigg]\frac{\xi\xi^T}{|\xi|^{2}}\hat{\mathbf{m}}_{0}^{-}\!+\cdots\nonumber\\ :=&(G_{21}^l,G_{22}^l,G_{23}^l,G_{24}^l)\cdot(\hat{n}_{0}^+,\hat{\mathbf{m}}_0^+,\hat{n}_{0}^-,\hat{\mathbf{m}}_{0}^-), \end{align}\tag{68}\] \[\begin{align} \label{3463} \hat{n}^{-,l}=&\Big[\frac{\beta_2(e^{\lambda_1t}+e^{\lambda_2t})}{2(\beta_1+\beta_4)}+\Big(\frac{\beta_1\beta_2(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_4\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}+\Big(\frac{\beta_1\beta_4(\nu^+-\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_4\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big]\hat{n}_0^+\nonumber\\ \!\!&\!+\!\Big[\mathrm{i}\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{3}{2}}(e^{\lambda_1t}\!-\!e^{\lambda_2t})+\frac{\beta_2}{\mathcal{R}|\xi|}\Big(e^{\lambda_3t}+e^{\lambda_4t}\Big)\Big]\frac{\mathrm{i}\xi^T}{|\xi|}\hat{\mathbf{m}}_{0}^{+}\nonumber\\ &+\Big[\frac{\beta_4(e^{\lambda_1t}\!+\!e^{\lambda_2t})}{2(\beta_1+\beta_4)}+\Big(\frac{\beta_1\beta_4(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}+\Big(\frac{\beta_1\beta_4(\nu^--\nu^+)}{(\beta_1+\beta_4)\mathcal{R}}-\frac{\beta_1\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big]\hat{n}_0^{-} \nonumber\\[-0.5mm] &+\Big[\mathrm{i}\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{3}{2}}(e^{\lambda_1t}-e^{\lambda_2t})-\frac{\beta_1}{\mathcal{R}|\xi|}(e^{\lambda_3t}+e^{\lambda_4t})\Big]\frac{\mathrm{i}\xi^T}{|\xi|}\hat{\mathbf{m}}_{0}^{-}+\cdots\nonumber\\[-0.5mm] :=&(G_{11}^l,G_{12}^l,G_{13}^l,G_{14}^l)\cdot(\hat{n}_{0}^+,\hat{\mathbf{m}}_0^+,\hat{n}_{0}^-,\hat{\mathbf{m}}_{0}^-), \end{align}\tag{69}\] \[\begin{align} \label{3464} &\!\!\!\!\!\!\!\!\!\hat{\mathbf{m}}^{-,l}=-\widehat{\Lambda^{-1}\nabla \varphi^-}-\widehat{\Lambda ^{-1}{\rm div} \phi^-}\nonumber\\ =&\Big[-\frac{\beta_2}{2(\beta_1+\beta_4)^\frac{1}{2}}(e^{\lambda_1t}-e^{\lambda_2t})\Big]\frac{\mathrm{i}\xi\hat{n}_{0}^{+}}{|\xi|}\nonumber\\ &+\bigg[\Big(\frac{\beta_2(e^{\lambda_1t}\!+\!e^{\lambda_2t})}{2(\beta_1+\beta_4)}+\Big(\frac{\beta_2(\beta_4\nu^+\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}+\Big(\frac{\beta_2(\beta_4\nu^+\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big)\frac{\xi\xi^T}{|\xi|^{2}}\bigg]\hat{\mathbf{m}}_{0}^+\nonumber\\ &+\Big[-\frac{\beta_4}{2(\beta_1+\beta_4)^\frac{1}{2}}(e^{\lambda_1t}+e^{\lambda_2t})\mathrm{i}\xi\Big]\hat{n}_{0}^{-}\nonumber\\[-0.5mm] &+\bigg[\Big(\frac{\beta_4(e^{\lambda_1t}\!+\!e^{\lambda_2t})}{2(\beta_1+\beta_4)}\!-\!\Big(\frac{\beta_1(\beta_4\nu^+\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_4}{\mathcal{R}}\Big)e^{\lambda_3t}-\Big(\frac{\beta_1(\beta_4\nu^+\!+\!\beta_1\nu^-)}{(\beta_1+\beta_4)\mathcal{R}}+\frac{\beta_2\tilde{\lambda}_3}{\mathcal{R}}\Big)e^{\lambda_4t}\Big)\frac{\xi\xi^T}{|\xi|^{2}}\nonumber\\ &\;\;\;\;\;\;\;+e^{-\nu^-|\xi|^2t}\Big(I-\frac{\xi\xi^T}{|\xi|^{2}}\Big)\bigg]\hat{\mathbf{m}}_{0}^{-}+\cdots\nonumber\\ :=&(G_{21}^l,G_{22}^l,G_{23}^l,G_{24}^l)\cdot(\hat{n}_{0}^+,\hat{\mathbf{m}}_0^+,\hat{n}_{0}^-,\hat{\mathbf{m}}_{0}^-). \end{align}\tag{70}\] We can mainly deal with the leading term of each entry, since the rest terms (denoted by “\(\cdots\)") in (67 )-(70 )) are analytic and have faster temporal decay rates than the leading terms. The main difficulty focuses on the interaction of the singular Riesz operators (the first and second-order) in low frequency, the wave operators and the heat kernels. We merely take a few typical leading terms for instances to elaborate on the avoidance of the seeming singularity at \(\xi=0\) and derive different wave patterns in the low frequency based on suitable combinations and cancellations.

To facilitate explanation, we denote \(\hat{\mathbf{w}}_t=\cos({\rm c}|\xi|t)\) and \(\hat{\mathbf{w}}=\frac{\sin({\rm c}|\xi|t)}{|\xi|}\) are the Fourier transform of wave operators, where \({\rm c}\) is the propagation speed of the bulk in Lemma 1. Then, we first reformulate the following typical term in \(\hat{G}_{12}^l\): \[\label{3465} \begin{array}{rl} &\frac{e^{\lambda_1t}-e^{\lambda_2t}}{2}\frac{\xi^T}{|\xi|}=\frac{\xi^T}{|\xi|}e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t}\sin({\rm c}|\xi|t+|\xi|\beta(|\xi|^2)t)\\[2mm] =&\xi^T e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t}\Big\{\sin({\rm c}|\xi|t)\cos(|\xi|\beta(|\xi|^2)t)+\cos({\rm c}|\xi|t)\sin(|\xi|\beta(|\xi|^2)t)\Big\}\\[2mm] =&\xi^T\hat{\mathbf{w}}e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t}+\xi^T\hat{\mathbf{w}}(\cos(|\xi|\beta(|\xi|^2)t)-1)e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t}\\[2mm] &+\xi^T\hat{\mathbf{w}}_t\frac{\sin(|\xi|\beta(|\xi|^2)t)}{|\xi|}e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t}, \end{array}\tag{71}\] where \(b_1=\frac{\beta_1\nu^++\beta_4\nu^-}{2(\beta_1+\beta_4)}>0\) and \(\beta(\cdot)\) is an analytic function.

Recall the estimates given in Li [45] for any \(N>0\) that: \[\begin{array}{rl} &\Big|\mathcal{F}^{-1}\Big((\cos(|\xi|\beta(|\xi|^2)t)-1)e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t}\Big)\Big|+\Big|\mathcal{F}^{-1}(\frac{\sin(|\xi|\beta(|\xi|^2)t)}{|\xi|}e^{-b_1|\xi|^2t+\mathcal{O}(|\xi|^4)t})\Big|\\[2mm] \lesssim &(1+t)^{-\frac{3}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-N}, \end{array}\] which together with Kirchhoff formula in dimension three give that \[\label{346540041} \Big|\partial_x^\alpha\mathcal{F}^{-1}\Big(\frac{e^{\lambda_1t}-e^{\lambda_2t}}{2}\frac{\xi^T}{|\xi|}\Big)\Big|\lesssim (1+t)^{-\frac{4+|\alpha|}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}.\tag{72}\]

Next, we study the first typical entry \(G_{22}\), which is different from that of the one-phase fluid model (the compressible Navier-Stokes model) in [21], [22], where \[\label{346540141} \begin{array}{rl} \hat{G}_{22}^l=&\frac{e^{\lambda_1t}+e^{\lambda_2t}}{2}\frac{\xi\xi^T}{|\xi|^{2}}-e^{-|\xi|^2t}\frac{\xi\xi^T}{|\xi|^{2}}+\cdots\\[2mm] \sim &\underbrace{(\hat{\mathbf{w}}_t-1)\frac{\xi\xi^T}{|\xi|^{2}}e^{-b_1|\xi|^2t}}_{{\rm Riesz\;wave\mathbin{\text{-}}I}}+\underbrace{\frac{\xi\xi^T}{|\xi|^{2}}(e^{-b_1|\xi|^2t}-e^{-|\xi|^2t})}_{{\rm Riesz\;wave\mathbin{\text{-}}I\!I}}+\cdots. \end{array}\tag{73}\] We use \(\hat{R}_i\) with \(i=1,2\) to denote the Riesz wave-I and Riesz wave-II, respectively. As in Liu-Noh [22] and Li [45] for the isentropic and non-isentropic compressible Navier-Stokes system, their pointwise space-time descriptions hold for an arbitrarily large integer \(N\) that \[\begin{align} \label{3466} \begin{array}{rl} &|\partial_x^\alpha(\mathcal{F}^{-1}(\hat{R_1}))|\lesssim (1+t)^{-\frac{3+|\alpha|}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-\frac{3+|\alpha|}{2}}+(1+t)^{-\frac{4+|\alpha|}{2}}\big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\big)^{-N},\\ &|\partial_x^\alpha(\mathcal{F}^{-1}(\hat{R_2}))|\lesssim (1+t)^{-\frac{3+|\alpha|}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-N}. \end{array} \end{align}\tag{74}\]

Due to the complicity of the two-phase fluid model, its entry \(G_{22}\) in low frequency only can be rewritten as \[\label{3467} \begin{array}{rl} \hat{G}_{22}^l=&a_1\frac{e^{\lambda_1t}+e^{\lambda_2t}}{2}\frac{\xi\xi^T}{|\xi|^{2}}-a_2e^{-|\xi|^2t}\frac{\xi\xi^T}{|\xi|^{2}}+\cdots\;\;\;\;\;\;\;\;\;\;(a_1\neq a_2)\\[2mm] \sim &\underbrace{(\hat{\mathbf{w}}_t-1)\frac{\xi\xi^T}{|\xi|^{2}}e^{-b_1|\xi|^2t}}_{{\rm Riesz\;wave\mathbin{\text{-}}I}}+\underbrace{\frac{\xi\xi^T}{|\xi|^{2}}e^{-|\xi|^2t}}_{{\rm Riesz\;wave\mathbin{\text{-}}I\!I\!I}}+\cdots, \end{array}\tag{75}\] and we also use \(\hat{R}_3\) to denote the Riesz wave-III. Although Riesz wave-III is different from Riesz wave-II, we can also get the following according Lemma 12: \[\begin{align} \label{3468} \begin{array}{rl} &|\partial_x^\alpha(\mathcal{F}^{-1}(\hat{R_3}))|\lesssim (1+t)^{-\frac{3+|\alpha|}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-\frac{3+|\alpha|}{2}}. \end{array} \end{align}\tag{76}\]

Next, we turn to the Riesz wave-IV (like \(|\xi|^{-1}e^{-|\xi|^2t}\) in Fourier space) in \(G_{12}^l,G_{14}^l,G_{32}^l,G_{34}^l\), which is denoted by \(R_4\) and has been studied by the first author and the second author in the present paper for the unipolar compressible Navier-Stokes-Poisson system in [40]. In particular, \[\begin{align} \label{3469} \begin{array}{rl} &|\partial_x^\alpha(\mathcal{F}^{-1}(\hat{R_4}))|\lesssim (1+t)^{-\frac{2+|\alpha|}{2}}\big(1+\frac{|x|^2}{1+t}\big)^{-\frac{2+|\alpha|}{2}}. \end{array} \end{align}\tag{77}\] One can see that the Riesz wave-IV has the worst decay for both temporal and spacial variables. Different from the difficulty in [40], where its Green’s function does not contain the Huygens waves, that is, the interaction of the Riesz wave-IV and the Huygens wave will not happen. However, due to the complication of the two-phase fluid model (1 ), we have to deal with the nonlinear coupling of these two waves (at least). In fact, one of the large difficulty in the present paper is to derive its sharp estimate. See the details in the proof of \(K_4\) in Lemma 11.

The second difficulty is arising from the non-conservative structure of the system (1 ) due to two pressure terms and the presence of the Huygens waves in the low frequency of Green’s function. In fact, in existing research works in this direction rely heavily on the nonlinear convolutions in Lemma 10, and these estimates are sharp in some sense. On the other hand, one can hardly find two kinds of linear combination to compensate this non-conservation as the previous work in this direction for the other related compressible fluid models. To achieve the relatively sharp pointwise space-time description of the solution to the nonlinear problem, we have to construct a series of new nonlinear convolution estimates of different wave patterns. See the details in Lemma 11.

In a conclusion, we have the pointwise estimate for Green’s function in the low frequency:

We emphasize again that the estimate (?? ) is crucial for us to gain the desired, sharp estimate concerning the nonlinear convolution of the Riesz wave-IV and the Huygens wave in Lemma 11.

Combining Lemma 7, Remark 8 and Lemma 6, we conclude the following.

4 Pointwise estimates for nonlinear system↩︎

In this section, by using the representation of the solution, we derive the pointwise estimates for the solution of the nonlinear system based on the pointwise space-time description of Green’s function and the convolution estimates of nonlinear coupling.

Without loss of generality, we use \(\partial_x^k\) to denote \(\partial_x^\alpha\) with \(|\alpha|=k\). By using Duhamel principle, we can get the representation of the solution \((n^+,\mathbf{m}^+,n^-,\mathbf{m}^-)\) for the nonlinear problem (43 ): \[\begin{align} \label{5461} &&\partial_x^k(n^+,\mathbf{m}^+,n^-,\mathbf{m}^-)^T\nonumber\\ &=&\partial_x^k G(x,t)\ast_x(n^+_0,\mathbf{m}^+_0,n^-_0,\mathbf{m}^-_0)^T +\int_0^t\partial_x^k G(\cdot,t-\tau)\ast_x(0,F_1,0,F_2)^T(\cdot,\tau)d\tau, \end{align}\tag{78}\] where the nonlinear terms \(F_1\) and \(F_2\) are defined in (45 ).

Initial propagation. Use \((\breve{n}^+,\breve{\mathbf{m}}^+,\breve{n}^-,\breve{\mathbf{m}}^-)\) to denote the linear part of the solution in (78 ). According to the representation (78 ), Proposition 1, the initial condition (2 ) and the convolution estimate in Lemma 10 for the initial propagation, we can immediately get the following: \[\label{5462} \begin{array}{rl} &\left|\partial_x^k\left(\!\! \begin{array}{ccc} \breve{n}^+ \\ \breve{n}^- \end{array} \!\!\right)\right|=\!\left|\partial_x^k\!\left(\! \begin{array}{cccc} G_{11} & G_{12} & G_{13} &G_{14}\\ G_{31} & G_{32} & G_{33} &G_{34}\\ \end{array} \!\right)\!\!\ast_x(n^+_0,\mathbf{m}^+_0,n^-_0,\mathbf{m}^-_0)^T \right|\\ \leq\!\! & \left|\partial_x^k\!\bigg(\! \begin{array}{cccc} G_{11}-G_{S1} & G_{12}-G_{S1} & G_{13}-G_{S1} & G_{14}-G_{S1}\\ G_{31}-G_{S1} & G_{32}-G_{S1} & G_{33}-G_{S1} & G_{34}-G_{S2}\\ \end{array} \!\bigg)\!\!\ast_x (n^+_0,\mathbf{m}^+_0,n^-_0,\mathbf{m}^-_0)^T\right|\\ &+\left| \left(\! \begin{array}{cccc} G_{S1} &G_{S1} & G_{S1} & G_{S1}\\ G_{S1} & G_{S1} & G_{S1} & G_{S1}\\ \end{array} \!\right)\!\ast_x\partial_x^k (n^+_0,\mathbf{m}^+_0,n^-_0,\mathbf{m}^-_0)^T\right|\\ \leq &\displaystyle C\epsilon\bigg((1+t)^{-\frac{2+k}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{2+k}{2}} +(1+t)^{-\frac{4+k}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg), \end{array}\tag{79}\] \[\label{5463} \begin{array}{rl} &\left|\partial_x^k\left(\!\! \begin{array}{ccc} \breve{\mathbf{m}}^+ \\ \breve{\mathbf{m}}^- \end{array} \!\!\right)\right|=\!\left|\partial_x^k\!\left(\! \begin{array}{cccc} G_{21} & G_{22} & G_{23} &G_{24}\\ G_{41} & G_{42} & G_{43} &G_{44}\\ \end{array} \!\right)\!\!\ast_x (n^+_0,\mathbf{m}^+_0,n^-_0,\mathbf{m}^-_0)^T\right|\\ \leq\!\! & \left|\partial_x^k\!\left(\! \begin{array}{cccc} G_{21}-G_{S2} & G_{22}-G_{S1} & G_{23}-G_{S1} &G_{24}-G_{S1}\\ G_{41}-G_{S1} & G_{42}-G_{S1} & G_{43}-G_{S2} &G_{44}-G_{S1}\\ \end{array} \!\right)\!\!\ast_x (n^+_0,\mathbf{m}^+_0,n^-_0,\mathbf{m}^-_0)^T\right|\\ &+\left| \bigg(\! \begin{array}{cccc} \uwave{G_{S2}} & G_{S1} & G_{S1} & G_{S1}\\ G_{S1} & G_{S1} & \uwave{G_{S2}} &G_{S1} \end{array} \!\bigg)\!\ast_x(\uwave{\partial_x^{k+1}n_{0}^+},\partial_x^k \mathbf{m}_0^+,\uwave{\partial_x^{k+1}n_{0}^-},\partial_x^k \mathbf{m}_{0}^-)^T\right|\\ \leq &\displaystyle C\epsilon\bigg((1+t)^{-\frac{3+k}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}} +(1+t)^{-\frac{4+k}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg). \end{array}\tag{80}\] Note that it requires one higher order regularity assumption on the initial value of the fraction densities \((n_{0}^+,n_{0}^-)\), which is consistent with the standard energy estimates for the compressible Naiver-Stokes-Korteweg system.

Nonlinear Coupling. According to the above initial propagation, we give the ansatz for the nonlinear problem with \(k\leq3\): \[\begin{align} \label{5464} |(n^+,n^-)|\leq 2C\epsilon\bigg\{(1+t)^{-1}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-\frac{3}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg\},\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\; \end{align}\tag{81}\] \[\begin{align} \label{5465} &&|\bar\rho^-n^++\bar\rho^+n^-|+|\partial_x^k(n^+,n^-)|\nonumber\\ &&\;\;\;\;\;\;\;\;\;\;\;\leq 2C\epsilon\bigg\{(1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg\},\;\;k=1,2,3, \;\;\;\;\;\;\;\;\;\; \end{align}\tag{82}\] \[\label{5466} |\partial_x^k (\mathbf{m}^+,\mathbf{m}^-)|\leq 2C\epsilon\bigg\{(1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg\},\;k=0,1,2,3.\;\;\tag{83}\]

To resolve this difficulty, we have to find out a refined decay estimate of the nonlinear terms \(Q_1=(\beta_1\nabla n^++\beta_2\nabla n^-)-\alpha^+\nabla P\) and \(Q_2=(\beta_3\nabla n^++\beta_4\nabla n^-)-\alpha^-\nabla P\) in \(F_1\) and \(F_2\) arising from two pressure terms. To this end, after a careful observation and reformulation, we can have \[\begin{align} Q_1&=&\mathcal{C}^{2}(n^+\!+\!1,n^-\!+\!1)n^+\Big(\frac{\rho^-}{\rho^+}\nabla n^++\nabla n^-\Big)\nonumber\\ &&+\Big(\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)\rho^-}{\rho^+}-\frac{\mathcal{C}^2(1,1)\bar\rho^-}{\bar\rho^+}\Big)\nabla n^++(\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)-\mathcal{C}^2(1,1))\nabla n^-\nonumber\\ &=&\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)}{\rho^+}n^+\bigg[(\rho^--\bar\rho^-)\nabla n^+\!+\!(\rho^+-\bar\rho^+)\nabla n^-\bigg]\nonumber\\ &&+\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)}{\rho^+}n^+(\uwave{\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-})\nonumber\\ &&+\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)-\mathcal{C}^2(1,1)}{\rho^+}[(\rho^--\bar\rho^-)\nabla n^++(\rho^+-\bar\rho^+)\nabla n^-]\nonumber\\ &&+\frac{1}{\rho^+}(\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)-\mathcal{C}^2(1,1))(\uwave{\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-})\nonumber\\ &&+\mathcal{C}^2(1,1)\bigg(\frac{1}{\rho^+}-\frac{1}{\bar\rho^+}\bigg)(\uwave{\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-})\nonumber\\ &&+\mathcal{C}^2(1,1)\bigg(\frac{1}{\rho^+}-\frac{1}{\bar\rho^+}\bigg)[(\rho^--\bar\rho^-)\nabla n^++(\rho^+-\bar\rho^+)\nabla n^-]\nonumber\\ &&+\frac{\mathcal{C}^2(1,1)}{\bar\rho^+}[(\uwave{\rho^--\bar\rho^-})\nabla n^++(\uwave{\rho^+-\bar\rho^+})\nabla n^-]:=\sum\limits_{i=1}^7 \mathcal{S}_i,\tag{84} \tag{85} \end{align}\] where \(\rho^\pm=\rho^\pm(n^++1,n^-+1)\) and \(\bar\rho^\pm=\rho^\pm(1,1)\). It is obvious that \(\mathcal{S}_1\), \(\mathcal{S}_3\) and \(\mathcal{S}_6\) are the cube of the unknowns and naturally have faster decay from the ansatz. The term \(\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-\) in \(\mathcal{S}_2\), \(\mathcal{S}_4\) and \(\mathcal{S}_5\) has an additional \((1+t)^{-\frac{1}{2}}\)-decay rate than \(\nabla n^\pm\) according to the ansatz (82 ). Finally, since \(\rho^\pm-\bar\rho^\pm\sim \bar\rho^-\nabla n^++\bar\rho^+\nabla n^-\) due to the relation between the densities \(\rho^\pm\) and the fraction densities \(n^\pm\), it also has an additional \((1+t)^{-\frac{1}{2}}\)-decay rate than \(n^\pm\) due to the corresponding estimate in ?? .

Similarly, one has \[\begin{align} Q_2&=&\mathcal{C}^{2}(n^+\!+\!1,n^-\!+\!1)n^-\Big(\frac{\rho^+}{\rho^-}\nabla n^-+\nabla n^+\Big)\nonumber\\ &&+\Big(\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)\rho^+}{\rho^-}-\frac{\mathcal{C}^2(1,1)\bar\rho^+}{\bar\rho^-}\Big)\nabla n^-+(\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)-\mathcal{C}^2(1,1))\nabla n^+\nonumber\\ &=&\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)}{\rho^-}n^-\bigg[(\rho^--\bar\rho^-)\nabla n^-\!+\!(\rho^+-\bar\rho^+)\nabla n^+\bigg]\nonumber\\ &&+\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)}{\rho^-}n^-(\uwave{\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-})\nonumber\\ &&+\frac{\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)-\mathcal{C}^2(1,1)}{\rho^-}[(\rho^--\bar\rho^-)\nabla n^-+(\rho^+-\bar\rho^+)\nabla n^+]\nonumber\\ &&+\frac{1}{\rho^+}(\mathcal{C}^2(n^+\!+\!1,n^-\!+\!1)-\mathcal{C}^2(1,1))(\uwave{\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-})\nonumber\\ &&+\mathcal{C}^2(1,1)\bigg(\frac{1}{\rho^-}-\frac{1}{\bar\rho^-}\bigg)(\uwave{\bar\rho^-\nabla n^++\bar\rho^+\nabla n^-})\nonumber\\ &&+\mathcal{C}^2(1,1)\bigg(\frac{1}{\rho^-}-\frac{1}{\bar\rho^-}\bigg)[(\rho^--\bar\rho^-)\nabla n^-+(\rho^+-\bar\rho^+)\nabla n^+]\nonumber\\ &&+\frac{\mathcal{C}^2(1,1)}{\bar\rho^-}[(\uwave{\rho^--\bar\rho^-})\nabla n^-+(\uwave{\rho^+-\bar\rho^+})\nabla n^+]:=\sum\limits_{i=8}^{14} \mathcal{S}_i.\tag{86} \tag{87} \end{align}\]

As a result, the convolution of the Riesz wave-IV and the nonlinear terms about the pressure terms \(\alpha^\pm\nabla P\) has the following refine form instead of that in (?? ): \[\begin{align} \mathcal{N}_1= \int_0^t\!\int_{\mathbb{R}^3}(1+t-\tau)^{-1}\Big(1+\frac{|x-y|^2}{1+t-\tau}\Big)^{-1}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau.\label{54610} \end{align}\tag{88}\] Based on the estimate \(K_4\) established in Lemma 11, one has \[\begin{align} \mathcal{N}_1\lesssim (1+t)^{-1}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-\frac{3}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}.\label{54611} \end{align}\tag{89}\]

On the other hand, we consider another difficulty in dealing with the interaction of the Huygens wave in Green’s function and the Huygens waves in nonlinear terms for the non-conservative fluid models. The estimate is totally new, which is given in \(K_7\) in Lemma 11. In particular, we have for arbitrarily large constant \(N>0\) that \[\begin{align} \mathcal{N}_2=&\int_{0}^{t}\!\int_{\mathbb{R}^3}(1+t-\tau)^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau\nonumber\\[3mm] \lesssim &\;(1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}+(1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}.\label{54612} \end{align}\tag{90}\] This is different from the classical estimate for the conservative system given in [21], [22]: \[\begin{align} &&\int_{0}^{t}\int_{\mathbb{R}^3}\uwave{(1+t-\tau)^{-\frac{5}{2}}}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1\!+\!t\!-\!\tau}\Big)^{-N}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1\!+\!\tau}\Big)^{-3}dyd\tau\nonumber\\ &\lesssim& (1+t)^{-2}\Big(\big(1+\frac{|x|^2}{1+t}\big)^{-\frac{3}{2}}+\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-\frac{3}{2}}\Big),\label{5461240141} \end{align}\tag{91}\] where the extra \((1+t)^{-\frac{1}{2}}\)-decay of Green’s function (the underlined term in (91 )) is arising from the conservation of the nonlinear fluid model studied in the previous work [21], [22], [41].

Going back to the estimate (90 ), we have to overcome the difficulty in the subregion space-time estimates caused by the deteriorated time decay of Green’s function in the convolution due to the non-conservation, and simultaneously ensure that the obtained pointwise estimates are optimal in terms of time decay (consistent with the \(L^2\)-decay rate derived by combining the energy method with spectral analysis), we have carried out more elaborate calculations in the subregion estimates. In particular, we have fully exploited the interchange between temporal and spatial variables within specific regions, and thus obtained relatively sharp pointwise space-time estimates. See the details in the proof of \(K_7\) in Appendix.

Finally, we shall close the ansatz on the solution of the nonlinear problem by substituting (81 )-(83 ) into the representation of the solution \((n^+,\mathbf{m}^+,n^-,\mathbf{m}^-)\) in (67 )-(68 ) and using the aforementioned nonlinear convolution estimates. Note that the singular term \(G_{S2}\) in \(G_{21}\) and \(G_{43}\) correspond to the nonlinear term “0" in the two continuity equations of the system (1 ). Hence, in the following, we only focus on the unique singular term \(G_{S1}\) when dealing with nonlinear coupling. In particular, we have to split Green’s function into the regular term \(G_{ij}-G_{S1}\) and the singular term \(G_{S1}\). Then, one can put all of the derivatives on \(G_{ij}-G_{S1}\) to deal with the convolution between the regular part \(G_{ij}-G_{S1}\) and the nonlinear terms, however, one have to put all of the derivatives on the nonlinear terms to deal with the convolution of the singular part and the nonlinear terms.

We just take nonlinear coupling of the fraction densities \(n^+\) for example, since that of \(n^-\) is the same and those of \(\mathbf{m}^\pm\) can be estimated similarly. To facilitate the description, we use \(\tilde{n}^+\) to denote the nonlinear part of \(n^+\) in (78 ). Then by using Proposition 1 and the nonlinear convolution estimates in Lemma 11, one has for \(0\leq k\leq2\) that \[\begin{align} \label{xuqtpvem} |\partial_x^k\tilde{n}^+|\leq &\int_0^t|\partial_x^k(G_{12}-G_{S1})(\cdot,t-\tau)\ast F_1(\cdot,\tau)|d\tau+\int_0^t|\partial_x^k(G_{14}-G_{S1})(\cdot,t-\tau)\ast F_2(\cdot,\tau)|d\tau\nonumber\\ &+\int_0^t|G_{S1}(\cdot,t-\tau)\ast \partial_x^k (F_1(\cdot,\tau),F_2(\cdot,\tau))|d\tau\nonumber\\ \lesssim &\;2C\epsilon\bigg((1+t)^{-1}\Big(1+\frac{|x|^2}{1+t}\Big)^{-1} +(1+t)^{-\frac{3}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}\bigg)\nonumber\\ &+\int_0^t|G_{S1}(\cdot,t-\tau)\ast \partial_x^k (n^+\partial_x^3n^++\partial_x^2(n^+\mathbf{m}^+)+\cdots)(\cdot,\tau)|d\tau\nonumber. \end{align}\tag{92}\] When \(k=0\), we need the pointwise ansatz of \(\partial_x^3 n^\pm\) and \(\partial_x^2\mathbf{m}^\pm\) according to the nonlinear terms and the nonlinear convolution in Lemma 11. Then when \(k=3\), it also requires the \(L^\infty\)-estimate of \(\partial_x^6n^\pm\) and \(\partial_x^5\mathbf{m}^\pm\), which results that we can close the ansatz in \(H^8\times H^7\)-framework based on the Sobolev embedding theorem. Finally, by using the smallness of \(\epsilon\), one can close the ansatz (82 )-(83 ) and hence proves Theorem 1.

5 Appendix: Analytic tools↩︎

Some useful lemmas are given here. The first one is used to derive the pointwise estimates of Green’s function in the low frequency. Note that we have replaced \(t^{-(n+|\alpha|)}B_N(|x|,t)\) in Wang-Yang[46] with \((1+t)^{-(n+|\alpha|)}B_N(|x|,t)\) in the following lemma since it is in the low frequency part.

The next two lemmas are often used to deal with initial propagation and nonlinear coupling, respectively. We also state several typical cases for completeness.

Due to the non-conservation of the system (1 ), and the slower decay rate of the two fraction densities of the system (1 ), we have to construct some refined nonlinear convolution estimates, especially for the convolutions containing the Huygens waves or the Rieszwave-IV (from Riesz operator of negative first order) and the Rieszwave-III (from Riesz operator of zero order). In fact, Lin, Wang in [47], [48] investigated the Boltzmann equation and proved that \[\begin{align} &\int_0^t\!\int_{\mathbb{R}^3}\uwave{\mathbf{1}_{|x|\leq {\rm c}t}(1+t-\tau)^{-2}}\Big(1+\frac{|x-y|^2}{1+t-\tau}\Big)^{-\frac{3}{2}}(1+\tau)^{-4}\Big(1+\frac{(|y|-{\rm c}\tau)^2}{1+\tau}\Big)^{-2}dyd\tau\nonumber\\ \lesssim &\;(1+t)^{-2}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}+(1+t)^{-\frac{5}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-1}.\label{646140141} \end{align}\tag{93}\] Here the underlined is arising from the aforementioned Riesz wave-I and the conservative structure. Different from (93 ), we obtained a similar result in the following \(K_5\) without the factor \(\mathbf{1}_{|x|\leq {\rm c}t}\). Moreover, we also achieve a satisfactory pointwise space-time description of the Riesz wave-IV convolved with the Huygens wave in \(K_4\).

We mainly focus on the proof of (?? ) due to the slowest decay on both the temporal and spacial variables of Riesz wave-IV in Green’s function, and the other estimate could be derived similarly. To facilitate the estimates, we first decompose space-time domain into the following \(5\) regions: \[\begin{align} D_{1}& =\left \{ |x|\leq \sqrt{1+t}\right \} \,, \\ \;D_{2}& =\left \{ {\rm c}t-\sqrt{1+t}\leq |x|\leq {\rm c}t+\sqrt{1+t}\right \} \,, \\ D_{3}& =\left \{ |x|\geq {\rm c}t+\sqrt{1+t}\right \} \,, \\ D_{4}& =\Big \{ \sqrt{1+t}\leq |x|\leq \frac{1}{2}{\rm c}t\Big \} \,, \\ D_{5}& =\Big \{ \frac{1}{2}{\rm c}t\leq |x|\leq {\rm c}t-\sqrt{1+t}\Big \} \,. \end{align}\]

Proof of \(K_4\). (Riesz wave-IV convolved with Huygens wave)

Case 1: \(\left( x,t\right) \in D_{1}\). Direct computation gives\[\begin{align} K_4 &=&\left( \int_{0}^{\frac{t}{2}}+\int_{\frac{t}{2}}^{t}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-1}\int_{0}^{\frac{t}{2}}\int_{\mathbb{R}^{3}}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &&+\left( 1+t\right) ^{-4}\int_{\frac{t}{2}}^{t}\left( 1+t-\tau \right) ^{-1}(1+t-\tau)^{\frac{3}{4}}(1+\tau)^{\frac{5}{4}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-1}+\left( 1+t\right) ^{-2} \lesssim ( 1+t) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}, \end{align}\] where we have used Young inequality and the fact \[\label{6468} \int_{\mathbb{R}^{3}}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dy\lesssim (1+\tau)^{\frac{5}{2}}.\tag{94}\]

Case 2: \(\left( x,t\right) \in D_{2}\). We split the integral \(I\) into two parts\[K_4=\left( \int_{\frac{t}{4}}^{t}+\int_{0}^{\frac{t}{4}}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau =:I_{1}+I_{2}.\]

For \(I_{1}\), one can see that \[\begin{align} I_{1} &\lesssim &\left( 1+t\right) ^{-4}\int_{\frac{t}{4}}^{t}\left( 1+t-\tau \right) ^{-1}(1+t-\tau)^{\frac{3}{4}}(1+\tau)^{\frac{5}{4}}d\tau \\ &\lesssim & \left( 1+t\right) ^{-2}\lesssim \left( 1+t\right) ^{-2}\left( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\right) ^{-1}. \end{align}\]

For \(I_{2}\), we further decompose \(\mathbb{R}^{3}\) into two parts\[I_{2}=\left( \int_{0}^{\frac{t}{4}}\int_{\left \vert y\right \vert \leq \frac{2}{3}\left \vert x\right \vert }+\int_{0}^{\frac{t}{4}}\int_{\left \vert y\right \vert >\frac{2}{3}\left \vert x\right \vert }\right) \left( \cdots \right) dyd\tau =I_{11}+I_{12}.\]If \(\left \vert y\right \vert \leq \frac{2}{3}\left \vert x\right \vert\), then we have \[\left \vert x-y\right \vert \geq \left \vert x\right \vert -\left \vert y\right \vert \geq \frac{\left \vert x\right \vert }{3},\]and thus \[\begin{align} I_{11} &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}\int_{0}^{\frac{t}{4}}\int_{\mathbb{R}^{3}}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}\int_{0}^{\frac{t}{4}}\left( 1+\tau \right) ^{-4}\left( 1+\tau \right) ^{\frac{5}{2}}d\tau \lesssim \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]If \(\left \vert y\right \vert >\frac{2}{3}\left \vert x\right \vert\) and \(0\leq \tau \leq \frac{t}{4}\), \[\left \vert y\right \vert -{\rm c}\tau \geq \frac{2}{3}\left \vert x\right \vert -{\rm c}\tau \geq \frac{\left \vert x\right \vert }{3}+\frac{1}{3}\left( {\rm c}t-\sqrt{1+t}\right) -\frac{{\rm c}t}{4}\geq \frac{\left \vert x\right \vert }{3}+\frac{1}{24}{\rm c}t\]$ \(for\)t$. Note that \(\left( x,t\right) \in D_{2}\) with \(t\geq \frac{2+2\sqrt{1+{\rm c}^{2}}}{{\rm c}^{2}}\) implies that \(\left \vert x\right \vert \geq \sqrt{1+t}\). Hence for \(t\geq \frac{32}{{\rm c}^2}\), \[\begin{align} I_{12} &\lesssim &(1+t)^{-1}\left( 1+\left \vert x\right \vert ^{2}\right) ^{-1}\int_{0}^{\frac{t}{4}}\int_{\left \vert y\right \vert >\frac{2}{3}\left \vert x\right \vert }\left( 1+t-\tau \right) ^{-1}\bigg( 1+\frac{\left \vert x-y\right \vert ^{2}}{ 1+t-\tau }\bigg) ^{-1} \\ &&\cdot \left( 1+\tau \right) ^{-1}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-1}dyd\tau \\ &\lesssim &\left( 1+t \right) ^{-1}\left( 1+\left \vert x\right \vert ^{2}\right) ^{-1}\int_{0}^{\frac{t}{4}}(1+t-\tau)^{-1}(1+t-\tau)^{\frac{3}{4}}(1+\tau)^{-1}(1+\tau)^{-\frac{5}{4}}d\tau \\ &\lesssim &\left( 1+t \right) ^{-1}\left( 1+\left \vert x\right \vert ^{2}\right) ^{-1}(1+t)\\ &\lesssim& (1+t)^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}, \end{align}\]and for \(0\leq t\leq 10\), \[\begin{align} I_{12} &\lesssim &\left( 1+t\right) ^{-1}\int_{0}^{\frac{t}{4}}\int_{\mathbb{R}^{3}}\left( 1+\tau \right) ^{-4}\left( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\right) ^{-2}dyd\tau \\ &\lesssim &C\lesssim \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]Thereupon, \[K_4\lesssim \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}+\left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}.\]
Case 3: \(\left( x,t\right) \in D_{3}\). We split the integral \(I\) into two parts\[K_4=\left( \int_{0}^{\frac{t}{2}}+\int_{\frac{t}{2}}^{t}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau =:I_{1}+I_{2}.\]

For \(I_{1}\), we decompose \(\mathbb{R}^{3}\) into two parts as follows: \[I_{1}=\int_{0}^{\frac{t}{2}}\left( \int_{\left \vert y\right \vert \leq \frac{2}{3}\left \vert x\right \vert }+\int_{\left \vert y\right \vert >\frac{2}{3}\left \vert x\right \vert }\right) \left( \ldots \right) dyd\tau =:I_{11}+I_{12}.\]If \(\left \vert y\right \vert \leq \frac{2}{3}\left \vert x\right \vert\), then \[\left \vert x-y\right \vert \geq \frac{\left \vert x\right \vert }{3},\]and thus\[\begin{align} I_{11} &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}\int_{0}^{\frac{t}{2}}\int_{\mathbb{R}^{3}}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]If \(\left \vert y\right \vert >\frac{2}{3}\left \vert x\right \vert\) and \(0\leq \tau \leq \frac{t}{2}\), we have\[\left \vert y\right \vert -{\rm c}\tau \geq \frac{2}{3}\left \vert x\right \vert -\frac{{\rm c}t}{2}=\frac{1}{6}\left \vert x\right \vert +\frac{1}{2}\left( \left \vert x\right \vert -{\rm c}t\right) \geq \frac{1}{6}{\rm c}t,\]so that\[\begin{align} I_{12} &\lesssim &\big( 1+\left \vert x\right \vert ^{2}\big) ^{\!-1}\!\int_{0}^{\frac{t}{2}}\!\int_{\left \vert y\right \vert >\frac{2}{3}\left \vert x\right \vert }\left( 1\!+\!t\!-\!\tau \right) ^{-1}\left( 1+\frac{\left \vert x-y\right \vert ^{2}}{ 1\!+\!t\!-\!\tau }\right) ^{-1} \left( 1+\tau \right) ^{-3}\left( 1\!+\!\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\right) ^{\!-1}\!\!dyd\tau \\ &\lesssim &\big( 1+\left \vert x\right \vert ^{2}\big) ^{-1}\int_{0}^{\frac{t}{2}}(1+t-\tau)^{-1}(1+t-\tau)^{\frac{3}{4}}(1+\tau)^{-3}(1+\tau)^{\frac{5}{4}}d\tau \\ &\lesssim &(1+t)^{-\frac{1}{4}}\big( 1+\left \vert x\right \vert ^{2}\big) ^{-1}\lesssim (1+t)^{-\frac{5}{4}}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}, \end{align}\]since \(\left \vert x\right \vert \geq \sqrt{1+t}\). Therefore,\[I_{1}\lesssim \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}.\]

For \(I_{2}\), we decompose \(\mathbb{R}^{3}\) into two parts \[I_{2}=\int_{\frac{t}{2}}^{t}\left( \int_{\left \vert y\right \vert -\mathbf{c}\tau \leq \frac{\left \vert x\right \vert -ct}{2}}+\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -ct}{2}}\right) \left( \cdots \right) dyd\tau =I_{21}+I_{22}. \]If \(\left \vert y\right \vert -{\rm c}\tau \leq \frac{\left \vert x\right \vert -{\rm c}t}{2}\), we have \[\left \vert x-y\right \vert \geq \left \vert x\right \vert - \frac{\left \vert x\right \vert -{\rm c}t}{2} -{\rm c}\tau \geq \frac{\left \vert x\right \vert -{\rm c}t}{2},\]and thus\[\begin{align} I_{21} &\lesssim &\left( 1+\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}\right) ^{-1}\int_{\frac{t}{2}}^{t}\int_{\mathbb{R}^{3}}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}\right) ^{-1}\int_{\frac{t}{2}}^t(1+\tau)^{-4}(1+\tau)^{\frac{5}{2}}d\tau \\ &\lesssim &(1+t)^{-\frac{1}{2}}\big( 1+\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}\big) ^{-1}\lesssim (1+t)^{-\frac{3}{2}}\bigg(1+\frac{(|x|-ct)^2}{1+t}\bigg)^{-1}, \end{align}\]since \(|x|-ct>\sqrt{1+t}\).
If \(\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}\), then \[\begin{align} I_{22} &\lesssim &\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}\\ && \times\int_{\frac{t}{2}}^{t}\int_{\mathbb{R}^3}\left( 1+t-\tau \right) ^{-1}\bigg( 1+\frac{\left \vert x-y\right \vert ^{2}}{ 1+t-\tau }\bigg) ^{-1} \bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-1} dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}\int_{\frac{t}{2}}^{t}(1+t-\tau)^{-1} \left( 1+t-\tau \right)^{\frac{3}{4}}(1+\tau)^{\frac{5}{4}} d\tau \\ & \lesssim& \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]Therefore, \[I_{2}\lesssim \left( 1+t\right) ^{-3/2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}.\]Combining all the estimates, we get desired estimate\[K_4\lesssim \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}+\left( 1+t\right) ^{-3/2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}.\]
Case 4: \(\left( x,t\right) \in D_{4}\). In this region \(\sqrt{1+t}\leq \left \vert x\right \vert \leq \frac{{\rm c}t}{2}\). We split the integral \(K_4\) into three parts:\[\begin{align} K_4&=&\int_{0}^{t}\int_{\mathbb{R}^{3}}\left( 1+t-\tau \right) ^{-4}\left( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\right) ^{-2}\left( 1+\tau \right) ^{-1}\left( 1+\frac{\left \vert y\right \vert ^{2}}{ 1+\tau }\right) ^{-1}dyd\tau \\ &=&\left( \int_{0}^{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}+\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}+\int_{t-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau =:I_{1}+I_{2}+I_{3}. \end{align}\]

For \(I_{1}\), we decompose \(\mathbb{R}^{3}\) into two parts\[I_{1}=\int_{0}^{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}\left( \int_{\left \vert y\right \vert \leq \frac{{\rm c}\left( t-\tau \right) -\left \vert x\right \vert }{2}}+\int_{\left \vert y\right \vert >\frac{{\rm c}\left( t-\tau \right) -\left \vert x\right \vert }{2}}\right) =:I_{11}+I_{12}.\]If \(0\leq \tau \leq \frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}\) and \(\left \vert y\right \vert \leq \frac{{\rm c}\left( t-\tau \right) -\left \vert x\right \vert }{2}\), then\[{\rm c}\left( t-\tau \right) -\left \vert x-y\right \vert \geq {\rm c}\left( t-\tau \right) -\left \vert x\right \vert -\left \vert y\right \vert \geq \frac{{\rm c}t-\left \vert x\right \vert }{4},\]so that \[\begin{align} I_{11} &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}\\ &&\times\int_{0}^{\frac{t}{2}}\int_{\mathbb{R}^3}(1+t-\tau)^{-3}\bigg( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\bigg) ^{-1}\left( 1+\tau \right) ^{-1}\left( 1+\frac{\left \vert y\right \vert ^{2}}{ 1+\tau }\right) ^{-1}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}\int_0^{\frac{t}{2}}(1+t-\tau)^{-3}(1+t-\tau)^{-\frac{5}{4}}(1+\tau)^{-1}(1+\tau)^{\frac{3}{4}}d\tau\\ &\lesssim& \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}(1+t)^{-1}\lesssim \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg)^{-1}. \end{align}\]If \(0\leq \tau \leq \frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}\) and \(\left \vert y\right \vert >\frac{{\rm c}\left( t-\tau \right) -\left \vert x\right \vert }{2}\), then\[\left \vert y\right \vert >\frac{{\rm c}\left( t-\tau \right) -\left \vert x\right \vert }{2}\geq \frac{{\rm c}t-\left \vert x\right \vert }{4},\]and thus\[\begin{align} I_{12} &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}\\ &&\times\int_{0}^{\frac{t}{2}-\frac{|x|}{c}}\int_{\mathbb{R}^3}(1+t-\tau)^{-4}\bigg( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim& \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}(1+t)^{-\frac{1}{2}}\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg)^{-1}. \end{align}\]Therefore,\[I_{1}\lesssim \left( 1+t\right) ^{-\frac{5}{2}}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}.\]

For \(I_{2}\), we have to we decompose \(\mathbb{R}^{3}\) into two parts\[I_{2}=\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\left( \int_{\left \vert y\right \vert \leq {\rm c}\tau}+\int_{\left \vert y\right \vert >{\rm c}\tau}\right) =:I_{21}+I_{22}.\]We use the spherical coordinates to obtain \[\begin{align} I_{21} &\lesssim &\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\int_{0}^{c\tau}\int_{0}^{\pi }\left( 1+t-\tau \right) ^{-4}\bigg( 1+\frac{\big( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) \big) ^{2}}{1+t-\tau }\bigg) ^{-2} \\ &&\cdot \left( 1+\tau \right) ^{-1}\bigg( 1+\frac{\left \vert r\right \vert ^{2}}{ 1+\tau }\bigg) ^{-1}r^{2}\sin \theta d\theta drd\tau \\ &\lesssim &\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\int_{0}^{c\tau}\int_{\left \vert \left \vert x\right \vert -r\right \vert }^{\left \vert x\right \vert +r}\left( 1+t-\tau \right) ^{-4}\bigg( 1+\frac{\left( z-{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\bigg) ^{-2}z( 1+\tau ) ^{-1} \\ &&\cdot \left( 1+\frac{r^{2}}{ 1+\tau }\right) ^{-1}\frac{r}{\left \vert x\right \vert }dzdrd\tau \\ &\lesssim &\left \vert x\right \vert ^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\int_{0}^{c\tau}\left( 1+t-\tau \right) ^{-4+\frac{3}{2}}\left( 1+\tau \right) ^{-1}\cdot r\left( 1+\frac{r^{2}}{ 1+\tau }\right) ^{-1}drd\tau \\ &\lesssim &\left \vert x\right \vert ^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\left( 1+t-\tau \right) ^{-\frac{5}{2}}\ln( 1+{\rm c}\tau)d\tau \\ &\lesssim &\left \vert x\right \vert ^{-1}\big(1+\frac{|x|}{2c}\big)^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\left( 1+t-\tau \right) ^{-\frac{3}{2}}\ln( 1+{\rm c}\tau)d\tau \\ &\lesssim &\left( 1+\left \vert x\right \vert \right) ^{-2}\lesssim \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}, \end{align}\]where we used \(z=\sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }\) and \(\sin \theta d\theta =\frac{z}{r\left \vert x\right \vert }dz\).

For \(I_{22}\), since \(|y|>{\rm c}\tau\), one has \[\begin{align} |y|>{\rm c}\tau>{\rm c}\big(\frac{t}{2}-\frac{|x|}{2{\rm c}}\big)=\frac{{\rm c}t-|x|}{2}\geq\frac{|x|}{2}, \end{align}\] which gives that \[\begin{align} I_{22}&\lesssim& (1+t)^{-1}\big(1+\frac{|x|^2}{1+t}\big)^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}(1+t-\tau)^{-4}\left( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\right) ^{-2}dyd\tau\nonumber\\ &\lesssim& (1+t)^{-1}\big(1+\frac{|x|^2}{1+t}\big)^{-1}. \end{align}\]

For \(I_{3}\), we decompose \(\mathbb{R}^{3}\) into two parts\[I_{3}=\int_{t-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t}\left( \int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}}+\int_{\left \vert y\right \vert >\frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}}\right) =:I_{31}+I_{32}.\]If \(t-\frac{\left \vert x\right \vert }{2{\rm c}}\leq \tau \leq t\), \(\left \vert y\right \vert \leq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\), then \[\left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \geq \left \vert x\right \vert -\left \vert y\right \vert -{\rm c}\left( t-\tau \right) \geq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\geq \frac{\left \vert x\right \vert }{4}.\]If \(t-\frac{\left \vert x\right \vert }{2{\rm c}}\leq \tau \leq t\), \(\left \vert y\right \vert >\frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\), then\[\left \vert y\right \vert >\frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\geq \frac{\left \vert x\right \vert }{4}.\]Hence, since \(\frac{{\rm c}t}{2}\geq\left \vert x\right \vert \geq \sqrt{1+t}\), it holds that \[\begin{align} I_{31} &\lesssim &(1+t)^{-1}\left( 1+\frac{|x|^2}{1+t}\right) ^{-1}\\ &&\cdot\int_{t-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t}\int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert \!-\!{\rm c}\left( t\!-\!\tau \right) }{2}}\left( 1\!+\!t\!-\!\tau \right) ^{-3}\left( 1+\frac{\left( \left \vert x\!-\!y\right \vert \!-\!{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\right) ^{\!-1}\!\!\!(1\!+\!\tau)^{-1}\big(1+\frac{|y|^2}{1+\tau}\big)^{-1}\!dyd\tau \\ &\lesssim &(1+t)^{-1}\left( 1+\frac{|x|^2}{1+t}\right) ^{-1}\int_{t-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t}(1+t-\tau)^{-\frac{7}{4}}(1+\tau)^{-\frac{1}{4}}d\tau\\ &\lesssim& \left( 1+t\right) ^{-\frac{5}{4}}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}, \end{align}\] and\[\begin{align} I_{32} &\lesssim &\left( 1+t\right) ^{-1}\left( 1+\frac{\left \vert x\right \vert ^{2}}{ 1+t }\right) ^{-1}\int_{t-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t}\int_{\mathbb{R}^{3}}\left( 1+t-\tau \right) ^{-4}\left( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\right) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]Combining all the estimates yields \[K_4\lesssim \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t }\bigg) ^{-\frac{3}{2}}+\left( 1+t\right) ^{-5/2}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}.\]
Case 5: \(\left( x,t\right) \in D_{5}\). In this region \(\frac{{\rm c}t}{2}\leq \left \vert x\right \vert \leq {\rm c}t-\sqrt{1+t}\). We split \(K_4\) into three parts\[\begin{align} K_4&=&\int_{0}^{t}\int_{\mathbb{R}^{3}}\left( 1+t-\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\bigg) ^{-2}\left( 1+\tau \right) ^{-2 }\bigg( 1+\frac{\left \vert y\right \vert ^{2}}{ 1+\tau}\bigg) ^{-1}dyd\tau \\ &=&\left( \int_{0}^{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}+\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }+\int_{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }^{t}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau =:I_{1}+I_{2}+I_{3}. \end{align}\]

For \(I_{1}\), the estimate is the same as the \(I_{1}\) of Case 4, so \[I_{1}\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}.\]

For \(I_{2}\), we have to we decompose \(\mathbb{R}^{3}\) into two parts\[I_{2}=\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }\left( \int_{\left \vert y\right \vert \leq {\rm c}(t-\tau)}+\int_{\left \vert y\right \vert >{\rm c}(t-\tau)}\right) =:I_{21}+I_{22}.\]In terms of \(I_{21}\), similar to the \(I_{21}\) of Case 4, we use the spherical coordinates to obtain\[\begin{align} I_{21} &\lesssim &\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }\int_{0}^{{\rm c}(t-\tau)}\int_{0}^{\pi }\left( 1+t-\tau \right) ^{-4}\Bigg( 1+\frac{\big( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) \big)^{2}}{1+t-\tau }\Bigg)^{-2} \\ &&\cdot \left( 1+\tau \right) ^{-1}\bigg( 1+\frac{\left \vert r\right \vert ^{2}}{ 1+\tau }\bigg) ^{-1}r^{2}\sin \theta d\theta drd\tau \\ &\lesssim &\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }\int_{0}^{{\rm c}(t-\tau)}\int_{\left \vert \left \vert x\right \vert -r\right \vert }^{\left \vert x\right \vert +r}\left( 1+t-\tau \right) ^{-4}\left( 1+\frac{\left( z-{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\right) ^{-2}z\left( 1+\tau \right) ^{-1} \\ &&\cdot \bigg( 1+\frac{r^{2}}{ 1+\tau }\bigg) ^{-1}\frac{r}{\left \vert x\right \vert }dzdrd\tau \\ &\lesssim &\left \vert x\right \vert ^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{t-\frac{\left \vert x\right \vert }{2{\rm c}}}\int_{0}^{{\rm c}(t-\tau)}\left( 1+t-\tau \right) ^{-4+\frac{3}{2}}\left( 1+\tau \right) ^{-1}\cdot r\left( 1+\frac{r^{2}}{ 1+\tau }\right) ^{-1}drd\tau \\ &\lesssim &\left \vert x\right \vert ^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }\left( 1+t-\tau \right) ^{-\frac{5}{2}}\ln( 1+{\rm c}(t-\tau))d\tau \\ &\lesssim &\left \vert x\right \vert ^{-1}\big(1+\frac{|x|}{2{\rm c}}\big)^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }\left( 1+t-\tau \right) ^{-\frac{3}{2}}\ln( 1+{\rm c}(t-\tau))d\tau \\ &\lesssim &\left( 1+\left \vert x\right \vert \right) ^{-2}\lesssim \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]For \(I_{22}\), since \(|y|>{\rm c}(t-\tau)\), one has \[\begin{align} |y|>{\rm c}(t-\tau)\geq {\rm c}\bigg(\frac{t}{2}-\frac{3}{4}\bigg(t-\frac{|x|}{{\rm c}}\bigg)\bigg)\geq {\rm c}\bigg(\frac{t}{2}-\frac{3}{4}\bigg(t-\frac{{\rm c}t/2}{2}\bigg)\bigg)=\frac{{\rm c}t}{8}>\frac{|x|}{8}, \end{align}\] which gives that \[\begin{align} I_{22}&\lesssim& (1+t)^{-1}\bigg(1+\frac{|x|^2}{1+t}\bigg)^{-1}\int_{\frac{t}{2}-\frac{\left \vert x\right \vert }{2{\rm c}}}^{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }(1+t-\tau)^{-4}\bigg( 1+\frac{\left( \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\bigg) ^{-2}dyd\tau\nonumber\\ &\lesssim& (1+t)^{-1}\bigg(1+\frac{|x|^2}{1+t}\bigg)^{-1}. \end{align}\]

For \(I_{3}\), we decompose \(\mathbb{R}^{3}\) into two parts\[I_{3}=\int_{\frac{1}{2}\big( t+\frac{3}{2}\big( t-\frac{\left \vert x\right \vert }{{\rm c}}\big) \big) }^{t}\left( \int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}}+\int_{\left \vert y\right \vert >\frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}}\right) \left( \cdots \right) dyd\tau =:I_{31}+I_{32}.\]If \(\frac{1}{2}\big( t+\frac{3}{2}\big( t-\frac{\left \vert x\right \vert }{{\rm c}}\big) \big) \leq \tau \leq t\) and \(\left \vert y\right \vert \leq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\), then \[\left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \geq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\geq \frac{{\rm c}t+\left \vert x\right \vert }{8}= \frac{{\rm c}t-\left \vert x\right \vert }{8}+\frac{|x|}{4}.\]If \(\frac{1}{2}\big( t+\frac{3}{2}\big( t-\frac{\left \vert x\right \vert }{{\rm c}}\big) \big) \leq \tau \leq t\) and \(\left \vert y\right \vert >\frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\), then \[\left \vert y\right \vert >\frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}\geq \frac{{\rm c}t-\left \vert x\right \vert }{8}+\frac{|x|}{4} .\]Hence, \[\begin{align} I_{31} \!\!\!\!\!\!\!\!&&\lesssim(1+t)^{-1}\left( 1+\frac{|x|^2}{1+t}\right) ^{-1}\int_{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }^{t}\int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert -{\rm c}\left( t-\tau \right) }{2}}\left( 1+t-\tau \right) ^{-3}\\ &&\cdot\bigg( 1+\frac{\left( \left \vert x\!-\!y\right \vert \!-\!{\rm c}\left( t-\tau \right) \right) ^{2}}{1\!+\!t\!-\!\tau }\bigg) ^{\!\!-1}\!\!\!\!(1\!+\!\tau)^{-1}\big(1+\frac{|y|^2}{1+\tau}\big)^{-1}\!dyd\tau \\ &&\lesssim(1+t)^{-1}\bigg( 1+\frac{|x|^2}{1+t}\bigg) ^{-1}\int_{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }^{t}(1+t-\tau)^{-\frac{7}{4}}(1+\tau)^{-\frac{1}{4}}d\tau\\ &&\lesssim \left( 1+t\right) ^{-\frac{5}{4}}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]Similarly, \(I_{32}\) can be estimated as follows: \[\begin{align} I_{32} &\lesssim &(1+t)^{-1}\bigg( 1+\frac{|x|^2}{1+t}\bigg) ^{-1}\\ &&\cdot\int_{\frac{1}{2}\left( t+\frac{3}{2}\left( t-\frac{\left \vert x\right \vert }{{\rm c}}\right) \right) }^{t}\int_{\left \vert y\right \vert \geq \frac{\left \vert x\right \vert \!-\!{\rm c}\left( t\!-\!\tau \right) }{2}}\left( 1\!+\!t\!-\!\tau \right) ^{-4}\left( 1+\frac{\left( \left \vert x\!-\!y\right \vert \!-\!{\rm c}\left( t-\tau \right) \right) ^{2}}{1+t-\tau }\right) ^{\!-2}dyd\tau \\ &\lesssim &(1+t)^{-1}\bigg( 1+\frac{|x|^2}{1+t}\bigg) ^{-1}. \end{align}\]By summing up the estimates in five domains \(D_i\) with \(i=1,2,3,4,5\), one has \[K_4\leq \left( 1+t\right) ^{-1}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{ 1+t }\bigg) ^{-\frac{3}{2}}+\left( 1+t\right) ^{-1}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1},\] and has completed the proof of (?? ) about Riesz wave-\(I\) convolved with Huygens wave.

Finally, we shall establish two nonlinear convolution estimates on the Huygens wave convolved with nonlinear Huygens wave or nonlinear diffusion wave for the non-conservation of the system (1 ), which is completely different from the previous results due to the slower decay rate of Huygens wave in Green’s function: \((1+t)^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}\) compared with \(K_3\) in Lemma 10.

We only focus on the proof of \(K_7\) since the proof of \(K_6\) is relatively easier.

Proof of \(K_7\). (Huygens wave convolved with Huygens wave for non-conservative system)

Case 1: \(\left( x,t\right) \in D_{1}\cup D_{2}\). Direct computation gives \[\begin{align} K_7 &\lesssim &\int_{0}^{\frac{t}{2}}\int_{\mathbb{R}^{3}}( 1+t) ^{-2}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &&+\int_{\frac{t}{2}}^{t}\int_{\mathbb{R}^{3}}\left( 1+t-\tau \right) ^{-2}\bigg(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\bigg)^{-N}\left( 1+t\right) ^{-4}dyd\tau \\ &\lesssim &\left( 1+t\right)^{-2}\int_{0}^{\frac{t}{2}}\left( 1+\tau \right)^{-4}\left( 1+\tau \right) ^{\frac{5}{2}}d\tau +\left( 1+t\right)^{-4}\int_{\frac{t}{2}}^{t}\left( 1+t-\tau \right) ^{-2}\left( 1+t-\tau \right) ^{\frac{5}{2}}d\tau \\ &\lesssim &(1+t)^{-2}, \end{align}\] which implies that \[\begin{align} K_7 \lesssim ( 1+t)^{-2}\bigg(1+\frac{|x|^2}{1+t}\bigg)^{-N},\;\;\;\;\;\;{\rm in}\;D_1,\\[2mm] K_7 \lesssim ( 1+t)^{-2}\bigg(1+\frac{(|x|-{\rm c}t)^2}{1+t}\bigg)^{-N},\;\;\;\;\;{\rm in}\;D_2. \end{align}\]

Case 2: \(\left( x,t\right) \in D_{3}\). We split the integral \(K_7\) into four parts\[\begin{align} J \!&\!=\!&\!\int_{0}^{\frac{t}{2}}\bigg( \int_{\left \vert y\right \vert -{\rm c}\tau \leq \frac{\left \vert x\right \vert -{\rm c}t}{2}}+\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\bigg) \left( \cdots \right) dyd\tau\\ && +\int_{\frac{t}{2}}^{t}\bigg( \int_{\left \vert y\right \vert -{\rm c}\tau \leq \frac{\left \vert x\right \vert -{\rm c}t}{2}}+\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\bigg) \left( \cdots \right) dyd\tau\\[2mm] \!&\!=\!&\!J_{11}+J_{12}+J_{21}+J_{22}. \end{align}\]

Note that \(\left \vert y\right \vert -{\rm c}\tau \leq \frac{\left \vert x\right \vert -{\rm c}t}{2}\) implies \(\left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) \geq \left \vert x\right \vert -\left \vert y\right \vert -{\rm c}\left( t-\tau \right) \geq \frac{\left \vert x\right \vert -{\rm c}t}{2}\). Hence, \[\begin{align} J_{11} &\lesssim &\left( 1+t\right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}\int_{0}^{\frac{t}{2}}\int_{\left \vert y\right \vert -{\rm c}\tau \leq \frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}\int_{0}^{\frac{t}{2}}\left( 1+\tau \right) ^{-4}\left( 1+\tau \right) ^{\frac{5}{2}}d\tau\\ &\lesssim& \left( 1+t\right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}. \end{align}\]

For \(J_{12}\), by using Young inequality, we have \[\begin{align} J_{12} &\lesssim &\int_{0}^{\frac{t}{2}}\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+t-\tau \right)^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\\ &&\;\;\;\;\;\;\;\cdot \left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+\tau }\bigg) ^{-1}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-1}dyd\tau \\ &\lesssim &(1+t)^{-2}\Big( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\Big) ^{-1}\int_{0}^{\frac{t}{2}}\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-1}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\\ &&\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\cdot \left( 1+\tau \right) ^{-3}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right)^{2}}{1+\tau }\bigg) ^{-1}dyd\tau \\ &\lesssim &(1+t)^{-2}\Big( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\Big) ^{-1}\int_{0}^{\frac{t}{2}}(1+t-\tau)^{-1}(1+t-\tau)(1+\tau)^{-3}(1+\tau)^{\frac{3}{2}}d\tau \\ &\lesssim &(1+t)^{-2}\Big( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\Big) ^{-1}. \end{align}\] For \(J_{21}\) and \(J_{22}\), it follows that \[\begin{align} J_{21} &\lesssim &\int_{\frac{t}{2}}^{t}\int_{\left \vert y\right \vert -{\rm c}\tau \leq \frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t-\tau}\Big)^{-N/2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N/2}\\ && \;\;\;\;\;\;\;\;\cdot\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N/2}\int_{\frac{t}{2}}^{t}\left( 1+t-\tau \right) ^{-2}\left( 1+t-\tau \right) ^{\frac{5}{2}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{5}{2}}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N/2}, \end{align}\]\[\begin{align} J_{22} &\lesssim &\int_{\frac{t}{2}}^{t}\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\\ &&\cdot\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\int_{\frac{t}{2}}^{t}\left( 1+t-\tau \right) ^{-2}\left( 1+t-\tau \right) ^{\frac{5}{2}}\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-2}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{5}{2}}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-2}. \end{align}\]

Gathering all the estimates, we can conclude that \[K_7\lesssim \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t\right) ^{2}}{1+t}\bigg) ^{-1}.\]
Case 3: \(\left( x,t\right) \in D_{4}\). We split the integral into four parts\[\begin{align} K_7&=&\bigg( \int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4\mathbf{c}}}+\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}+\int_{\frac{t}{2}}^{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }+\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\bigg) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau \\ &=:&J_{1}+J_{2}+J_{3}+J_{4}. \end{align}\]

For \(J_{1}\), we decompose \(\mathbb{R}^{3}\) into two parts\[J_{1}=\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\left( \int_{\left \vert y\right \vert \leq \frac{{\rm c}t-\left \vert x\right \vert }{2}}+\int_{\left \vert y\right \vert >\frac{{\rm c}t-\left \vert x\right \vert }{2}}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau =:J_{11}+J_{12}.\]If \(0\leq \tau \leq \frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}\) and \(\left \vert y\right \vert \leq \frac{{\rm c}t-\left \vert x\right \vert }{2}\), then \[{\rm c}\left( t-\tau \right) -\left \vert x-y\right \vert \geq {\rm c}\left( t-\tau \right) -\left \vert x\right \vert -\left \vert y\right \vert \geq \frac{{\rm c}t-\left \vert x\right \vert }{4}\geq \frac{{\rm c}t}{8}.\]Hence, it holds that \[\begin{align} J_{11} &\lesssim &\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\!\!\int_{\left \vert y\right \vert \leq \frac{{\rm c}t-\left \vert x\right \vert }{2}}\left( 1\!+\!t\!-\!\tau\right) ^{-2}\Big(1+\frac{(|x\!-\!y|\!-\!{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{\!-N}\left( 1\!+\!\tau \right) ^{-4}\bigg( 1\!+\!\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{\!-2}\!dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-2}\Big(1\!+\!\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}\!\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\!\left( 1+\tau \right) ^{-4}\left( 1+\tau \right) ^{\frac{5}{2}}d\tau \lesssim \left( 1+t\right) ^{-2}\Big(1\!+\!\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}, \end{align}\]and\[\begin{align} J_{12} &\lesssim &\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\!\!\!\int_{\left \vert y\right \vert >\frac{{\rm c}t-\left \vert x\right \vert }{2}}\left( 1\!+\!t\!-\!\tau \right) ^{-2}\Big(1\!+\!\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{\!-N}\!\left( 1\!+\!\tau \right) ^{-2}\left( 1\!+\!\left( {\rm c}t\!-\!\left \vert x\right \vert \right)^{2}\right) ^{\!-2}\!dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\left( 1+t-\tau \right) ^{-2}\left( 1+t-\tau \right) ^{\frac{5}{2}}\left( 1+\tau \right) ^{-2}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{7}{2}}\lesssim \left( 1+t\right) ^{-\frac{1}{2}}\left( 1+\left \vert x\right \vert \right) ^{-3}\lesssim \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-\frac{3}{2}}, \end{align}\]since \(\sqrt{1+t}\leq \left \vert x\right \vert \leq \frac{{\rm c}t}{2}\). Thus, \[\begin{align} J_{1}\lesssim \left( 1+t\right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N}+\left( 1+t\right) ^{-2}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-\frac{3}{2}}. \end{align}\]

For \(J_{2}\), \(\tau\geq\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}\) implies that \[\begin{align} J_{2} &\lesssim &\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}\int_{\mathbb{R}^{3}}\left( 1+t\right) ^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-2}\left( 1+\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}\right) ^{-4}\int_{0}^{t}\int_{\mathbb{R}^{3}}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-2}\left( 1+t\right) ^{-4}\left( 1+t\right) ^{3}\lesssim \left( 1+t\right) ^{-3}\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-\frac{3}{2}}, \end{align}\]and for \(J_3\), the relation \(t-\tau\geq t-\frac{t}{2}-\frac{1}{4}(t+\frac{|x|}{{\rm c}})=\frac{{\rm c}t-|x|}{4c}>\frac{|x|}{2{\rm c}}\) gives that \[\begin{align} J_{3} &\lesssim &\int_{\frac{t}{2}}^{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }\!\!\!\int_{\mathbb{R}^{3}}\left( 1\!+\!t\!-\!\tau \right) ^{-2}\Big(1\!+\!\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{\!-N}\!\left( 1\!+\!t\right) ^{-4}\bigg( 1\!+\!\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}\!dyd\tau \\ &\lesssim &\Big( 1+\frac{|x|}{2{\rm c}}\Big)^{-2}\left( 1+t\right) ^{-4}\int_{0}^{t}\int_{\mathbb{R}^{3}}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\Big( 1+\frac{|x|}{2c}\Big)^{-2}\left( 1+t\right) ^{-1}\lesssim (1+t)^{-\frac{3}{2}}\Big(1+\frac{|x|^2}{1+t}\Big)^{-\frac{3}{2}}. \end{align}\]

As for \(J_{4}\), we decompose \(\mathbb{R}^{3}\) into two parts\[J_{4}=\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\left( \int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert +{\rm c}t}{2}}+\int_{\left \vert y\right \vert >\frac{\left \vert x\right \vert +{\rm c}t}{2}}\right) \left( \cdots \right) dyd\tau =:J_{41}+J_{42}.\]If \(\frac{t}{2}+\frac{1}{4}\big( t+\frac{\left \vert x\right \vert }{{\rm c}}\big) \leq \tau \leq t\) and \(\left \vert y\right \vert \leq \frac{\left \vert x\right \vert +{\rm c}t}{2}\), then\[{\rm c}\tau -\left \vert y\right \vert \geq \frac{{\rm c}t}{2}+\frac{1}{4}\left( {\rm c}t+\left \vert x\right \vert \right) -\frac{\left \vert x\right \vert +{\rm c}t}{2}\geq \frac{{\rm c}t-\left \vert x\right \vert }{4}\geq \frac{{\rm c}t}{8}.\]If \(\frac{t}{2}+\frac{1}{4}\big( t+\frac{\left \vert x\right \vert }{{\rm c}}\big) \leq \tau \leq t\) and \(\left \vert y\right \vert >\frac{\left \vert x\right \vert +{\rm c}t}{2}\), then\[\begin{align} \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) &\geq &\left \vert y\right \vert -\left \vert x\right \vert -{\rm c}\left( t-\tau \right) \geq \frac{\left \vert x\right \vert +{\rm c}t}{2}-\left \vert x\right \vert -{\rm c}t+\frac{{\rm c}t}{2}+\frac{1}{4}\left( {\rm c}t+\left \vert x\right \vert \right) \\ &\geq &\frac{{\rm c}t-\left \vert x\right \vert }{4}\geq \frac{{\rm c}t}{8}. \end{align}\]Hence, \[\begin{align} J_{41} &\lesssim &\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert +{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\\ &&\cdot\left( 1+t\right) ^{-4}\bigg( 1+\frac{( {\rm c}t-\left \vert x\right \vert )^{2}}{1+t}\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-2}(1+t)^{\frac{3}{2}} \\ &\lesssim &\left( 1+t\right) ^{-\frac{9}{2}}\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\left( 1+\left \vert x\right \vert \right) ^{-3}\lesssim \left( 1+t\right) ^{-3}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-\frac{3}{2}}, \end{align}\]\[\begin{align} J_{42} &\lesssim &\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\int_{\left \vert y\right \vert >\frac{\left \vert x\right \vert +{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N/2}\\ && \cdot \Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N/2}\left( 1+t\right) ^{-4}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N/2}\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\left( 1+t-\tau \right) ^{-\frac{5}{2}+\frac{5}{2}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-3}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N/2}. \end{align}\] As a result, \[\begin{align} K_7\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-\frac{3}{2}}+\left( 1+t\right) ^{-2}\Big(1+\frac{(|x|-{\rm c}t)^2}{1+t}\Big)^{-N/2} \end{align}\] for \(\left( x,t\right) \in D_{4}\).

Case 4: \(\left( x,t\right) \in D_{5}\). We split the integral \(K_7\) into five parts\[K_7=\left( \int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}+\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}+\int_{\frac{t}{2}}^{\frac{|x|}{2{\rm c}}+\frac{1}{4}(t+\frac{|x|}{{\rm c}})} +\int_{\frac{|x|}{2{\rm c}}+\frac{1}{4}(t+\frac{|x|}{{\rm c}})}^{\frac{t}{2}+\frac{1}{4}(t+\frac{|x|}{{\rm c}})} +\int^{t}_{\frac{t}{2}+\frac{1}{4}(t+\frac{|x|}{{\rm c}})}\right) \int_{\mathbb{R}^{3}}\left( \cdots \right) dyd\tau =:\sum_{i=1}^{5}J_{i}.\]

For \(J_{1}\), we decompose \(\mathbb{R}^{3}\) into two parts\[J_{1}=\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\left( \int_{\left \vert y\right \vert \leq \frac{{\rm c}t-\left \vert x\right \vert }{2}}+\int_{\left \vert y\right \vert >\frac{{\rm c}t-\left \vert x\right \vert }{2}}\right) \left( \cdots \right) dyd\tau =:J_{11}+J_{12}.\]If \(0\leq \tau \leq \frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}\) and \(\left \vert y\right \vert \leq \frac{{\rm c}t-\left \vert x\right \vert }{2}\), then \[{\rm c}\left( t-\tau \right) -\left \vert x-y\right \vert \geq {\rm c}\left( t-\tau \right) -\left( \left \vert x\right \vert +\left \vert y\right \vert \right) \geq \frac{{\rm c}t-\left \vert x\right \vert }{4}.\]If \(0\leq \tau \leq \frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}\) and \(\left \vert y\right \vert >\frac{{\rm c}t-\left \vert x\right \vert }{2}\), then \[\left \vert y\right \vert -{\rm c}\tau \geq \frac{{\rm c}t-\left \vert x\right \vert }{4}.\] Hence,\[\begin{align} J_{11} &\lesssim &\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\int_{\left \vert y\right \vert \leq \frac{{\rm c}t-\left \vert x\right \vert }{2}}\left( 1+t\right) ^{-2}\Big(1+\frac{({\rm c}t-|x|)^2}{1+t}\Big)^{-N}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-2}\Big(1+\frac{({\rm c}t-|x|)^2}{1+t}\Big)^{-N}\int_{0}^{\frac{{\rm c}t-\left \vert x\right \vert }{4{\rm c}}}\left( 1+\tau \right) ^{-4+\frac{5}{2}}d\tau \lesssim \left( 1+t\right) ^{-2}\Big(1+\frac{({\rm c}t-|x|)^2}{1+t}\Big)^{-N}. \end{align}\]For \(J_{12}\), by using Young inequality, we have \[\begin{align} J_{12} &\lesssim &\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\\ &&\cdot \left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+\tau }\bigg) ^{-1}\bigg( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-1}dyd\tau \\ &\lesssim &(1+t)^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\bigg) ^{-1}\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}\int_{\left \vert y\right \vert -{\rm c}\tau >\frac{\left \vert x\right \vert -{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-1}\\ &&\cdot \Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\left( 1+\tau \right) ^{-3}\left( 1+\frac{\left( \left \vert y\right \vert -{\rm c}\tau \right) ^{2}}{1+\tau }\right) ^{-1}dyd\tau \\ &\lesssim &(1+t)^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\bigg) ^{-1}\int_{0}^{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}(1+t-\tau)^{-1}(1+t-\tau)(1+\tau)^{-3}(1+\tau)^{\frac{3}{2}}d\tau \\ &\lesssim &(1+t)^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\bigg) ^{-1}. \end{align}\]Therefore, we get\[J_{1}\lesssim (1+t)^{-2}\bigg( 1+\frac{\left( \left \vert x\right \vert -{\rm c}t \right) ^{2}}{1+t}\bigg) ^{-1}.\]

For \(J_{2}\), we use the spherical coordinates to obtain

\[\begin{align} J_{2} &\lesssim &\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}\int_{0}^{\infty }\int_{0}^{\pi }\left( 1+t-\tau \right) ^{-2}\bigg(1+\frac{( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) ) ^{2}}{ 1+t-\tau}\bigg)^{-N}\\ &&\cdot\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}r^{2}\sin \theta d\theta drd\tau \\ &\lesssim &\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}\int_{0}^{\infty }\int_{\left \vert \left \vert x\right \vert -r\right \vert }^{\left \vert x\right \vert +r}\left( 1+t-\tau \right) ^{-2}\bigg(1+\frac{( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) ) ^{2}}{ 1+t-\tau}\bigg)^{-N}\\ &&\cdot\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}rz\frac{1}{\left \vert x\right \vert }dzdrd\tau \\ &\lesssim &\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}\int_{0}^{\infty }\int_{0}^{\infty }\left( 1+t-\tau \right) ^{-2}\bigg(1+\frac{( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) ) ^{2}}{ 1+t-\tau}\bigg)^{-N}\\ &&\cdot\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}rz\frac{1}{\left \vert x\right \vert }dzdrd\tau \\ &\lesssim &\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}\int_{0}^{\infty }\left( 1+t-\tau \right) ^{-2+\frac{3}{2}}\left( 1+\tau \right) ^{-4}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}\frac{r}{\left \vert x\right \vert }drd\tau \\ &\lesssim &\left( 1+t\right) ^{-2+\frac{3}{2}}\left \vert x\right \vert ^{-1}\int_{\frac{t}{4}-\frac{\left \vert x\right \vert }{4{\rm c}}}^{\frac{t}{2}}\left( 1+\tau \right) ^{-4+\frac{3}{2}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{3}{2}}\left( 1+{\rm c}t-\left \vert x\right \vert \right) ^{-\frac{3}{2}}\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\left( 1+{\rm c}t-\left \vert x\right \vert \right) ^{\frac{1}{2}}\left( 1+{\rm c}t-\left \vert x\right \vert \right) ^{-2} \\ &\lesssim &\left( 1+t\right)^{-2}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}, \end{align}\]where we have used the fact \({\rm c}t-|x|\geq\sqrt{1+t}\) in \(D_5\).

For \(J_{3}\), we use the spherical coordinates again to obtain\[\begin{align} J_{3} &\lesssim &\int_{\frac{t}{2}}^{\frac{3|x|}{4{\rm c}}+\frac{t}{4}}\int_{0}^{\infty }\int_{0}^{\pi }\left( 1+t-\tau \right) ^{-2}\bigg(1+\frac{( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) ) ^{2}}{ 1+t-\tau}\bigg)^{-N}\\ &&\cdot\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}r^{2}\sin \theta d\theta drd\tau \\ &\lesssim &\int_{\frac{t}{2}}^{\frac{3|x|}{4{\rm c}}+\frac{t}{4}}\int_{0}^{\infty }\int_{\left \vert \left \vert x\right \vert -r\right \vert }^{\left \vert x\right \vert +r}\left( 1+t-\tau \right) ^{-2}\bigg(1+\frac{( \sqrt{\left \vert x\right \vert ^{2}+r^{2}-2r\left \vert x\right \vert \cos \theta }-{\rm c}\left( t-\tau \right) ) ^{2}}{ 1+t-\tau}\bigg)^{-N}\\ &&\cdot\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}rz\frac{1}{\left \vert x\right \vert }dzdrd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\left \vert x\right \vert ^{-1}\int_{\frac{t}{2}}^{\frac{3|x|}{4{\rm c}}+\frac{t}{4}} \int_{0}^{\infty }\left( 1+t-\tau \right) ^{-2+\frac{3}{2}}\bigg( 1+\frac{\left( r-{\rm c}\tau \right) ^{2}}{1+\tau }\bigg) ^{-2}rdrd\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{7}{2}}\int_{\frac{t}{2}}^{\frac{3|x|}{4{\rm c}}+\frac{t}{4}}\left( 1+t-\tau \right) ^{-\frac{1}{2}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{7}{2}}\bigg[\bigg(1+\frac{t}{2}\bigg)^{\frac{1}{2}}- \bigg(1+\frac{{\rm c}t-|x|}{4{\rm c}}\bigg)^{\frac{1}{2}}\bigg] \\ &\lesssim& \left( 1+t\right) ^{-3}\lesssim\left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}. \end{align}\] For \(J_{4}\), similar to \(J_{3}\), we have \[\begin{align} J_{4} &\lesssim &\left( 1+t\right) ^{-\frac{7}{2}}\int_{\frac{3|x|}{4{\rm c}}+\frac{t}{4}}^{\frac{|x|}{4{\rm c}}+\frac{3t}{4}}\left( 1+t-\tau \right) ^{-\frac{1}{2}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{7}{2}}\bigg[\left(1+\frac{3({\rm c}t-|x|)}{4{\rm c}}\right)^{\frac{1}{2}} -\left(1+\frac{({\rm c}t-|x|)}{2{\rm c}}\right)^{\frac{1}{2}}\bigg]\lesssim (1+t)^{-3}\\ &\lesssim& \left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}. \end{align}\] For \(J_{5}\), we decompose \(\mathbb{R}^{3}\) into two parts\[J_{5}=\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\left( \int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert +{\rm c}t}{2}}+\int_{\left \vert y\right \vert >\frac{\left \vert x\right \vert +{\rm c}t}{2}}\right) \left( \cdots \right) dyd\tau =:J_{51}+J_{52}.\] If \(\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) \leq \tau \leq t\) and \(\left \vert y\right \vert \leq \frac{\left \vert x\right \vert +{\rm c}t}{2}\), then\[{\rm c}\tau -\left \vert y\right \vert \geq \frac{{\rm c}t}{2}+\frac{1}{4}\left( {\rm c}t+\left \vert x\right \vert \right) -\frac{\left \vert x\right \vert +{\rm c}t}{2}\geq \frac{{\rm c}t-\left \vert x\right \vert }{4}.\]If \(\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) \leq \tau \leq t\) and \(\left \vert y\right \vert >\frac{\left \vert x\right \vert +{\rm c}t}{2}\), then\[\begin{align} \left \vert x-y\right \vert -{\rm c}\left( t-\tau \right) &\geq &\left \vert y\right \vert -\left \vert x\right \vert -{\rm c}\left( t-\tau \right) \geq \frac{\left \vert x\right \vert +{\rm c}t}{2}-\left \vert x\right \vert -{\rm c}t+\frac{{\rm c}t}{2}+\frac{1}{4}\left( {\rm c}t+\left \vert x\right \vert \right) \geq \frac{{\rm c}t-\left \vert x\right \vert }{4}. \end{align}\]Hence, \[\begin{align} J_{51} &\lesssim &\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\int_{\left \vert y\right \vert \leq \frac{\left \vert x\right \vert +{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N}\\ &&\cdot\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}\left( 1+\frac{\left( {\rm c}\tau-\left \vert y\right \vert \right) ^{2}}{1+t}\right) ^{-1}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}(1+t-\tau)^{-2}(1+t-\tau)^{\frac{5}{4}}(1+\tau)^{\frac{3}{4}}d\tau\\ &\lesssim &\left( 1+t\right) ^{-3}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1}, \end{align}\]\[\begin{align} J_{52} &\lesssim &\int_{\frac{t}{2}+\frac{1}{4}\big( t+\frac{\left \vert x\right \vert }{{\rm c}}\big) }^{t}\int_{\left \vert y\right \vert >\frac{\left \vert x\right \vert +{\rm c}t}{2}}\left( 1+t-\tau \right) ^{-2}\Big(1+\frac{({\rm c}t-|x|)^2}{1+t}\Big)^{-N/2}\\ &&\cdot\Big(1+\frac{(|x-y|-{\rm c}(t-\tau))^2}{1+t-\tau}\Big)^{-N/2}\left( 1+t\right) ^{-4}dyd\tau \\ &\lesssim &\left( 1+t\right) ^{-4}\Big(1+\frac{({\rm c}t-|x|)^2}{1+t}\Big)^{-N/2}\int_{\frac{t}{2}+\frac{1}{4}\left( t+\frac{\left \vert x\right \vert }{{\rm c}}\right) }^{t}\left( 1+t-\tau \right) ^{-2+\frac{5}{2}}d\tau \\ &\lesssim &\left( 1+t\right) ^{-\frac{5}{2}}\Big(1+\frac{({\rm c}t-|x|)^2}{1+t}\Big)^{-N/2}. \end{align}\]Gathering all the estimates, we have \[K_7\lesssim \left( 1+t\right) ^{-\frac{3}{2}}\bigg( 1+\frac{\left \vert x\right \vert ^{2}}{1+t}\bigg) ^{-\frac{3}{2}}+\left( 1+t\right) ^{-2}\bigg( 1+\frac{\left( {\rm c}t-\left \vert x\right \vert \right) ^{2}}{1+t}\bigg) ^{-1},\] and have completed the proof of (?? ).

The following two lemmas are usually employed to derive space-time estimates of the fluid models. In fact, Lemma 12 is used to deal with the term with -1-order Riesz operator in low frequency of Green funcion, Lemma 14 is used to deal with the terms containing double Riesz operator with the symbol \(\frac{\xi\xi^T}{|\xi|^2}\) in low frequency part of Green’s function. For readers’ convenience, we write down the proof since it was just stated in the previous works without a detailed proof.

Acknowledgments↩︎

Z.G. Wu was supported by Guangxi Natural Science Foundation \(\#\)​2025JJA110146. W.K. Wang was supported by National Natural Science Foundation of China \(\#\)​12271357 and Shanghai Science and Technology Innovation Action Plan \(\#\)​21JC1403600. Y.H. Zhang was supported by Guangxi Natural Science Foundation \(\#\)​2025JJF110001, \(\#\)​2024GXNSFDA010071, National Natural Science Foundation of China \(\#\)​12271114, Center for Applied Mathematics of Guangxi (Guangxi Normal University), and the Key Laboratory of Mathematical Model and Application (Guangxi Normal University), Education Department of Guangxi Zhuang Autonomous Region.

References↩︎