WELLPOSEDNESS AND DECAY RATES FOR THE CAUCHY PROBLEM OF THE MOORE-GIBSON-THOMPSON EQUATION ARISING IN HIGH INTENSITY ULTRASOUND

In this paper, we study the Moore–Gibson–Thompson equation in R , which is a third order in time equation that arises in viscous thermally relaxing fluids and also in viscoelastic materials (then under the name of standard linear viscoelastic model). First, we use some Lyapunov functionals in the Fourier space to show that, under certain assumptions on some parameters in the equation, an energy norm related with the solution decays with a rate (1 + t)−N/4. But this does not give the decay rate of the solution itself. Hence, in the second part of the paper, we show an explicit representation of the solution in the frequency domain by analyzing the eigenvalues of the Fourier image of the solution and writing the solution accordingly. We use this eigenvalues expansion method to give the decay rate of the solution (and its derivatives), which (for the solution) results in (1 + t)1−N/4 for N = 1, 2 and (1 + t)1/2−N/4 when N ≥ 3.


Introduction, derivation of the model and well-posedness.
Acoustic is an active field of research which is concerned with the generation and spacetime evolution of small mechanical perturbations in fluid (sound waves) or in solid (elastic waves). One of the important equations in nonlinear acoustics is the Kuznetsov equation: where u represents the acoustic velocity potential, and b, c, and B/A are the diffusivity and speed of sound, and the nonlinearity parameter, respectively. The derivation of equation (1.1) can be obtained from the general equations of fluid mechanics but we include a brief summary of its derivation, that can be found with more details in [12], [22] or [4] and the references therein, for instance. It is based on the relations coming from the conservation of mass, conservation of momentum and entropy balance in the model of thermo-viscous flow in a compressible fluid: • the equation of conservation of momentum (Newton's second law): • the equation of conservation of energy (first law of thermodynamics): (1.4) θ(η t + (v · ∇)η) = −∇ · q + T : D.
In the previous equations we are considering v, p and η as the velocity, pressure and specific entropy of the acoustic particle, and and q as mass density and the heat flux vector, respectively. Also, D is the deformation or strain tensor given by and T is the Cauchy-Poisson stress tensor given by where I is the identity matrix, µ is the shear viscosity (the first coefficient of viscosity) and λ = ζ − 2 3 µ, where ζ is the second coefficient of viscosity (the bulk viscosity). The components of T : D are T ij D ij , where T ij , D ij are the components of T and D, respectively.
It can be seen (see the above references for the details) that equations (1.3) and (1.4) can be rewritten as (1.5) (v t + (v · ∇)v) = −∇p + µ∆v + (ζ + µ/3)∇(∇ · v) and (1.6) θ(η t + (v · ∇)η) = 2µD : D + λ(∇ · v) 2 − ∇ · q, respectively. The previous equations, together with (1.2) and the following equations of state (1.7) p = p( , η), θ = θ( , η) are the Navier Stokes equations. First, we assume that the deviation of , p, η and θ from their equilibrium values 0 , p 0 , η 0 and θ 0 is small. By taking the Taylor series expansion of (1.7) around values at rest 0 and η 0 and ignoring the higher order terms, we get p( , η) = p( 0 , η 0 ) + ∂p ∂ ( 0 , η 0 ) ( − 0 ) + 1 2 We put where ∇p 0 = 0, χ is the coefficient of volume expansion and γ = c p /c v is the ratio of specific heat, as c p and c v are the specific heat capacities at constant pressure and constant volume. Then, the pressure p is given by By assuming that the flow is rotation free, that is ∇ × v = 0, and introducing the acoustic velocity potential v = −∇u, then it has been shown in [14] and [4], that equation (1.1), can be derived from the above set of equations by assuming the Fourier law of heat conduction (1.9) q = −K∇θ, where K is the thermal conductivity and θ is the absolute temperature. It is known that by modeling heat conduction with the Fourier law (1.9), which assumes the flux q to be proportional to the gradient of the temperature ∇θ at the same time t, leads to the paradox of infinite heat propagation speed (that is, any thermal disturbance at a single point has an instantaneous effect everywhere in the medium) and also fails when q increases or ∇θ decreases (see [10]). To overcome this drawback, a number of modifications of the basic assumption on the relation between the heat flux and the temperature have been made, such as the Maxwell-Cattaneo law, the Gurtin-Pipkin theory, the Jeffreys law, the Green-Naghdi theory and others. The common feature of these theories is that all permit transmission of heat flow as thermal waves at finite speed. See [2,11] for more details. One of these laws is the Maxwell-Cattaneo law, that assumes the following relation between heat flux and temperature: where τ is the relaxation time of the heat flux (usually small). By considering (1.10), instead of (1.9) and combining it with the equations of fluid mechanics, we get, instead of (1.1), the Jordan-Moore-Gibson-Thompson equation (see [10]) where b = δ + τ c 2 , with δ being the diffusivity of sound. In this paper we consider the linearized version of equation (1.11), known as the Moore-Gibson-Thompson equation in the acoustics theory: This linear equation is an active field of research. It also arises in viscous thermally relaxing fluids and has applications in medical and industrial use of high intensity ultrasound such as lithotripsy, thermotherapy or ultrasound cleaning (see [16]). But it also appears in viscoelasticity theory under the name of standard linear model of vicoelasticity (sometimes also called Kelvin or Zener model) to explain the behaviour of certain viscoelastic materials (that is, that exhibit both a viscous fluid and an elastic solid response) such as, for instance, fluids with complex microestructure (see [15]). In this context, u represents the linear deformations of a viscoelastic solid with an approach that is considered to be more realistic than the usual Kelvin-Voigt model. Actually, this model seems to be the simplest one that reflects both creeping and stress relaxation effects in viscoelastic materials (see [6] or [21] for more details). The derivation of equation (1.12) in R in the context of viscoelasticity theory can be obtained in the following way. Let us recall that in viscoelasticity theory, springs and dashpots represent the elastic and viscous components of the materials, respectively. According to [7], [1] or [6], among others, in the one dimensional case this equation represents a linear spring connected in series with a Kelvin-Voigt system, that is, another linear spring connected in parallel with a dashpot. This is a common way of approaching viscoelastic systems using a rheological point of view. Using this formulation (see again [7], [1] or [6] for more details), it is easy to see that such a system is governed by the following relation where σ is the stress and e is the strain. According to [6], τ is the stress relaxation time under constant strain, β is the strain relaxation time under constant stress and E is the relaxed elastic modulus. In any case, all of them are parameters involving the elastic and viscous constants of the material. More concretely, if η stands for the dashpot viscosity coefficient, and E 1 , E 2 represent the Young modulus of the first and second elastic springs respectively, one has that (1.14) However, there are a few references in which the standard linear model is described as a linear spring connected in parallel with a Maxwell model, that is, a spring and a dashpot connected in series (see for instance [18], where this model is also called the 3-parameter model). In this case, the relation between the stress and strain of the system is also (1.13), but the parameters are In both descriptions of the standard viscoelastic model, as relation (1.13) is the same, the corresponding equation would be where ρ represents the longitudinal density of the material. This equation is obtained thinking our material as a sequence of increasingly series-coupled systems of (1.13) type, that is, a continuous model with (1.13) as single component (see Chapter 6 of [5] or Section 2 of [19] for a similar deduction on different models). In [7], [1] or [20] it is assumed that 0 < τ < β (dissipative system), with τ, β being small constants. Observe that in both descriptions of the parameters (1.14) and (1.15) this is a natural assumption: in (1.14) it is fulfilled unless η = 0 (no dashpot), E 1 = 0 (only the Kelvin-Voigt sytem) or E 2 = ∞ (infinitely rigid second spring) ,while in (1.14) it is fulfilled unless η = 0 (no dashpot) or E 1 = ∞ (infinitely rigid first spring). In all these particular cases, we would obtain the conservative case τ = β. The case τ > β is treated in [3], where the authors show the chaotic behaviour of the corresponding solution.
The initial boundary value problem associated to (1.12) has been studied recently by many authors in bounded domains. In [12] (see also [13]), the authors considered the linearized equation (1.16) τ u ttt + αu tt + c 2 Au + bAu t = 0 where A is a positive self-adjoint operator, and showed that by neglecting diffusivity of the sound coefficient (b = 0) there arises a lack of existence of a semigroup associated with the linear dynamics. On the other hand, they showed that when the diffusivity of the sound is strictly positive (b > 0), the linear dynamics are described by a strongly continuous semigroup, which is exponentially stable provided that γ = α − τ c 2 /b > 0, while if γ = 0 the energy is conserved (the same type of results are obtained in [1] or [7] using energy methods, or in [17] using the analysis of the spectrum of the operator). The exponential decay rate results in [17] are completed in [20], where the obtention of an explicit scalar product where the operator is normal allows the authors to obtain the optimal exponential decay rate of the solutions. Finally, in [3], the authors show the caotic behaviour of the system when γ > 0, as we have mentioned above.
Observe that in this third order in time equation, the strong damping term bAu t is responsible for the well-posedness of the problem, while in the wave equation with strong damping (τ = 0) the strong damping term is responsible for the analiticity of the semigroup.
We also mention the recent paper [15] where the authors consider (1.16) with a memory damping term and show an exponential decay of the energy provided that the kernel is exponentially decaying. This result is generalized in [16], where it is shown that the memory kernel decay determines the solution decay.
To the best of our knowledge, these equations have not been studied yet in an unbounded domain. So, the goal of this paper is to show the well-posedness and investigate the decay rate of the solutions of the Moore-Gibson-Thompson equation in an unbounded domain. Namely, we consider the equation with the following initial data In order to state and prove our results, let us first and without loss of generality, take c = 1. In addition, we assume that 0 < τ < β, which corresponds to the dissipative case.
We first observe the well-posedness of (1.17) with initial data in H s (R N ), s ≥ 1. This can be seen using the result of [12] (also [17]) after applying the following change of variables to the equation (1.17): u = e −γt v with γ > 0. After this change of variables we obtain the equation Observe that γ i > 0 for i = 1, 2, 3 if γ is small enough. After rearranging terms, equation (1.19) can be written as Observe that, if γ is small enough, A = − ∆ − γ 3 γ 1 Id is a positive and self-adjoint operator in H 1 (R N ) and γ 1 , γ 2 > 0. Observe also that γ 4 − γ 2 γ 3 γ 1 v t can be seen as a bounded perturbation term. So, the hypotheses of [12] are satisfied and, hence, (1.17) is a well-posed problem in H 1 (R N ) (and, by linearity, in H s (R N ) for s ≥ 1). The two main results we obtain for the decay rate of this equation can be seen in Theorems 2.6 and 4.1 (and 4.2) below. First, in Section 2 and using the energy method in the Fourier space, we show that, if τ < β, the energy norm of V and of its higher-order derivatives ∂ j x V , with V = (τ u tt + u t , ∂ x (τ u t + u), ∂ x u t ), decay as for a certain c > 0 (see Theorem 2.6). The decay rate in (1.21) is a direct consequence of the estimate of the Fourier imageV (x, t): which we derive using the Lyapunov functional method (see Proposition 2.1). This method is a very powerful tool in proving the decay rate of the energy norm (see [8,9]). But, unfortunately, it is obvious that the estimate (1.21) does not give us any information about the decay of the solution u(x, t) itself. It only gives us the decay rate of the energy norm which eventually follows the decay rate of the slowest component of the vector V . So, to overcome this limitation of the Lyapunov functional method, we use the eigenvalues expansion method, which is based essentially on the behavior of the eigenvalues of the equation in the Fourier space. In Section 3 we give a complete description of the solutions of the characteristic equation of the corresponding operator and use it to divide the frequency domain into three main parts (low frequency, middle frequency and high frequency regions) and estimate the Fourier image of the solution in each region. In Section 4 and using this method, we are able to prove our second main result, which is the following decay rate of the solution of (1.17)-(1.18) when 0 < τ < β with initial data in (for certain c, C > 0 and 0 ≤ j ≤ s) and even better estimates if N + j ≥ 3 (see Theorem 4.2): For initial data in the weighted space L 1,1 (R N ) ∩ H s (R N ), s ≥ 1, the above estimates will be improved (see Theorem 4.4 for more details).
To summarize, the remaining part of this paper is organized as follows. In Section 2 we use the energy method in the Fourier space to build an appropriate Lyapunov functional, that is used to derive the decay rate of the energy norm explained above. Section 3 is devoted to the eigenvalues expansion method, that is used in Section 4 to derive the decay rate of the solution and its spacial derivatives.

Energy method in the Fourier space
In this section, we apply the energy method in the Fourier space to show the decay rate of the energy- First, we can write the problem in the Fourier space taking the Fourier transform of equation (1.17) and the initial data (1.18). We then obtain the following ODE initial value problem: the previous ODE can be rewritten as the following first order system We can write the previous system in a matrix form as with the initial dataÛ Thus, the pointwise estimate of the Fourier image of V reads as follows.
Proposition 2.1. Letû be the solution of (2.1)-(2.2). Assume that τ < β. Then, the Fourier image of the above vector V satisfies the estimate for all t ≥ 0 and certain c, C > 0, where The proof of Proposition 2.1 will be given through some lemmas, where a certain Lyapunov functional is obtained and used. First, we may rewrite system (2.3) as Lemma 2.2. The energy functional associated to system (2.8) is and satisfies, for all t ≥ 0, the identity Proof. Summing up the second and the third equation in (2.8) we get Multiplying (2.11) byv + τw and taking the real parts, we obtain, Next, multiplying the second equation in (2.8) by τ (β − τ )v and taking the real part, we get Now, multiplying the second equation in (2.8) by τ and adding the result to the first equation, we obtain (2.14) (û + τv) t = τŵ +v.
Next, we define the functional F 2 (ξ, t) as Proof. Multiplying the second equation in (2.8) by −τ (v + τw) and (2.11) by −τv, we obtain, respectively, Summing up the above two equations and taking the real parts, we obtain Applying Young's inequality, we obtain the estimate (2.19) for any 1 , 2 > 0.
Proof of Proposition 2.1. We define the Lyapunov functional L(ξ, t) as where γ 0 and γ 1 are positive numbers that will be fixed later on.
Taking the derivative of (2.20) with respect to t and making use of (2.10), (2.17) and (2.19), where we used the fact that |ξ| 2 /(1+|ξ| 2 ) ≤ 1. In the above estimate, we can fix our constants in such a way that the previous coefficients are positive. This can be achieved as follows: we pick 0 and 1 small enough such that 0 < 1 and 1 < 1. After that, we take γ 1 large enough such that Once γ 1 and 0 are fixed, we select 2 small enough such that Finally, and recalling that τ < β, we may choose γ 0 large enough such that Consequently, we deduce that there exists a positive constant γ 2 such that for all t ≥ 0, On the other hand, it is not difficult to see that from (2.20), (2.9), (2.16) and (2.18) and for γ 0 , large enough, that there exists two positive constants γ 3 and γ 4 such that Combining (2.22) and (2.23), we deduce that there exists a positive constant γ 5 such that for all t ≥ 0, A simple application of Gronwall's lemma, leads to the estimate (2.6), as L and the norm of V are equivalent.
In order to prove our first main result, we also need the first inequality of the following lemma. The rest of it will be used to prove the decay results in the other sections.
Lemma 2.5. For all t ≥ 0 and for all j ≥ 0, c > 0, the following estimates hold: Also, Moreover, Proof. First, to prove inequality (2.25) we will first prove that for given c > 0 and k ≥ 0, we have for all t ≥ 0, where C is a positive constant independent of t. To see this, observe first that for 0 ≤ t ≤ 1, the estimate (2.29) is obvious. On the other hand, for t ≥ 1, we have Now, by using (2.30) and the change of variables z = cr 2 t, where Γ is the gamma function. This yields (2.29 |ξ| j e −c|ξ| 2 t t 2 dξ where we have used (2.25). This inequality holds true for all N ≥ 1, in particular for N = 1, 2. But now, for N ≥ 3, we can improve the above estimate and get (2.28). Indeed, by taking the change of variable r = |ξ| and dξ = |ξ| N −1 dr, we write Now, for j + N − 3 ≥ 0 (which holds for all j ≥ 0 and N ≥ 3 or for all j + N ≥ 3 and N ≥ 1), we have Consequently, we have from above that , which is exactly (2.28).
We can now proceed to give and prove our first main result, which reads as follows.
Applying the Plancherel theorem and using the estimate in (2.6), we obtain Exploiting (2.32), we infer that (2.34) where we have used the inequality (2.25). In the high-frequency region (|ξ| ≥ 1), we have Collecting the above two estimates, we obtain (2.31). This finishes the proof of Theorem 2.6.
Remark 2.7. The estimate (2.31) does not give the decay rate of the solution u. In fact, it gives the decay rates of the norms τ u tt + u t L 2 , ∂ x (τ u t + u) L 2 and ∂ x u t L 2 (and also of the corresponding derivatives). These three norms are expected to decay with different rates. In order to know the decay rate of the solution and all its spatial derivatives, we need to use the explicit form of the solution and the eigenvalues expansion. This will be done in the following section.

Eigenvalues expansion
In this section, we use the eigenvalues expansion and the explicit form of the Fourier image of the solution in order to find the decay rates of the solution and its spacial derivatives.
The characteritic equation associated to (2.4) is The solutions λ i , i = 1, 2, 3 of the previous equation are the eigenvalues of Φ(ξ). We will use either λ i (ξ) or λ i (|ξ|) to denote them (depending on which of both notations is more convenient and when no confusion is possible) during the text below.
The following proposition on the description of these eigenvalues is an adaptation of some of the results of Proposition 4 of [20] (some part also in [17]), that we summarize and adapt here for a better comprehension and to be used later in the present work.
Remark 3.2. Observe that with the previous labelling of the eigenvalues in part 1 of Proposition 3.1, each λ i (|ξ|), i = 1, 2, 3, is a continuous function in ξ. However, during the rest of the paper and for simplicity in the notation we may call λ 1 the real eigenvalue and λ 2,3 the complex conjugate ones for all ξ in some proofs (it will be mentioned when this is done).
From the first equation, we have two possibilites: α = 0 or α = ± β τ |ξ|. In the first case and using now the second equation, the only possibility is that α = |ξ| = 0, which actually means that λ 0 = 0 is a (double) real eigenvalue when |ξ| = 0. The second one would be fulfilled only if |ξ| = 0 or if β = τ . If |ξ| = 0, we would again obtain λ 0 = 0 as a (double) real eigenvalue. The case β = τ will not be considered since we assumed that 0 < τ < β. Hence, the characteristic equation associated to (2.1) has no pure imaginary solutions in the dissipative case.
In order to give the decay rate of the solution in the next section, we now proceed to give asymptotic approximations of the eigenvalues of Φ(ξ) when |ξ| → 0 and |ξ| → ∞. For this purpose, it will be more convenient to apply the change of variables ζ = i|ξ| in the characteristic equation (3.1), that now becomes: Recall that λ j (ζ), j = 1, 2, 3, are the roots of (3.6), that we write j ζ 2 + ..., j = 1, 2, 3.
Remark 3.4. The behavior of the solution of (2.4) depends on the behavior of the function e Re λ j (ξ)t , j = 1, 2, 3. Since in most cases λ j (|ξ|) is a power series of |ξ|, so it is its real part.
Observe that as Re λ j (|ξ|) < 0, the frequencies that give the dominant part of all e Re λ j (|ξ|)t are those corresponding to small frequencies |ξ|. For this reason, the behavior of the real part near |ξ| = 0 determines the decay rate of the solution. For large frequencies, and again as Re λ j (|ξ|) < 0, it is clear that e Re λ j (|ξ|)t can be always estimated by e −ct if the powers in the Taylor series expansion of Re λ j (|ξ|) near infinity are positive or by |ξ| m e −c|ξ| −ĉ t for a certain m > 0 if one of the powers in the Taylor series expansion is negative. In both cases and using Plancherel's theorem, we see that the integral in the high frequencies is bounded if and only if some derivatives of the solution are bounded, which means that the solution should be in some Sobolev spaces and this gives the regularity of initial data needed for the desired decay rate.
Let us now divide the frequency space into three regions: low frequency, high frequency and middle frequency region, that is The choice of ν 1 and ν 2 will be discussed in the proofs of Proposition 3.5 and of Lemma 3.8. For the moment, we need ν 1 and ν 2 sufficiently small and large, respectively, such that the asymptotic expansions of Propositions 3.5 and 3.7 hold.
We write the solution of the system (2.1) in the above three regions. In the following Propositions 3.5, 3.7 and 3.9 we give bounds of the solution on each of this three regions using the previous asymptotic expansions of the eigenvalues. These bounds will be used in Theorem 4.1 to proof the decay estimate of the solution of problem (1.17).
Proposition 3.5. If 0 < τ < β, the solutionÛ (ξ, t) of (2.4) satisfies, for all ξ ∈ Υ L with |ξ| = 0, the estimates: where L 1,1 is the L 1 -weighted space defined by Remark 3.6. Observe that the estimate (3.13) is not satisfied if |ξ| = 0 but, as it is a set of measure zero, it will not affect the decay of the solution in Theorem 4.1.
Proof. According to part 1 of Proposition 3.1, if ξ ∈ Υ L is such that |ξ| is small enough (that is, |ξ| < ν 1 < √ m 1 ) we know that there exist one real root and two complex conjugate ones.
For simplicity, we are going to denote λ 1 the real root and λ 2,3 the complex conjugate ones when |ξ| → 0, both when 0 < τ /β < 1/9 and 1/9 < τ /β < 1 (see Remark 3.2 and Proposition 3.1). Hence, the solution of the equation (2.1) when ξ ∈ Υ L can be written in terms of the corresponding eigenvalues as We may use the initial values (2.2) to find the above constants by solving the system (we are omitting the ξ dependence in order to simplify the notation). By neglecting the small terms, we have from (3.9) that in Υ L the eigenvalues are: Solving the system (3.17) and using the asymptotic expressions of the eigenvalues when |ξ| → 0 in (3.18) we get that, in Υ L , We can now take the following approximate solution of (2.1) and (2.2), Observe that, at a first leading order, solving the system (3.17) is equivalent to solving which has the previous C 1 , C 2 , C 3 as exact solution.
where D i (ξ) can be written in terms of the initial data (2.2) and hence satisfy the same system as in Proposition 3.5, (3.17). From (3.12) and by neglecting the small terms, we also know that, when |ξ| → ∞, Hence, observe that at a first leading order, solving the corresponding system is equivalent to solving Solving the corresponding system and using the asymptotic expressions of the eigenvalues when |ξ| → ∞ in (3.22) we get that in Υ H In a similar way as in Proposition 3.5, we deduce that (3.20) holds.
This last inequality together with (3.23) and (3.28) leads to (3.25). Second, if the roots of (3.1) are real and distinct, then the solution of (2.1) is written as where C 1 , C 2 and C 3 are satisfying the system λ 2 1 C 1 + λ 2 2 C 2 + λ 2 3 C 3 =û 2 (that this system was obtained imposing thatû(ξ, t) must satisfy the initial conditions (2.2)). It is not hard to see that (3.29) also holds and therefore (as before) (3.25) is also satisfied.

Decay estimates
In this section, we show the decay estimates for the solution u(x, t) of the system (1.17)-(1.18) using the eigenvalues expansion results of Section 3. For simplicity in the notation, all the norms of this section will be in R N , although we are going to skip this in the notation.
Then for any t ≥ 0 the following decay estimates hold for all 0 ≤ j ≤ s and certain constants C, c > 0 independent of t and of the initial data: where c = min 1 β , β−τ 2βτ , | Re(λ 2,3 (ξ ν ))| . The above estimates can be improved for if N + j ≥ 3 and get the following Theorem.
Remark 4.3. It is clear that for N = 3, the estimate (4.2) improves the one in (4.1). Indeed, from (4.1), we deduce in this case that the L 2 -norm of the solution does not decay if j = 0. On the other hand and from (4.2) we deduce that, for N = 3, the L 2 -norm of the solution decays with the rate (1 + t) −1/4 . Also, for N = 1 and j ≥ 2 or N = 2 and j ≥ 1, the estimate (4.2) gives faster decay rate than (4.1).
Observe that, as we had already pointed out in the previous Remarks 3.6 and 3.10, the values of ξ where we have double or triple real roots are sets of measure zero and, hence, they do not affect the decay of the solution.
As we have said, the estimates in Υ M and Υ H remain the same as in the proof of Theorem 4.1.