Multiharmonic Algorithms for Contrast-Enhanced Ultrasound

TL;DR

This paper introduces new algorithms to efficiently model nonlinear effects in contrast-enhanced ultrasound using multiharmonic methods.

Contribution

The paper proposes multiharmonic algorithms that reduce computational costs in modeling microbubble dynamics in ultrasound.

Findings

The existence of time-periodic solutions for the Westervelt-ODE system is rigorously established.

Multiharmonic algorithms are derived for the system under time-periodic excitation.

Numerical investigations show how harmonics and microbubbles influence acoustic wave propagation.

Abstract

Harmonic generation plays a crucial role in contrast-enhanced ultrasound, both for imaging and therapeutic applications. However, accurately capturing these nonlinear effects is computationally demanding when using traditional time-domain approaches. To address this issue, we develop algorithms based on a time discretization that uses a multiharmonic Ansatz applied to a model that couples the Westervelt equation for acoustic pressure with a volume-based approximation of the Rayleigh–Plesset equation for the dynamics of microbubble contrast agents. We first rigorously establish the existence of time-periodic solutions for this Westervelt-ODE system. We then derive a multiharmonic representation of the system under time-periodic excitation and develop iterative algorithms that rely on the successive computation of higher harmonics assuming either real-valued or complex-valued solution fields. In the real-valued setting, we characterize the approximation error in terms of the number of harmonics and a contribution arising from the fixed-point iteration. Finally, we investigate these algorithms numerically and illustrate how the number of harmonics and the presence of microbubbles influence the propagation of acoustic waves.

Figures

Overview of multiharmonic algorithms in this work

Relative error (34) of the quantity of interest with respect to the mesh size $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }$$\end{document}$ on a log–log scale

Five-harmonic expansion $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {R}(p_{h_{FEM }}^{(5)}(x_1,x_{2,0},t_0))$$\end{document}$ plotted as a function of $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_1$$\end{document}$

Relative error (34) of the quantity of interest with respect to the truncation level *N* on a semi-log scale (logarithmic *y*-axis)

Pointwise differences $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|p_{h_{FEM }}^N(x,t_0) - p_{h_{FEM }}^{(10)}(x,t_0) |$$\end{document}$ for $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^5$$\end{document}$ at a fixed time $

Pressure field $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Re (p_{h_{FEM }}^{(5)}(x,t_0))$$\end{document}$ at a fixed reference time $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document}$ for $\documentclass[12pt]{m

Real part of the pressure from the five-harmonic expansion computed using (28) without bubbles ( $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0=0$$\end{document}$ , dashed green line) and with bubbles (blue line) over time at three spatial pointsFig. 8Real part of the pressure $\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy

Full text

Introduction

Contrast-enhanced ultrasound has become an important tool in biomedical applications, with gas-filled microbubbles being used to improve both diagnostic and therapeutic procedures. While ultrasound propagation can be described by linear acoustic models in the small-amplitude regime, nonlinear effects arise at higher acoustic pressures, leading to the generation of higher harmonics in the frequency domain. These nonlinearities are further amplified by microbubble contrast agents, which exhibit strongly nonlinear oscillatory behavior when exposed to ultrasound waves. This delicate back-and-forth interaction is beneficial for improving resolution in imaging but also for enhancing therapeutic treatments, such as targeted drug delivery; see, for example, [12, 18, 37] for details. With the rise in the number of applications of contrast-enhanced ultrasound, accurate modeling and efficient simulation in this context have also become prominent research topics; see, e.g., [6, 10, 28, 39], and the references given therein.

The present work builds upon [33], where time-domain mathematical models for ultrasound contrast imaging with microbubbles based on a nonlinear acoustic wave equation coupled to a Rayleigh–Plesset-type ODE in a continuum (effective-medium) description of the bubbly mixture have been derived and investigated in terms of local well-posedness and numerical simulations. As noted in [33], a major computational challenge when simulating such systems stems from different time scales on which the wave equation and ODE (nonlinearly) evolve. As a result, straightforward numerical approaches for solving such systems demand using prohibitively small time steps. In this work, we approach this issue by developing multiharmonic algorithms for time-periodic solution fields, which offer a potentially more efficient modeling alternative by making use of harmonic expansions; see, e.g., [3, 4], where multiharmonic ideas have been explored for problems arising in electromagnetism and, e.g., [15, 21, 22, 34], where they have been developed for single-physics acoustic models. Instead of resolving every oscillation in the space-time domain, these methods decompose the field into a sum of harmonics, which can be obtained as solutions of suitable (in our case) Helmholtz problems and algebraic equations. These methods are especially promising for simulating real-time ultrasound applications in which the number of harmonics used can be relatively low.

As the starting time-domain model of contrast-enhanced ultrasound we employ the following coupled nonlinear wave-ODE system:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \boxed { \begin{aligned}&\ \ p_{tt} - c^2 \Delta p - b \Delta p_t = \eta (p^2)_{tt} + c^2 \rho _0 n_0(x) v_{tt}+ h(x,t) \quad & \text {in } \Omega \times (0,T), \\&\ \ v_{tt} + \delta \omega _0 v_t + \omega _0^2 v = \zeta v^2 + \xi (2v v_{tt}+ v_t^2)- \mu p \quad & \text {in } \Omega \times (0,T), \\ \end{aligned} } \end{aligned}$$\end{document}

consisting of the damped Westervelt equation [40] for the acoustic pressure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p=p(x,t)$$\end{document} (that is, fluctuations in the background pressure) and an ODE pointwise a.e. in space for the volume variation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v=v(x,t)$$\end{document} of microbubbles. The system in (1) models the interaction of acoustic pressure waves with oscillations of microbubbles.

Modeling Background

In the Westervelt equation, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c>0$$\end{document} denotes the (constant) speed of sound in the medium, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b>0$$\end{document} the diffusivity of sound, so that the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-b \Delta p_t$$\end{document} introduces strong damping, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0>0$$\end{document} the mass density of the mixture at equilibrium. The nonlinearity coefficient is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta = \frac{\beta _a}{\rho _0 c^2}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _a$$\end{document} is the nonlinearity parameter in the medium. The source of the pressure waves is provided in part through a contribution due to bubble oscillations, modeled by the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^2 \rho _0 n_0 v_{tt}$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0=n_0(x)$$\end{document} is the bubble number density at equilibrium, and by the source function \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=h(x,t)$$\end{document} . We allow \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0$$\end{document} to vary in space (that is, we assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0 \in L^\infty (\Omega )$$\end{document} ), as this setting is relevant for imaging applications; see, e.g., [25]. Moreover, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0$$\end{document} and h regulate the strength of the acoustic source and will be assumed to be small enough in the well-posedness analysis; see Theorem 1 and the discussion in Section 2.3 on the physical meaning of the smallness assumption.

The ODE in (1) captures nonlinear harmonic oscillations driven by the acoustic pressure via the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$- \mu p$$\end{document} . The total microbubble volume is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$V= v_0 + v$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0$$\end{document} is the equilibrium volume. The coefficient \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta = \frac{4 \nu }{\omega _0 R_0^2}$$\end{document} is the viscous damping coefficient and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _0 = \sqrt{\frac{3 \kappa P_0}{\rho _0 R_0^2}}$$\end{document} the natural frequency, where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} is the bubble radius at equilibrium volume \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0$$\end{document} (that is, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0 = \frac{4 \pi }{3} R_0^3$$\end{document} ), \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} the adiabatic exponent and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_0$$\end{document} the ambient pressure in the mixture (so that the total pressure is \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_0 +p$$\end{document} ). The coefficients appearing on the right-hand side of the ODE are given in terms of equilibrium values by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = \frac{4 \pi R_0}{\rho _0}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta = \frac{(\kappa + 1) \omega _0^2}{2 v_0}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi = \frac{1}{6 v_0}$$\end{document} . As discussed in [17, Ch. 5], the ODE in (1) can be seen as a volume-based approximation of the following Rayleigh–Plesset equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \rho _0 \left[ R R_{tt} + \tfrac{3}{2} R_t^2 \right] = p_b - 4 \nu \frac{R_t}{R} - p, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} is the kinematic viscosity and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_b$$\end{document} a constant pressure contribution, which can be derived using the relation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ V= \frac{4 \pi }{3} R^3$$\end{document} and the adiabatic gas law \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{p_b}{P_0} = \left( \frac{v_0}{V}\right) ^{\kappa }$$\end{document} . Introducing the volume variable and expanding the resulting expression about the equilibrium radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} , while neglecting higher-order terms, yields a simplified ODE for the volume perturbation v. A detailed derivation of this approximation is provided in [35, Sec. 2.3]. Reformulating the dynamics in terms of v instead of R removes singular terms and results in a more manageable nonlinear structure. Since pressure perturbations couple more directly to volume oscillations, the volume-based formulation is also better suited for harmonic expansions. In contrast, a harmonic expansion in terms of R is expected to lead to higher-order interactions that complicate the isolation of harmonic components due to the cubic dependence of volume on radius.

We equip (1) with absorbing-type boundary conditions of the following form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \beta p_t + \gamma p + \nabla p \cdot \boldsymbol{n}= 0 \quad \text { on } \partial \Omega \times (0, T), \quad \beta ,\, \gamma >0, \end{aligned}$$\end{document}

and we are interested in time-periodic solutions that satisfy

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} p(0)=p(T), \qquad p_t(0)=p_t(T) \qquad \text { in } \Omega , \\ v(0)=v(T), \qquad v_t(0)=v_t(T) \qquad \text { in } \Omega , \end{aligned} \end{aligned}$$\end{document}

which corresponds to steady-state oscillatory (stable cavitation) regimes rather than transient cavitation. Although in the existence analysis we do not require the source h to be time periodic, for developing multiharmonic algorithms, we assume that h represents the second time derivative of a T-periodic function g; see Section 3 for details.

Main Contributions

The overall purpose of the present work is to establish multiharmonic approximation approaches for time-domain systems in the form of (1), that is, systems that incorporate both nonlinear acoustic propagation and nonlinear microbubble oscillations. To approximate the problem we employ a cut-off Fourier series for the pressure and volume fields:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} p \approx p^N= \mathfrak {R}\left\{ \sum _{m=0}^{N} \exp (\imath m \omega t) p_m^N(x)\right\} , \quad v \approx v^N= \mathfrak {R}\left\{ \sum _{m=0}^{N} \exp (\imath m \omega t) v_m^N(x)\right\} , \end{aligned} \end{aligned}$$\end{document}

with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = \dfrac{2 \pi }{T}$$\end{document} . In addition to rigorously considering the setting of real-valued pressure-volume fields, we also discuss multiharmonic algorithms derived from complex-valued fields (that is, from dropping the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {R}$$\end{document} operator in (2)). This approach is sometimes adopted in electromagnetism (see, e.g., [9, 41]) and it has been discussed in [21] for the de-coupled Westervelt equation. Here, it results in simplified algorithms, which are, however, only formally investigated; see Section 5 for details. Toward reaching the overall goal of the work, our main contributions pertain to

• rigorously establishing the existence of time-periodic solutions for the Westervelt-ODE system in (1);
• deriving cutoff multiharmonic approximations of the system under time-periodic excitation for computing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N_m, v^N_m)$$\end{document} ;
• developing and analyzing linearized multiharmonic cut-off algorithms for computing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N_m, v^N_m) $$\end{document} . In particular, we characterize the error of linearized multiharmonic algorithms for real-valued fields in terms of the number of harmonics and a contribution due to the fixed-point iteration; see Theorem 2. For convenience, an overview of algorithms investigated in this work is provided in Figure 1.Fig. 1. Overview of multiharmonic algorithms in this work

Novelty and Related Work

To the best of our knowledge, this is the first work to rigorously develop multiharmonic approaches for coupled nonlinear wave-ODE systems of this type, providing both well-posedness theory and numerical analysis. An approach to second harmonic generation for a linearization of (1) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta =0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=0$$\end{document} in the wave equation has been formally set up in [17, Ch. 5] without convergence guarantees.

Rigorous multiharmonic studies have been performed on single-physics equations. Existence of periodic solutions of nonlinear acoustic models, including the Westervelt equation, has been investigated rigorously in [7, 8, 21–23, 34]. An iterative multiharmonic algorithm for the Westervelt equation has been proposed and rigorously studied in [34]. A numerical algorithm based on using the complex Fourier Ansatz for the de-coupled Westervelt equation and a boundary element approach for the resulting Helmholtz problems has been developed and investigated in [15]. In the context of electromagnetism, a multiharmonic treatment of the quasi-stationary Maxwell problem has been developed and rigorously analyzed in [3]; numerical simulation aspects are discussed in [4]. These results, however, do not address the coupling to microbubble/ODE dynamics.

Periodic solutions of Rayleigh–Plesset-type equations that do not incorporate modeling of the acoustic propagation have been rigorously investigated in the literature; we refer to, for example, [16, 42] and [38, Ch. 9], and the references provided therein. We mention in passing that the mathematical literature on non-periodic models in nonlinear acoustics of non-bubbly media is quite rich; see, e.g., [1, 11, 20, 24] and the references provided therein.

Organization of the Paper

The remainder of the exposition is organized as follows. In Section 2, we determine sufficient conditions for the well-posedness of the time-periodic boundary value problem for the Westervelt-ODE system (1). In Section 3, we derive a multiharmonic cut-off representation of the system and propose a linearization which results in a simplified setting for computing Fourier coefficients in (2). In Section 4, we characterize the error and prove the convergence of the proposed linearized multiharmonic scheme to the solution of the Westervelt-ODE system as the number of harmonics N tends to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\infty $$\end{document} . The main theoretical result of the section is contained in Theorem 2. In Section 5, we discuss multiharmonic algorithms resulting from complex solutions fields. Finally, in Section 6, we investigate and compare the introduced algorithms.

Existence of Time-Periodic Solutions

In this section, we establish the basis for the numerical analysis by investigating the existence and uniqueness of periodic solutions for the coupled PDE-ODE system. More precisely, given \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0$$\end{document} and a bounded domain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega \subset \mathbb {R}^d$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \in \{2,3\}$$\end{document} , we study the problem:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} \ &p_{tt}- c^2 \Delta p - b \Delta p_t= \eta (p^2)_{tt}+c^2 \rho _0 n_0(x) v_{tt}+ h \quad & \text {in } \Omega \times (0,T),\\&\beta p_t+ \gamma p + \nabla p \cdot \boldsymbol{n}= 0 & \text {on } \partial \Omega \times (0,T), \\&p(0) = p(T), \ p_t(0) = p_t(T) & \text {in } \Omega , \\&v_{tt}+ \delta \omega _0 v_t+ \omega _0^2 v = \zeta v^2 + \xi (2v v_{tt}+ v_t^2)- \mu p \quad & \text {in } \Omega \times (0,T),\\&v(0) = v(T), \ v_t(0) = v_t(T) & \text {in } \Omega . \end{aligned} \right. \end{aligned}$$\end{document}

The analysis will be based on successive approximations of the system, for which knowledge of linear time-periodic ODE and wave problems will be very helpful. We thus discuss those results next.

Notation. Below we use \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$lhs \lesssim rhs $$\end{document} to denote \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$lhs \le C \cdot rhs $$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C>0$$\end{document} is a generic constant. When writing norms in Bochner spaces, we omit the temporal domain (0, T). For example, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \cdot \Vert _{L^p(L^q(\Omega ))}$$\end{document} denotes the norm in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^p(0,T; L^q(\Omega ))$$\end{document} .

Auxiliary Existence Results for Linear Time-Periodic Problems

To set up the local well-posedness analysis of (3), we first state two separate results on the well-posedness of an ODE (describing damped oscillations with forcing) and a linear time-periodic wave problem.

Lemma 1

Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f \in L^2(0,T; L^\infty (\Omega ))$$\end{document} . Furthermore, let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _0>0$$\end{document} . Then, the periodic ODE problem

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned}&v_{tt}+ \delta \omega _0 v_t+ \omega _0^2 v = f & \text {a.e.\ in } \Omega \times (0,T),\\&v(0) = v(T), \ v_t(0) = v_t(T) & \text {a.e.\ in } \Omega , \end{aligned} \right. \end{aligned}$$\end{document}

has a unique solution

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} v \in \mathcal {X}_v=\{v \in H^2(0,T; L^\infty (\Omega )): \, v(0)=v(T), \, v_t(0)= v_t(T) \ \text {a.e.}\} \end{aligned} \end{aligned}$$\end{document}

that satisfies

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert v\Vert _{\mathcal {X}_v} :=\Vert v\Vert _{L^\infty (L^\infty (\Omega ))}+ \Vert v_t\Vert _{L^\infty (L^\infty (\Omega ))}+\Vert v_{tt}\Vert _{L^2(L^\infty (\Omega ))} \lesssim \Vert f\Vert _{L^2(L^\infty (\Omega ))}. \end{aligned} \end{aligned}$$\end{document}

Proof

The existence and uniqueness can be established using Floquet theory (see, for example [36] and [14, Ch. 3]) as the only solution of the homogeneous problem is zero due to the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta \omega _0>0$$\end{document} . The details are provided in Appendix A for completeness. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Proposition 1

(see [34, Theorem 2.1]) Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega \subset \mathbb {R}^d$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \in \{2,3\}$$\end{document} , be a bounded domain with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^{1,1}$$\end{document} boundary, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta ,\,\gamma >0$$\end{document} , c, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b>0$$\end{document} . Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h\in L^2(0,T;L^2(\Omega ))$$\end{document} . Then there exists a unique (weak) solution

$$[eqn]$$

of the time-periodic boundary-value problem

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned}&p_{tt}- c^2 \Delta p - b \Delta p_t= h\quad & \text {in } \Omega \times (0,T),\\&\beta p_t+ \gamma p + \nabla p \cdot \boldsymbol{n}= 0 & \text {on } {\partial \Omega \times (0,T)}, \\&p(0) = p(T), \ p_t(0) = p_t(T) & \text {in } \Omega . \end{aligned} \right. \end{aligned}$$\end{document}

The solution satisfies

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert p\Vert _{\mathcal {X}_p}:=&\, ||p||_{L^2(H^2(\Omega ))} + ||p_t||_{L^2(H^{3/2}(\Omega ))}+\Vert p_{tt}\Vert _{L^2(L^2(\Omega ))} + \Vert \Delta p_t\Vert _{L^2(L^2(\Omega ))} \\&\quad + \Vert \nabla p \cdot \boldsymbol{n}\Vert _{H^1(L^2(\partial \Omega ))} \\ \le&\,{ C(b, \gamma , c, \beta , T, \Omega ) ||h||_{ L^2(0,T;L^2(\Omega ))}}, \end{aligned} \end{aligned}$$\end{document}

for some constant \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C=C(b, \gamma ,c,\beta , T, \Omega ) > 0$$\end{document} , independent of h.

We note that having \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma >0$$\end{document} in (4) is crucial for uniqueness. Furthermore, having strong damping (that is, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-b \Delta p_t$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b>0$$\end{document} ) in the wave equation contributes to its parabolic-like character, and is heavily exploited in the proof of the above statement. We also note that the functions f and h do not have to be time periodic for the well-posedness results in this section, although we will make that assumption on the acoustic source in Section 3 to develop (iterative) multiharmonic algorithms.

Analysis of the Westervelt-ODE System

Equipped with the previous results on linear problems, we are now ready to prove the well-posedness of the coupled pressure-volume problem. To state it, we introduce

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathbb {B}_{r_p, r_v}= \left\{ (p, v) \in \mathcal {X}_p\times \mathcal {X}_v:\, \Vert p\Vert _{\mathcal {X}_p} \le r_p, \ \Vert v\Vert _{\mathcal {X}_v} \le r_v\right\} , \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v>0$$\end{document} will be set as small enough in the course of the proof below. The proof will have constructive nature using successive approximations.

Theorem 1

Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$T>0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega \subseteq \mathbb {R}^d$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$d \in \{2,3\}$$\end{document} , be open, bounded, connected, with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C^{1,1}$$\end{document} boundary, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma $$\end{document} , c, b, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _0>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta $$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi \in \mathbb {R}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0 \in L^\infty (\Omega )$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h \in L^2(0,T; L^2(\Omega ))$$\end{document} . Then there exist \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _p>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta _v>0$$\end{document} , such that if

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} ||h||_{L^2(0,T;L^2(\Omega ))} \le \delta _p\quad \text {and} \quad \Vert {n_0}\Vert _{L^\infty (\Omega )} \le \delta _v, \end{aligned}$$\end{document}

then there exist radii \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v>0$$\end{document} , such that problem (3) has a unique solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p, v) \in \mathbb {B}_{r_p, r_v}$$\end{document} .

Proof

As announced, the proof follows by constructing the solutions using successive approximations. Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v>0$$\end{document} . We take \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^{(0)}, v^{(0)}) \in \mathbb {B}_{r_p, r_v}$$\end{document} and set up successive approximations using the system given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} \ &p^{(N)}_{tt} - c^2 \Delta p^{(N)}- b \Delta p^{(N)}_t = \eta ((p^{(N-1)})^2)_{tt}+c^2 \rho _0 n_0 v^{(N)}_{tt} + h \, & \text {in } \Omega \times (0,T),\\&\beta p^{(N)}_t + \gamma p^{(N)}+ \nabla p^{(N)}\cdot \boldsymbol{n}= 0 & \text {on } \partial \Omega \times (0,T), \hspace{-0.28cm}\\&p^{(N)}(0) = p^{(N)}(T), \ p^{(N)}_t(0) = p^{(N)}_t(T) & \text {in } \Omega , \end{aligned} \right. \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned}&v^{(N)}_{tt} + \delta \omega _0 v^{(N)}_t + \omega _0^2 v^{(N)} & \\ =&\, - \mu p^{(N-1)}+\zeta (v^{(N-1)})^2 - \xi \left( 2 v^{(N-1)}v^{(N-1)}_{tt} + (v^{(N-1)}_t)^2 \right) & \text {in } \Omega \times (0,T),\\&v^{(N)}(0) = v^{(N)}(T), \ v^{(N)}_t(0) = v^{(N)}_t(T) & \text {in } \Omega . \end{aligned} \right. \end{aligned}$$\end{document}

We proceed through several steps.

(i) Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=1$$\end{document} . Note that we can first solve (6b) and use its solution as the input for solving (6a). By Lemma 1, we have a unique \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^{(1)}\in \mathcal {X}_v$$\end{document} that solves (6b), such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert v^{(1)}\Vert _{\mathcal {X}_v} \lesssim \Vert p^{(0)}\Vert _{L^2(L^\infty (\Omega ))}+ \Vert v^{(0)}\Vert ^2_{\mathcal {X}_v}. \end{aligned} \end{aligned}$$\end{document}

From here, we find that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert v^{(1)}\Vert _{\mathcal {X}_v} \lesssim r_p+ r_v^2, \end{aligned} \end{aligned}$$\end{document}

and we can conclude that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert v^{(1)}\Vert _{\mathcal {X}_v} \le r_v$$\end{document} provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} are small enough. By Proposition 1, the linear problem in (6a) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=1$$\end{document} has a unique solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{(1)}\in \mathcal {X}_p$$\end{document} , such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert p^{(1)}\Vert _{\mathcal {X}_p} \lesssim&\, \Vert -\eta ((p^{(0)})^2)_{tt} - c^2 \rho _0 n_0 {v^{(1)}_{tt}} + h\Vert _{L^2(L^2(\Omega ))} \\ \lesssim&\, \Vert p^{(0)}\Vert ^2_{\mathcal {X}_p} + \Vert n_0\Vert _{L^\infty (\Omega )}\Vert v^{(1)}\Vert _{\mathcal {X}_v} + \Vert h\Vert _{L^2(L^2(\Omega ))}. \end{aligned} \end{aligned}$$\end{document}

By using the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert v^{(1)}\Vert _{\mathcal {X}_v} \le r_v$$\end{document} , we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert p^{(1)}\Vert _{\mathcal {X}_p} \lesssim r_p^2 + \Vert n_0\Vert _{L^\infty (\Omega )} r_v+ \Vert h\Vert _{L^2(L^2(\Omega ))}. \end{aligned} \end{aligned}$$\end{document}

From here, we can conclude that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert p^{(1)}\Vert _{\mathcal {X}_p} \le r_p$$\end{document} for small enough \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert h\Vert _{L^2(L^2(\Omega ))}$$\end{document} by reducing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert n_0\Vert _{L^\infty (\Omega )}$$\end{document} .

(ii) Let now \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^{(N-1)}, v^{(N-1)}) \in \mathbb {B}_{r_p, r_v}$$\end{document} . By Lemma 1, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert v^{(N)}\Vert _{\mathcal {X}_v} \lesssim \Vert p^{(N-1)}\Vert _{L^2(L^\infty (\Omega ))}+ \Vert v^{(N-1)}\Vert ^2_{\mathcal {X}_v} \end{aligned} \end{aligned}$$\end{document}

and by Proposition 1, the linear problem in (6a) with periodic time conditions has a unique solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{(N)}\in \mathcal {X}_p$$\end{document} , such that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert p^{(N)}\Vert _{\mathcal {X}_p} \le&\, C \Vert -\eta ((p^{(N-1)})^2)_{tt} - c^2 \rho _0 n_0 v^{(N)}_{tt} + h\Vert _{L^2(L^2(\Omega ))}\\ \lesssim&\, \Vert p^{(N-1)}\Vert ^2_{\mathcal {X}_p} + \Vert n_0\Vert _{L^\infty (\Omega )}\Vert v^{(N)}\Vert _{\mathcal {X}_v} + \Vert h\Vert _{L^2(L^2(\Omega ))}. \end{aligned} \end{aligned}$$\end{document}

By then proceeding analogously to the case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=1$$\end{document} , we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert p^{(N)}\Vert _{\mathcal {X}_p} \le r_p, \quad \Vert v^{(N)}\Vert _{\mathcal {X}_v} \le r_v, \end{aligned} \end{aligned}$$\end{document}

provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert n_0\Vert _{L^\infty (\Omega )}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert h\Vert _{L^2(L^2(\Omega ))}$$\end{document} are small enough. We thus conclude that for small enough \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert h\Vert _{L^2(L^2(\Omega ))}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert n_0\Vert _{L^\infty (\Omega )}$$\end{document} , there exist \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p>0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v>0$$\end{document} , such that if \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^{(0)}, v^{(0)}) \in \mathbb {B}_{r_p, r_v}$$\end{document} , then for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , problem (6) has a unique solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^{(N)}, v^{(N)}) \in \mathbb {B}_{r_p, r_v}$$\end{document} .

(iii) We next wish to show that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{(p^{(N)}, v^{(N)})\}_{N \ge 1}$$\end{document} is a Cauchy sequence in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {X}_p\times \mathcal {X}_v$$\end{document} . We note that for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N_2 \in \mathbb {N}$$\end{document} , the difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\overline{p}, \overline{v})=(p^{(N_1)}-p^{(N_2)}, v^{(N_1)}-v^{(N_2)})$$\end{document} solves

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} \ &\overline{p}_{tt} - c^2 \Delta \overline{p}- b \Delta \overline{p}_t & \\ =&\, -\eta ((p^{(N_1-1)}-p^{(N_2-1)})(p^{(N_1-1)}+p^{(N_2-1)}))_{tt} - c^2 \rho _0 n_0 (v^{(N_1)}_{tt}-v^{(N_2)}_{tt}) \, & \text {in } \Omega ,\\&\beta \overline{p}_t + \gamma \overline{p}+ \nabla \overline{p}\cdot \boldsymbol{n}= 0 & \text {on } \partial \Omega , \hspace{-0.28cm}\\ \end{aligned} \right. \end{aligned}$$\end{document}

and

$$[eqn]$$

From here, similarly to before, using Proposition 1 and Lemma 1, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert p^{(N_1)}-p^{(N_2)}\Vert _{\mathcal {X}_p} \lesssim r_p\Vert p^{(N_1-1)}-p^{(N_2-1)}\Vert _{\mathcal {X}_p} + \Vert n_0\Vert _{L^\infty (\Omega )}\Vert v^{(N_1)}-v^{(N_2)}\Vert _{\mathcal {X}_v} \end{aligned} \end{aligned}$$\end{document}

and

$$[eqn]$$

By adding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(8)+ \lambda \cdot (9)$$\end{document} , with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda >0$$\end{document} , we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\Vert p^{(N_1)}-p^{(N_2)}\Vert _{\mathcal {X}_p} + (\lambda - C \Vert n_0\Vert _{L^\infty (\Omega )})\Vert v^{(N_1)}-v^{(N_2)}\Vert _{\mathcal {X}_v}\\ \le&\, C ( \lambda + r_p) \Vert p^{(N_1-1)}-p^{(N_2-1)}\Vert _{\mathcal {X}_p} + C \lambda r_v\Vert v^{(N_1-1)}-v^{(N_2-1)}\Vert _{\mathcal {X}_v} \end{aligned} \end{aligned}$$\end{document}

for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C>0$$\end{document} . We then impose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C \Vert n_0\Vert _{L^\infty (\Omega )} < \lambda $$\end{document} , and choose \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda $$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} sufficiently small to conclude that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{(p^{(N)}, v^{(N)})\}_{N \ge 1}$$\end{document} is a Cauchy sequence. Thus, there exists (p, v) such that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{(p^{(N)}, v^{(N)})\}_{N \ge 1}$$\end{document} converges to it in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {X}_p\times \mathcal {X}_v$$\end{document} . Passing to the limit in (6) proves that (p, v) solves the pressure-volume problem.

(iv) Finally, to show that such a constructed solution is unique, we take

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ (p^{<1>}, v^{<1>}), \, (p^{<2>}, v^{<2>}) \in \mathbb {B}_{r_p, r_v}, $$\end{document}

and note that the difference \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(\overline{p}, \overline{v})=(p^{<1>}-p^{<2>}, v^{<1>}-v^{<2>})$$\end{document} solves

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} \ &\overline{p}_{tt} - c^2 \Delta \overline{p}- b \Delta \overline{p}_t = -\eta (\overline{p}(p^{<1>}+p^{<2>}))_{tt} - c^2 \rho _0 n_0 \overline{v}_{tt} \, & \text {in } \Omega \times (0,T),\\&\beta \overline{p}_t + \gamma \overline{p}+ \nabla \overline{p}\cdot \boldsymbol{n}= 0 & \text {on } \partial \Omega \times (0,T), \hspace{-0.28cm}\\ \end{aligned} \right. \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\overline{v}_{tt} + \delta \omega _0 \overline{v}_t + \omega _0^2 \overline{v}=- \mu \overline{p}+ \zeta \overline{v}\left( v^{<1>}+v^{<2>}\right) \\&\quad - \xi \left( 2 (\overline{v}v^{<1>}_{tt}+v^{<2>}\overline{v}_{tt} ) + \overline{v}_t(v^{<1>}_t+v^{<2>}_t)\right) . \end{aligned} \end{aligned}$$\end{document}

We obtain analogously to before

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert \overline{p}\Vert _{\mathcal {X}_p} \lesssim r_p\Vert \overline{p}\Vert _{\mathcal {X}_p}+ \Vert n_0\Vert _{L^\infty (\Omega )} \Vert \overline{v}\Vert _{\mathcal {X}_v} \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Vert \overline{v}\Vert _{\mathcal {X}_v} \lesssim \Vert \overline{p}\Vert _{\mathcal {X}_p}+ r_v\Vert \overline{v}\Vert _{\mathcal {X}_v}. \end{aligned} \end{aligned}$$\end{document}

For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert n_0\Vert _{L^\infty (\Omega )}$$\end{document} small enough, these inequalities allow us to conclude that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \overline{p}\Vert _{\mathcal {X}_p}=\Vert \overline{v}\Vert _{\mathcal {X}_v}=0$$\end{document} . This step completes the proof. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

The small-solution setting of Theorem 1 effectively forces the nonlinear ODE to retain the damped harmonic oscillator structure with T-periodic coefficients. Indeed, under the assumptions of Theorem 1, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} small enough, so that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \max \, \left\{ 2| \xi |, \frac{1}{\delta \omega _0}| \xi |, \frac{1}{\omega _0^2} |\zeta | \right\} \cdot r_v<1, \end{aligned}$$\end{document}

we can rewrite the ODE in the following form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} v_{tt}+ \frac{\delta \omega _0-\xi v_t}{1-2 \xi v} v_t+ \frac{\omega _0^2-\zeta v}{1-2 \xi v} v =- \frac{\mu }{1-2\xi v} p. \end{aligned} \end{aligned}$$\end{document}

Setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell = \dfrac{\delta \omega _0-\xi v_t}{1-2 \xi v}>0$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q =\dfrac{\omega _0^2-\zeta v}{1-2 \xi v}>0$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f = - \dfrac{\mu }{1-2\xi v} p$$\end{document} , we see that the volume fluctuation v can be represented as the solution of a second-order ODE (pointwise a.e. in space) with positive T-periodic coefficients and a T-periodic right-hand side, namely

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&v_{tt}+\ell (t)v_t+q(t) v =f(t), \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell (0)=\ell (T)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(0)=q(T)$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f(0)=f(T)$$\end{document} .

On the Smallness Assumption

The assumption (5) on the smallness of the microbubble number density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0$$\end{document} in Theorem 1 regulates the strength of the source of acoustic waves, along with the smallness of the external source h. To assess how meaningful this assumption is, we consider a nondimensional system obtained by introducing the scaling

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} x = L \tilde{x}, \qquad t = \frac{L}{c} \tilde{t}, \quad p = p^{ref }\, \tilde{p}, \qquad v = v^{ref }\, \tilde{v}, \end{aligned} \end{aligned}$$\end{document}

where L is a characteristic length scale, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{ref }$$\end{document} is a reference pressure amplitude, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^{ref }$$\end{document} a reference volume. After the change of variables and division by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{c^2 p^{ref }}{L^2}$$\end{document} , we obtain the dimensionless pressure equation:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned}&\ \ \tilde{p}_{tt} - \Delta \tilde{p} - \tilde{b} \Delta \tilde{p}_t = \tilde{\eta } (\tilde{p}^2)_{tt} + \tilde{\kappa } {n}_0 \tilde{v}_{tt} +\tilde{ \alpha } h \quad & \text {in } \Omega \times (0,T), \end{aligned} \end{aligned}$$\end{document}

where the relevant transformed coefficient is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tilde{\kappa } {n}_0 = \frac{ \rho _0 c^2 v^{ref }}{p^{ref }} {n}_0. \end{aligned}$$\end{document}

Since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^{ref }n_0$$\end{document} can be understood as the gas volume fraction (see [33, Sec. 2]), the smallness assumption corresponds to requiring moderate concentrations of microbubbles. For typical parameters in medical imaging, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c \sim 1500\hbox {m/s}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho \sim 1000\, \hbox {kg/m}^3$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^{ref }\in [10^5, 10^6]\, $$\end{document} Pa, and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0 \sim 10^{-6}\,$$\end{document} m, \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^{ref }= \frac{4}{3} R_0^3 \pi $$\end{document} , we obtain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\kappa } n_0 \ll 1$$\end{document} for microbubble densities in the range \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0 \in [10^9, 10^{12}]\,$$\end{document} bubbles \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$/m^3$$\end{document} .

Time Discretization Via a Multiharmonic Ansatz for Real Fields

In this section, we discuss the time discretization of the system by using a multiharmonic Ansatz, in the general spirit of [21]. We focus on a special case of having T-periodic acoustic excitation. That is, we assume that the acoustic source term has the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=g_{tt}$$\end{document} , where the function g is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g(x,t)=\mathfrak {R}\left\{ \sum _{m=0}^M \exp (\imath m \omega t)h_m^M(x)\right\} , \quad \omega = \dfrac{2 \pi }{T}, \quad h_m^M \in L^2(\Omega ; \mathbb {C}). \end{aligned}$$\end{document}

Under the assumptions of Theorem 1, a unique time-periodic solution (p, v) of the system given in (3) exists. For numerical approximations, we thus employ the following multiharmonic Ansatz:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u^N(x,t)=&\, \frac{1}{2} \sum _{m=0}^{N} \left( \exp (\imath m \omega t) u_m^N(x) + \exp (-\imath m \omega t) \overline{u^N_m(x)} \right) \\ =&\, \mathfrak {R}\left\{ \sum _{m=0}^{N} \exp (\imath m \omega t) u_m^N(x)\right\} . \end{aligned} \end{aligned}$$\end{document}

Alternatively, the multiharmonic Ansatz can be written in a more standard Fourier series form:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} u^N(x,t) = \sum _{m=0}^N \left[ u^c_m(x) \cos (m \omega t)+ u^s_m(x)\sin (m\omega t)\right] \end{aligned}$$\end{document}

with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^N_m(x) = u_m^c(x)- \imath u_m^s(x)$$\end{document} . We then look for the approximate pressure-volume field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N, v^N) \in X_N\times X_N$$\end{document} , with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} X_N = \left\{ \mathfrak {R}\left\{ \sum _{m=0}^N \exp ( \imath m \omega t) u_m^N(x) \right\} \, : \, u_m^N \in H^2(\Omega ; \mathbb {C}) \right\} , \end{aligned}$$\end{document}

where the coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N_m, v^N_m)$$\end{document} in the expansion are determined from the following system:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} \ &p_{tt}^N- c^2 \Delta p^N- b \Delta p_t^N & \\ =&\, \eta \,Proj _{X_N}\left[ ((p^N)^2)_{tt}\right] +c^2 \rho _0 n_0 v^N+ Proj _{X_N}h & \text {in } \Omega \times (0,T),\\&\beta p_t^N+ \gamma p^N+ \nabla p^N\cdot \boldsymbol{n}= 0 & \ \text {on } \partial \Omega \times (0,T),\\&v_{tt}^N+ \delta \omega _0 v_t^N+ \omega _0^2 v^N & \\ =&\, - \mu p^N+ \,Proj _{X_N}\left[ \zeta (v^N)^2 + \xi \left( 2 v^Nv_{tt}^N+ (v_t^N)^2 \right) \right] & \text {in } \Omega \times (0,T). \end{aligned} \right. \end{aligned}$$\end{document}

We next show that (13) can be equivalently rewritten as a system of Helmholtz problems and algebraic equations for computing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N_m, v^N_m)$$\end{document} .

A Multiharmonic Cut-Off Algorithm

Below we skip writing the dependencies of Fourier coefficients on x for readability.

Proposition 2

Let the assumptions of Theorem 1 hold with the acoustic source term assumed to have the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=g_{tt}$$\end{document} , where g is given in (11). The problem in (13) yields the following coupled system for computing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N_m, v^N_m)$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0 \le m \le N$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m = 0:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad p_0^N =0, \\ \text {(ii)} & \quad \omega _0^2 v_0^N = \zeta \left( v_0^N \right) ^2 + \sum _{j=1}^N \left( \frac{\zeta }{2} - \frac{\xi }{2} \omega ^2 j^2 \right) \left| v_j^N \right| ^2, \\ \end{aligned} \end{array}\right. } \\ m= 1:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_1^N - \frac{c^2 + \imath b \omega }{\omega ^2} \Delta p_1^N+ c^2 n_0\rho _0 v_1^N = - h_1^N - \sum _{k=3:2}^{2N-1} \eta \overline{p_{\frac{k-1}{2}}^N} p_{\frac{k+1}{2}}^N, \\ \text {(ii)} & \quad \left( - \omega ^2 + \imath \delta \omega _0 \omega + \omega _0^2 \right) v_1^N + \mu p_1^N\\ & \quad = (\zeta -\xi \omega ^2) v_0^N v_1^N + \sum _{k=1:2}^{2N-1} \left( \zeta - \xi \omega ^2 \frac{k^2 + 3}{4}\right) \overline{v_{\frac{k-1}{2}}^N} v_{\frac{k+1}{2}}^N, \\ \end{aligned} \end{array}\right. } \\ m = \left\{ 2, \ldots , N \right\} :&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_m^N- \frac{ c^2 + \imath m b \omega }{m^2 \omega ^2} \Delta p_m^N + c^2 \rho _0 n_0 v_m^N\\ & \quad = - h_m^N -\eta \left( \sum _{l=1}^{m-1} p_l^N p_{m-l}^N- 2 \sum _{k=m+2:2}^{2N-m} \overline{p_{\frac{k-m}{2}}^N} p_{\frac{k+m}{2}}^N \right) , \\ \text {(ii)}& \quad \left( - \omega ^2m^2 + \imath \delta \omega _0 \omega m+\omega _0^2\right) v_m^N + \mu p_m^N, \\ & \quad = \sum _{l=0}^m \left( \frac{\zeta }{2} - \frac{\xi \omega ^2}{2} (m-l)(2m-l) \right) v_l^N v_{m-l}^N \\ & \quad \quad + \sum _{k=m:2}^{2N-m} \left( \zeta - \xi \omega ^2 \frac{k^2 + 3m^2}{4}\right) \overline{v_{\frac{k-m}{2}}^N} v_{\frac{k+m}{2}}^N. \\ \end{aligned} \end{array}\right. } \end{aligned}$$\end{document}

Additionally, each of the individual functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_m^N$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 1, \ldots , N \right\} $$\end{document} satisfies the following boundary conditions:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} (\imath \omega m \beta + \gamma ) p_m^N + \nabla p_m^N \cdot \boldsymbol{n}= 0 \qquad \text { on }\, \partial \Omega . \end{aligned}$$\end{document}

Proof

We begin by considering the derivation of the multiharmonic algorithm for the Westervelt’s equation and the ODE separately. For the PDE, one can follow the steps outlined in [21], where the de-coupled Westervelt equation is considered, by setting the right-hand side to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$g = h + c^2 \rho _0 n_0 v$$\end{document} . Projections onto \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_N$$\end{document} of product terms appearing on the right-hand side of (13) can be expressed using [21, Lemma 3.1].

For the ODE in (13), the Ansatz given in (12), together with [21, Lemma 3.1], yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 0 =&\mathfrak {R}\left\{ - \omega ^2 \sum _{m=0}^N v_m^N m^2 \exp (\imath m \omega t) + \delta \omega _0 \imath \omega \sum _{m=0}^N v_m^Nm \exp ( \imath m \omega t) + \omega _0^2 \sum _{m=0}^N v_m^N \exp (\imath m \omega t) \right. \\&\left. + \mu \sum _{m=0}^N p_m^N \exp (\imath m \omega t)- \frac{\zeta }{2} \left( (v_0^N)^2 + \sum _{j=0}^N \left| v_j^N\right| ^2 \right. \right. \\&\left. \left. + \sum _{m=1}^N \left[ \sum _{l=0}^m v_l^N v_{m-l}^N + 2 \sum _{k=m:2}^{2N-m} \overline{v_{\frac{k-m}{2}}^N} v_{\frac{k+m}{2}}^N \right] \exp (\imath m \omega t) \right) + \xi \frac{\omega ^2}{2}\left( \sum _{j=0}^N j^2 \left| v_j^N\right| ^2 \right. \right. \\&\left. \left. + \sum _{m=1}^N \left[ \sum _{l=0}^m (m-l)(2m-l) v_l^N v_{m-l}^N + \sum _{k=m:2}^{2N-m} \frac{k^2 +3m^2}{2} \overline{v_{\frac{k-m}{2}}^N} v_{\frac{k+m}{2}}^N \right] \exp (\imath m \omega t) \right) \right\} . \end{aligned}$$\end{document}

From here, we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} 0 =&\mathfrak {R}\left\{ \sum _{m=0}^N \left[ \left( - \omega ^2 m^2 + \delta \omega _0 \imath \omega m + \omega _0^2 \right) v_m^N + \mu p_m^N \right] \exp (\imath m \omega t) \right. \\&\left. - \frac{\zeta }{2} (v_0^N)^2 + \sum _{j=0}^N \left( - \frac{\zeta }{2} + \xi \frac{\omega ^2}{2} j^2 \right) \left| v_j^N\right| ^2 + \sum _{m=1}^N \left[ \sum _{l=0}^m \left( -\frac{\zeta }{2} + \xi \frac{\omega ^2}{2} (m-l)(2m-l) \right) v_l^N v_{m-l}^N \right. \right. \\&\left. \left. + \sum _{k=m:2}^{2N-m} \left( -\zeta + \xi \omega ^2 \frac{k^2 + 3m^2}{4}\right) \overline{v_{\frac{k-m}{2}}^N} v_{\frac{k+m}{2}}^N \right] \exp (\imath m \omega t) \right\} . \end{aligned}$$\end{document}

Combining our derivations above and using linear independence of the functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ t \mapsto \exp ( \imath m \omega t)$$\end{document} yields a coupled nonlinear system for the functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\left\{ p_0^N, v_0^N, \cdots , p_N^N, v_N^N \right\} $$\end{document} . Additionally, each of the individual functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_m^N$$\end{document} must satisfy the boundary conditions in (14). The 0-th equation is given by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Delta p^N_0 = 0$$\end{document} , which together with the 0-th boundary condition (and the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma >0$$\end{document} ) implies that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_0^N$$\end{document} vanishes. By dividing the resulting PDEs by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m^2\omega ^2$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \ge 1$$\end{document} , we arrive at the claimed coupled system. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

We observe that without the sums over k on the right-hand side of the system in Proposition 2, the system would be triangular and could be solved by substitution. One way of obtaining a triangular form is by linearizing the right-hand side as we show next.

A Linearized Multiharmonic Cut-Off Algorithm

Next, we propose a linearized multilevel method. That is, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , we consider the sequence of the following linearized equations where the quadratic terms are approximated by taking into account \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N-1$$\end{document} harmonics:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned} \ &p_{tt}^N- c^2 \Delta p^N- b \Delta p_t^N & \\ =&\, \eta \,Proj _{X_N}\left[ ((p^{N-1})^2)_{tt}\right] +c^2 \rho _0 n_0 v_{tt}^N+ Proj _{X_N}h & \text {in } \Omega \times (0,T),\\&\beta p_t^N+ \gamma p^N+ \nabla p^N\cdot \boldsymbol{n}= 0 & \text {on } \partial \Omega \times (0,T),\\&v_{tt}^N+ \delta \omega _0 v_t^N+ \omega _0^2 v^N & \\ =&\, - \mu p^N+ \,Proj _{X_N}\left[ \zeta (v^{N-1})^2 + \xi \left( 2 v^{N-1}v_{tt}^{N-1}+ (v_t^{N-1})^2 \right) \right] & \text {in } \Omega \times (0,T), \end{aligned} \right. \end{aligned}$$\end{document}

where for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^N \in X_N$$\end{document} , we formally define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_m^N = 0$$\end{document} for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m>N$$\end{document} . To start the algorithm, we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^0 =0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^0 = 0$$\end{document} . We next show that this problem can be seen as a triangular system of Helmholtz problems and algebraic equations that can be solved by successive substitution.

Proposition 3

Let the assumptions of Theorem 1 hold with the acoustic source term assumed to have the form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=g_{tt}$$\end{document} , where g is given in (11). Then, an equivalent formulation of the problem in (15) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^0=v^0=0$$\end{document} is given by the following coupled system:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m = 0:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad p_0^N =0, \\ \text {(ii)} & \quad \omega _0^2 v_0^N= \frac{\zeta }{2} \left( v_0^{N-1} \right) ^2 + \sum _{j=0}^{N-1} \left( \frac{\zeta }{2} - \frac{\xi \omega ^2}{2}j^2 \right) \left| v_j^{N-1} \right| ^2,\\ \end{aligned} \end{array}\right. } \\ m= 1:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_1^N - \frac{ c^2 + \imath b \omega }{\omega ^2} \Delta p_1^N+ c^2 \rho _0 n_0 v_1^{N} = -h_1^N - \eta \sum _{k=3:2}^{2N-3} \overline{p_{\frac{k-1}{2}}^{N-1}} p_{\frac{k+1}{2}}^{N-1}, \\ \text {(ii)} & \quad \frac{1}{\alpha _1}v_1^N + \mu p_1^{N} = \left( \zeta - \xi \omega ^2\right) v_0^{N-1}v_1^{N-1} \\ & \hspace{2.5cm} + \sum _{k=1:2}^{2N-3} \left( \zeta - \xi \omega ^2 \frac{k^2 + 3}{4}\right) \overline{v_{\frac{k-1}{2}}^{N-1}} v_{\frac{k+1}{2}}^{N-1}, \\ \end{aligned} \end{array}\right. } \\ m = \left\{ 2, \ldots , N \right\} :&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_m^N- \frac{\left( c^2 + \imath m b \omega \right) }{ \omega ^2 m^2}\Delta p_m^N + c^2 \rho _0 n_0 v_m^{N}\\ & \quad = - h_m^N - \frac{ \eta }{2} \left( \sum _{l=1}^{m-1} p_l^{N-1} p_{m-l}^{N-1} -2 \sum _{k = m+2:2}^{2(N-1)-m} \overline{p_{\frac{k-m}{2}}^{N-1}} p_{\frac{k+m}{2}}^{N-1} \right) , \\ \text {(ii)}& \quad \frac{1}{\alpha _m} v_m^N + \mu p_m^{N} = \sum _{l=0}^m \left( \frac{\zeta }{2} - \frac{\xi \omega ^2}{2} (m-l)(2m-l) \right) v_l^{N-1} v_{m-l}^{N-1} \\ & \hspace{2.5cm} + \sum _{k=m:2}^{2(N-1)-m} \left( \zeta - \frac{\xi \omega ^2 (k^2 + 3 m^2)}{4} \right) \overline{v_{\frac{k-m}{2}}^{N-1}} v_{\frac{k+m}{2}}^{N-1} \end{aligned} \end{array}\right. } \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , where for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \ge 1$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _m = \left( - m^2\omega ^2 + \imath m \delta \omega _0 \omega +\omega _0^2\right) ^{-1}$$\end{document} . Additionally, each of the individual functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_m^N$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 1, \ldots , N \right\} $$\end{document} , should satisfy boundary conditions (14).

Proof

The statement follows analogously to before by applying the Ansatz given in (12) and making use of the expressions provided in [21, Lemma 3.1] for projections onto \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_N$$\end{document} of products of two functions in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_N$$\end{document} . \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Looking at the system in Proposition 3, we see that we can express the volume harmonics as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} v^N_m = -\alpha _m \mu p^N_m + \alpha _m f_v^{N-1} \text { for } \ m \ge 1, \end{aligned}$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f^{N-1}_v&= \sum _{l=0}^m \left( \frac{\zeta }{2} - \frac{\xi \omega ^2}{2} (m-l)(2m-l) \right) v_l^{N-1} v_{m-l}^{N-1}\\&\quad + \sum _{k=m:2}^{2(N-1)-m} \left( \zeta - \frac{\xi \omega ^2 (k^2 + 3 m^2)}{4} \right) \overline{v_{\frac{k-m}{2}}^{N-1}} v_{\frac{k+m}{2}}^{N-1}. \end{aligned}$$\end{document}

By substituting these values into the equations for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N_m$$\end{document} , the coupled system in Proposition 3 can be rewritten as a system of (16) and Helmholtz equations given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned}&-\left( 1+\imath \omega \frac{m b}{c^2} \right) \Delta p^N_m - (k^2+\mathfrak {a}) p^N_m +\imath \mathfrak {b}p^N_m = -k^2 h_m^N - f_{p,v}^{N-1} \quad \text { in }\, \Omega ,\quad \text {with } k= \frac{m \omega }{c},\\&(\imath \omega m \beta + \gamma ) p^N_m + \nabla p^N_m \cdot \boldsymbol{n}= 0 \quad \text { on }\, \partial \Omega , \end{aligned} \right. \end{aligned}$$\end{document}

where the source is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} f_{p,v}^{N-1} = \alpha _m \omega ^2 m^2 \rho _0 n_0 f^{N-1}_v + k^2 \frac{ \eta }{2} f^{N-1}_p \end{aligned}$$\end{document}

with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ f^{N-1}_p = \sum _{l=1}^{m-1} p_l^{N-1} p_{m-l}^{N-1} -2 \sum _{k = m+2:2}^{2(N-1)-m} \overline{p_{\frac{k-m}{2}}^{N-1}} p_{\frac{k+m}{2}}^{N-1}. $$\end{document}

In (17), we have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathfrak {a}= \mu \rho _0 n_0 \frac{m^2 \omega ^2 (\omega _0^2 - m^2 \omega ^2)}{(\omega _0^2-m^2\omega ^2)^2 + (m \delta \omega _0 \omega )^2}\ \text { and } \ \mathfrak {b}= \mu \rho _0 n_0 \frac{m^3 \omega ^3 \delta \omega _0 }{(\omega _0^2-m^2\omega ^2)^2 + (m \delta \omega _0 \omega )^2}. \end{aligned}$$\end{document}

Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {b}>0$$\end{document} , while the sign of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {a}$$\end{document} depends on the relation between \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega _0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \omega $$\end{document} . The form of the Helmholtz equation in (17) reveals more clearly the influence of microbubbles on the wave propagation, through the dissipation signaled by the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$i \mathfrak {b}p^N_m$$\end{document} term and modification of the wave number via the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-(k^2+\mathfrak {a}) p^N_m$$\end{document} term. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Convergence of the Linearized Multiharmonic Algorithm for Real Fields

We next wish to establish convergence of the iterative scheme in (15) as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \rightarrow \infty $$\end{document} in a suitable norm. We can mimic the proof of Theorem 1 to show by induction that for each \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , problem (15) has a unique solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N, v^N) \in X_N\times X_N$$\end{document} . For fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , the existence and uniqueness of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N$$\end{document} follow by noting that, thanks to Proposition 3, harmonics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N_m$$\end{document} are obtained as solutions of the Helmholtz problems in (17). Furthermore, an analogous testing procedure to the one in Theorem 1 applied on (15) yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \Vert p^N\Vert _{\mathcal {X}_p} \le r_p, \quad \Vert v^N\Vert _{\mathcal {X}_v} \le r_v, $$\end{document}

provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \mu \Vert _{L^\infty (\Omega )}$$\end{document} are sufficiently small. We omit these details here and focus on establishing convergence. The convergence will be shown in the norms of the spaces \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y^0_p$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Y^0_{v}$$\end{document} , where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} Y_p^\ell =&\, H^{\ell +1}(0,T; H^1(\Omega )) \cap H^{\ell +2}(0,T; L^2(\Omega )) \cap H^{\ell +2}(0,T; L^2(\partial \Omega )),\\ Y_{v}^\ell =&\, H^{\ell +2}(0,T; L^2(\Omega )) \end{aligned} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell \ge 0$$\end{document} . To make writing more compact, we introduce the short-hand norm notation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} |\!|\!|(p,v)|\!|\!|_{Y^\ell _p\times Y^\ell _{v}} :=\Vert p\Vert _{Y^\ell _p}+\Vert v\Vert _{Y^\ell _{v}} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p,v) \in Y^\ell _p\times Y^\ell _{v}$$\end{document} , and we also adopt the operator notation

$$[eqn]$$

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} <\mathcal {L}_2 v^N, \phi ^N> :=\int _0^T \int _{\Omega }\left( v_{tt}^N+ \delta \omega _0 v_t^N+ \omega _0^2 v^N\right) \phi ^N\, d xd s, \quad \phi ^N\in X_N. \end{aligned}$$\end{document}

Then the algorithm can be restated as follows. Given the previous iterate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^{N-1}, v^{N-1}) \in X_N\times X_N$$\end{document} , we compute \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N, v^N)$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} as the solution of

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} <\mathcal {L}_1 p^N, \phi ^N> =&\, \int _0^T \int _{\Omega }\left( {\eta } ((p^{N-1})^2)_{tt} + c^2 \rho _0 n_0 v_{tt}^N+h\right) \phi ^N\, d xd s, \end{aligned} \end{aligned}$$\end{document}

and

$$[eqn]$$

for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi ^N\in X_N$$\end{document} .

In the convergence proof, we will involve the truncated Fourier series of a continuous-in-time function u given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \tilde{u}^N(x,t)= \sum _{m=0}^N \left[ u^{c}_m (x) \cos (m \omega t)+ u_m^{s}(x)\sin (m \omega t)\right] \in X_N, \quad N \ge 0, \end{aligned} \end{aligned}$$\end{document}

with the coefficients

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} u_m^{c}(x) = \frac{2}{T} \int _0^T u(x,t) \cos (m \omega t) \, d t, \quad u_m^{s}(x) = \frac{2}{T} \int _0^T u(x,t) \sin (m \omega t) \, d t. \end{aligned} \end{aligned}$$\end{document}

By integrating the time-periodic Westervelt equation and acoustic boundary conditions in (3) over (0, T), we can conclude that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _0^T p(x,t) \, d t=0$$\end{document} in \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega $$\end{document} . Thus \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{p}^0=0$$\end{document} , which is reflected by having \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N_0=0$$\end{document} in the multiharmonic algorithms.

It can be shown analogously to [3, Lemma 12] that the following error bound holds for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\ell \ge 1$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} |\!|\!|(p-\tilde{p}^N, v-\tilde{v}^N)|\!|\!|_{(Y^0_p\cap H^2(H^1(\Omega ))) \times Y^0_{v}} \le&\, C N^{-\ell } (\Vert p\Vert _{Y_p^{\ell } \cap H^{\ell +2}(H^1(\Omega ))}+\Vert v\Vert _{Y_v^{\ell }}), \end{aligned} \end{aligned}$$\end{document}

which we can exploit to characterize the error of the scheme in the next statement. We include the proof of (19) in Appendix B for completeness.

Theorem 2

Let the assumptions of Theorem 1 hold with the acoustic source term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=g_{tt}$$\end{document} , where g is given in (11). Assume additionally that

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} (p,v) \in \left( \mathcal {X}_p\cap Y_p^\ell \cap H^{\ell +2}(0,T; H^1(\Omega )) \right) \times \left( \mathcal {X}_v\cap Y_{v}^\ell \right) , \quad \ell \ge 1. \end{aligned} \end{aligned}$$\end{document}

Let \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^0=v^0=0$$\end{document} and assume that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N, v^N)$$\end{document} are computed using (18) for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} . Then, provided \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \mu \Vert _{L^\infty (\Omega )}$$\end{document} are sufficiently small, there exists \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q=q(r_p, r_v)<1$$\end{document} , such that

$$[eqn]$$

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C>0$$\end{document} does not depend on N.

Proof

The proof follows by splitting the error as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} p -p^N=&\, (p- \tilde{p}^N)- (p^N- \tilde{p}^N) {=}{:}err (\tilde{p}^N)-e^{p,N}, \\ v -v^N=&\, (v-\tilde{v}^N) -(v^N-\tilde{v}^N) {=}{:}err (\tilde{v}^N)-e^{v,N}, \end{aligned} \end{aligned}$$\end{document}

and then representing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(e^{p,N}, e^{v,N})$$\end{document} as the solution of a suitable semi-discrete wave-volume system. The discrete error \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(e^{p,N}, e^{v,N})$$\end{document} can be seen as the solution of

$$[eqn]$$

and

$$[eqn]$$

for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\phi ^N\in X_N$$\end{document} . Note that in the testing procedure for this problem we can exploit the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _0^T \frac{d }{d t}(\cdot ) \, d t=0$$\end{document} for time-periodic functions. By testing (21a) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^{p,N}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_t^{p,N}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_{tt}^{p,N}$$\end{document} , employing Young’s and Hölder’s inequalities, and combining the resulting estimates, we can derive the following bound:

$$[eqn]$$

where we have also exploited the equivalence of the norms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert w\Vert _{H^1(\Omega )}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert \nabla w\Vert _{L^2(\Omega )}+\Vert w\Vert _{L^2(\partial \Omega )}$$\end{document} (see [32, Theorem 1.9]). The derivation of (22) is provided in Appendix C.

Then by using the rewriting

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} ((p^{N-1})^2)_{tt}-(p^2)_{tt}&=\, 2(p^{N-1}_t-p_t)(p^{N-1}_t+p_t)+2(p^{N-1}-p)p^{N-1}_{tt}\\&\quad +2 p(p^{N-1}_{tt}-p_{tt}), \end{aligned} \end{aligned}$$\end{document}

and the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert p\Vert _{\mathcal {X}_p}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert p^N\Vert \le r_p$$\end{document} for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , we obtain

$$[eqn]$$

From here via \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p-p^N=err (\tilde{p}^N)-e^{p,N}$$\end{document} and the triangle inequality, we arrive at

$$[eqn]$$

Similarly, testing (21b) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e^{v,N}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_t^{v,N}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$e_{tt}^{v,N}$$\end{document} yields, after standard manipulations,

$$[eqn]$$

By using the fact that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert v\Vert _{\mathcal {X}_v}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert v^N\Vert _{\mathcal {X}_v} \le r_v$$\end{document} for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 1$$\end{document} , and the triangle inequality, from (21b) we then have

$$[eqn]$$

Adding \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda \cdot $$\end{document} (24) and (23) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda >0$$\end{document} small enough, so that the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\lambda \Vert p^N-p\Vert _{L^2(L^2(\Omega ))}$$\end{document} can be absorbed by the left-hand side, and then possibly reducing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert n_0\Vert _{L^\infty (\Omega )}$$\end{document} so that the term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert n_0\Vert _{L^\infty (\Omega )} \Vert v_{tt}^N-v_{tt}\Vert _{L^2(L^2(\Omega ))}$$\end{document} can be absorbed, yields

$$[eqn]$$

By possibly further reducing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_p$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_v$$\end{document} , we obtain \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$q(r_p, r_v) <1$$\end{document} , such that

$$[eqn]$$

for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C_0>0$$\end{document} . The statement then follows by employing (19). \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Convergence of the scheme then follows in a straightforward manner from (22).

Corollary 1

(Convergence of the linearized scheme) Let the assumptions of Theorem 2 hold. Then the solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N, v^N) \in X_N\times X_N$$\end{document} of the iteration scheme (18) converges with respect to the norm \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|\!|\!|(\cdot , \cdot )|\!|\!|_{Y^0_p\times Y^0_{v}}$$\end{document} to the solution \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p, v) \in \mathbb {B}_{r_p, r_v}$$\end{document} of (3) as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \rightarrow \infty $$\end{document} .

Proof

From (20), by iteration we obtain

$$[eqn]$$

for some \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$C>0$$\end{document} , independent of N. The statement then follows analogously to [34, Theorem 3.1], where the multiharmonic discretization of the de-coupled Westervelt equation is studied. We thus omit the details here. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

Setting the Zeroth Microbubble Harmonic to Zero

Unlike with the Westervelt equation, from the equation for v, we cannot in general conclude by integrating from 0 to T that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\int _0^T v(t) \, d t=0$$\end{document} . Nevertheless, in numerical methods for solving multiharmonic problems it is often assumed that the zeroth harmonic does not contribute significantly to the dynamics; this assumption is also made in the formal two-harmonic generation study for the linearized wave-ODE model with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$b=\eta =0$$\end{document} in [17, Ch. 5.3.2].

With the assumption \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_0^N = v_0^N = 0$$\end{document} for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N \ge 0$$\end{document} , the system derived in Proposition 3 further simplifies to

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m= 1:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_1^{N} - \frac{ c^2 + \imath b \omega }{\omega ^2} \Delta p_1^N+ c^2 \rho _0 n_0 v_1^{N} = -h_1^N - \sum _{k=3:2}^{2N-3} \eta \overline{p_{\frac{k-1}{2}}^{N-1}} p_{\frac{k+1}{2}}^{N-1} \\ \text {(ii)} & \quad v_1^N = \alpha _1 \left( - \mu p_1^{N} + \sum _{k=3:2}^{2N-3} \left( \zeta - \xi \omega ^2 \frac{k^2 + 3}{4}\right) \overline{v_{\frac{k-1}{2}}^{N-1}} v_{\frac{k+1}{2}}^{N-1}\right) \\ \end{aligned} \end{array}\right. } \end{aligned}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m = \left\{ 2, \ldots , N \right\} :&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_m^N- \frac{\left( c^2 + \imath m b \omega \right) }{ \omega ^2 m^2}\Delta p_m^N + c^2 \rho _0 n_0 v_m^{N}\\ & \quad = - h_m^N - \frac{ \eta }{2} \left( \sum _{l=1}^{m-1} p_l^{N-1} p_{m-l}^{N-1} -2 \sum _{k = m+2:2}^{2(N-1)-m} \overline{p_{\frac{k-m}{2}}^{N-1}} p_{\frac{k+m}{2}}^{N-1} \right) \\ \text {(ii)}& \quad v_m^N = \alpha _m \left( -\mu p_m^{N} + \sum _{l=1}^{m-1} \left( \frac{\zeta }{2} - \frac{\xi \omega ^2}{2} (m-l)(2m-l) \right) v_l^{N-1} v_{m-l}^{N-1} \right. \\ & \quad \left. + \sum _{k=m+2:2}^{2(N-1)-m} \left( \zeta {-} \frac{\xi \omega ^2 (k^2 + 3 m^2)}{4} \right) \overline{v_{\frac{k-m}{2}}^{N-1}} v_{\frac{k+m}{2}}^{N-1} \right) \\ \end{aligned} \end{array}\right. } \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _m = \left( - m^2 \omega ^2 + \imath m \delta \omega _0 \omega + \omega _0^2 \right) ^{-1}$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 1, \ldots N\right\} $$\end{document} . Now \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_m^N$$\end{document} can be eliminated explicitly by substituting the second equation into the first, leading to a system of inhomogeneous Helmholtz equations.

Time Discretization Via a Multiharmonic Ansatz for Complex Fields

In this section, we formally investigate multiharmonic algorithms based on neglecting the fact that the excitation and solution fields should be real-valued. We look for approximate solutions in

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \tilde{X}_N := \left\{ \sum _{m=0}^N \exp ( \imath m \omega t) u_m^N(x) \, : \, u_m^N \in H^2(\Omega ; \mathbb {C}) \right\} . $$\end{document}

This approach, adopted also in the investigation of the (de-coupled) Westervelt equation in [21], will allow us to arrive at a simple triangular approximation scheme in the frequency domain. To this end, we simply set

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} g^M(x,t) = \displaystyle \sum _{m=0}^{M}\exp (\imath \omega m t) h_m^M(x), \quad \omega = \frac{2 \pi }{T}, \end{aligned}$$\end{document}

and make the following Ansatz:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} p^N(x,t)= \sum _{m=0}^{N} \exp (\imath m \omega t) p_m^N(x), \qquad v^N(x,t)= \sum _{m=0}^{N} \exp (\imath m \omega t) v_m^N(x), \end{aligned}$$\end{document}

where the coefficients \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N_m$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v^N_m \in H^2(\Omega ; \mathbb {C})$$\end{document} are determined by solving (13) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{X_N}(\cdot )$$\end{document} replaced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{\tilde{X}_N}(\cdot )$$\end{document} .

A Multiharmonic Cut-Off Algorithm

We next set up a simplified algorithm (compared to the one in Proposition 2) in the sense that the right-hand side in the resulting system will contain only one summation term, such that there is reduced coupling across harmonics, making it computational more efficient.

In this complex setting, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0^N$$\end{document} , we obtain the quadratic equation

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \omega _0^2 v_0^N - \zeta (v_0^N)^2 = 0 \end{aligned}$$\end{document}

and choose the zero solution, such that we have \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0^N=0$$\end{document} as before in (25). This choice is consistent with the formulation used throughout the numerical simulations and avoids introducing spurious constant components.

Proposition 4

Let the assumptions of Theorem 1 hold with the acoustic source term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=g_{tt}$$\end{document} , where g is given in (26). Under the Ansatz in (27) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N_0=v^N_0=0$$\end{document} , the problem in (13) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{X_N}(\cdot )$$\end{document} replaced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{\tilde{X}_N}(\cdot )$$\end{document} can be rewritten as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} m= 1:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_1^N- \frac{c^2 + \imath \omega b }{\omega ^2} \Delta p^N_1 + c^2 \rho _0 n_0 v_1^N = -h_1^N \\ \text {(ii)} & \quad \frac{1}{\alpha _1} v_1^N + \mu p_1^N =0\\ \end{aligned} \end{array}\right. } \\ m = \left\{ 2, \ldots , N \right\} :&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_m^N - \frac{(c^2 + \imath m \omega b)}{m^2 \omega ^2} \Delta p_m^N + c^2 \rho _0 n_0 v_m^N \\ & \qquad =- h_m^N- \eta \sum _{l=1}^{m-1} p_l^N p_{m-l}^N \\ \text {(ii)} & \quad \frac{1}{\alpha _m} v_m^N + \mu p_m^N\\ & \quad = \sum _{l=1}^{m-1} \left( \zeta - \xi \omega ^2 (m-l)(2m-l) \right) v_{l}^Nv^{N}_{m-l}, \end{aligned} \end{array}\right. } \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _m = \left( -m^2\omega ^2 + \imath m \delta \omega _0 \omega + \omega _0^2\right) ^{-1}$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 1, \ldots N\right\} $$\end{document} . Additionally, each of the individual functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_m^N$$\end{document} , where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 0, \ldots , N \right\} $$\end{document} , should satisfy boundary conditions (14).

Proof

Plugging the approximation given in (27) into (3) yields

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sum _{m=0}^{N}&\left[ - m^2 \omega ^2 p_m^N- (c^2 + \imath m \omega b) \Delta p_m^N + c^2 \rho _0 n_0 m^2 \omega ^2 v_m^N + m^2\omega ^2 h_m^N\right. \\&\left. \qquad \qquad + \sum _{l=0}^m \eta m^2 \omega ^2 p_l^N p_{m-l}^N \right] \exp (\imath m \omega t)= 0. \end{aligned}$$\end{document}

Concerning the nonlinear term on the right hand-side of the ODE, we obtain

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \text {Proj}_{\tilde{X}_N}(2 v^Nv^N_{tt} + (v^N_t)^2)&= - \sum _{m=0}^N \sum _{l=0}^m \left( 2(m-l)^2 + l (m-l) \right) \omega ^2 v_l^N v_{m-l}^N \exp (\imath m \omega t) \\&= - \sum _{m=0}^N \sum _{l=0}^m (m-l)(2m-l) \omega ^2 v_l^N v_{m-l}^N \exp (\imath m \omega t). \end{aligned}$$\end{document}

Therefore, the approximation for the ODE is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \sum _{m=0}^{N}&\left[ \phantom {\sum _{l=0}^m} \left( -m^2 \omega ^2 + \delta \omega _0 \imath m \omega + \omega _0^2 \right) v_m^N + \mu p_m^N \right. \\&\left. + \sum _{l=0}^m \left( - \zeta + \xi \omega ^2 (m-l)(2m-l)\right) v_l^N v_{m-l}^N \right] \exp (\imath m \omega t) = 0. \end{aligned}$$\end{document}

Using linear independence of the functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t \mapsto \exp (\imath m \omega t)$$\end{document} , we immediately arrive at the iterative coupled system. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

A Linearized Multiharmonic Cut-Off Algorithm

We next also want to look at the linearized set-up given in (15) (with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{X_N}(\cdot )$$\end{document} replaced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{\tilde{X}_N}(\cdot )$$\end{document} ) in the setting of complex solution fields. In this setting, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u^N \in \tilde{X}_N$$\end{document} , we formally define \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$u_m^N = 0$$\end{document} for all \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m>N$$\end{document} .

Proposition 5

Let the assumptions of Theorem 1 hold with the acoustic source term \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=g_{tt}$$\end{document} , where g is given in (26). Under the Ansatz in (27) with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p^N_0=v^N_0=0$$\end{document} , the problem in (15), with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{X_N}(\cdot )$$\end{document} replaced by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$Proj _{\tilde{X}_N}(\cdot )$$\end{document} , can be rewritten as follows:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} m= 1:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_1^N- \frac{c^2 + \imath \omega b }{\omega ^2} \Delta p^N_1 + c^2 \rho _0 n_0 v_1^{N} = -h_1^N, \\ \text {(ii)} & \quad \frac{1}{\alpha _1} v_1^N + \mu p_1^{N} =0, \\ \end{aligned} \end{array}\right. } \\ m = \left\{ 2, \ldots , N \right\} :&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad - p_m^N - \frac{(c^2 + \imath m \omega b)}{m^2 \omega ^2} \Delta p_m^N + c^2 \rho _0 n_0 v_m^{N} \\ & \quad =- h_m^N-\eta \sum _{l=1}^{m-1} p_l^{N-1} p_{m-l}^{N-1}, \\ \text {(ii)} & \quad \frac{1}{\alpha _m} v_m^N + \mu p_m^{N} \\ & \quad = \sum _{l=1}^{m-1} \left( \zeta - \xi \omega ^2 (m-l)(2m-l) \right) v_{l}^{N-1}v^{N-1}_{m-l}, \end{aligned} \end{array}\right. } \end{aligned} \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _m = \left( -m^2\omega ^2 + \imath m \delta \omega _0 \omega + \omega _0^2\right) ^{-1}$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 1, \ldots N\right\} $$\end{document} . Additionally, each of the individual functions \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_m^N$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \left\{ 0, \ldots , N \right\} $$\end{document} is supposed to satisfy boundary conditions (14).

Proof

Similarly to the previous case, substituting the approximation provided in (27) into (3) immediately yields the iterative coupled system. \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\square $$\end{document}

We see that when working with complex fields from the beginning, we end up with a lower triangular system of inhomogeneous Helmholtz equations and algebraic equations that can be solved by substitution.

Remark 1

(On the assumption of complex fields) In this section, the approximations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$(p^N,v^N)$$\end{document} are complex-valued, so convergence could only be studied in the complex counterpart of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$X_N\times X_N$$\end{document} . This does not by itself imply that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Re (p^N,v^N)\rightarrow (p,v)$$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N\rightarrow \infty $$\end{document} ; such a conclusion would require enforcing the reality constraint (symmetric \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm m$$\end{document} spectrum with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}_{-m}=\overline{\hat{u}_m}$$\end{document} ). Following [3], one may formulate the multiharmonic Ansatz with complex Fourier coefficients (omitting the explicit \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Re $$\end{document} in the time reconstruction), but the resulting nonlinear maps are not complex-differentiable (holomorphic). Hence a Newton scheme in the sense of complex analysis is not applicable; in practice one uses the real formulation (or a split into real/imaginary parts) despite the algebraic brevity of the complex notation.

For clarity, eq. (25) (real basis) contains both sum- and difference-frequency couplings generated by the quadratic terms, visible in the two convolution sums \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sum _{l=1}^{m-1} p_l\,p_{m-l}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \sum _{k=m+2:2}^{2(N-1)-m} p_{(k-m)/2}\,p_{(k+m)/2}$$\end{document} . In contrast, the complex projection (28) is posed on the half-spectrum \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\{e^{\mathrm i m\omega _0 t}\}_{m=1}^N$$\end{document} without imposing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\hat{u}_{-m}=\overline{\hat{u}_m}$$\end{document} ; only the sum-frequency interactions remain and the block system becomes lower-triangular (solvable by substitution). If one augments the complex space to include \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\pm m$$\end{document} and enforces the reality constraint, the complex and real formulations are algebraically equivalent. In our parameter regime, the neglected difference-frequency contributions are numerically negligible, which explains the agreement observed in Fig. 3 in Section 6.

A Simple Two-Harmonic Approximation Scheme

In certain imaging applications, the fundamental frequency and the second harmonic are considered the most significant; see, for example, [31]. We thus compare the multiharmonic iterative scheme with a simplified two-harmonic scheme to evaluate the extent of information loss.

For two harmonics, that is, for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_0^{(2)}=0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0^{(2)}=0$$\end{document} , the approximations for the pressure p and the volume v are given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} p(x,t) \approx p^{(2)}(x,t)&= p_1^{(2)}(x) \exp (\imath \omega t) + p_2^{(2)}(x) \exp (2 \imath \omega t), \\ v(x,t) \approx v^{(2)}(x,t)&= v_1^{(2)}(x) \exp (\imath \omega t) + v_2^{(2)}(x) \exp (2 \imath \omega t). \end{aligned} \end{aligned}$$\end{document}

The corresponding scheme for computing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{1,2}^{(2)}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_{1,2}^{(2)}$$\end{document} can be found by setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2$$\end{document} in Proposition 4. This results in the system

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} m= 1:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad p_1^{(2)}+ \frac{c^2 + \imath \omega b }{\omega ^2} \Delta p^{(2)}_1 - c^2 \rho _0 n_0 v_1^{(2)} = h_1^{(2)}, \\ \text {(ii)} & \quad \frac{1}{\alpha _1} v_1^{(2)} + \mu p_1^{(2)} =0,\\ \end{aligned} \end{array}\right. } \\ m = 2:&{\left\{ \begin{array}{ll} \begin{aligned} \text {(i)} & \quad p_2^{(2)}+ \frac{c^2 + 2 \imath \omega b}{4 \omega ^2} \Delta p_2^{(2)} - c^2 \rho _0 n_0 v_2^{(2)} = h_2^{(2)} +\frac{ \beta _a }{ \rho _0 c^2} (p_1^{(2)})^2, \\ \text {(ii)} & \quad \frac{1}{\alpha _2} v_2^{(2)} + \mu p_2^{(2)} = \left( \zeta - 3 \xi \omega ^2 \right) (v_1^{(2)})^2, \end{aligned} \end{array}\right. } \end{aligned} \end{aligned}$$\end{document}

where, as before,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \alpha _m = \frac{1}{- m^2 \omega ^2 + \imath m \delta \omega _0 \omega + \omega _0^2} = \frac{ \omega _0^2- m^2 \omega ^2}{ (\omega _0^2- m^2 \omega ^2 )^2+ (m \delta \omega _0 \omega )^2} -\imath \frac{ m \delta \omega _0 \omega }{ (\omega _0^2- m^2 \omega ^2 )^2+ (m \delta \omega _0 \omega )^2} \end{aligned}$$\end{document}

for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \in \{1,2\}$$\end{document} . Above, we have rewritten \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta $$\end{document} as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\eta = \dfrac{ \beta _a }{ \rho _0 c^2}$$\end{document} to make the influence of the nonlinearity parameter in the medium \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _a$$\end{document} on the system explicit. In this two-harmonic setting, we can express \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_1^{(2)}$$\end{document} in terms of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_1^{(2)}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_2^{(2)}$$\end{document} in terms of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_2^{(2)}$$\end{document} :

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} v_1^{(2)}&= - \alpha _1 \mu \, p_1^{(2)}, \\ v_2^{(2)}&= \alpha _2 \left[ - \mu p_2^{(2)} + (\zeta -3\xi \omega ^2) {\alpha }_1^2 \mu ^2 \left( p_1^{(2)}\right) ^2 \right] . \end{aligned} \end{aligned}$$\end{document}

We can then use these expressions to eliminate \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_1^{(2)}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_2^{(2)}$$\end{document} for the first and third equations in (29). Multiplying the resulting equations by \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m^2 \omega ^2/c^2$$\end{document} yields the following system for the two pressure coefficients:

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} m = 1: \qquad&\frac{ \omega ^2}{\tilde{c}_1^2} p_1^{(2)} +\left( 1 + \imath \frac{ b \omega }{c^2} \right) \Delta p_1^{(2)} = \frac{\omega ^2}{c^2} h_1^{(2)},\\ m = 2: \qquad&\frac{4 \omega ^2}{\tilde{c}_2^2} p_2^{(2)} + \left( 1 +\imath \frac{2 b \omega }{c^2} \right) \Delta p_2^{(2)} = \frac{4 \omega ^2}{c^2} h_2^{(2)} + \frac{4 \omega ^2 }{\rho _0 c^4 } \left( \beta _a + \tilde{\beta }_a\right) \left( p_1^{(2)}\right) ^2, \end{aligned} \end{aligned}$$\end{document}

where

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{1}{\tilde{c}_{m}^2} = \frac{1}{c^2} + \rho _0 n_0 \mu \alpha _{m}, \qquad \tilde{\beta }_a= c^4 \rho _0^2 n_0(\zeta -3\xi \omega ^2)\mu ^2 \alpha _1^2 \alpha _2. \end{aligned}$$\end{document}

The above system clearly displays the influence of microbubbles on wave propagation, particularly through their impact on the effective wave speed and the nonlinearity parameter.

• From equations (30), we see that the effective wave number in the presence of microbubbles in a two-harmonic setting is

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \tilde{k}_m= \frac{m \omega }{\tilde{c}_m}, \quad \text {with } \tilde{c}_m \ \text { given in } (31). \end{aligned}$$\end{document}

Note that \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{k}_m$$\end{document} is in general complex since \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\alpha _m$$\end{document} is complex, and thus the attenuation in the wave propagation is due both to the acoustic dissipation via the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\imath \frac{m b \omega }{c^2} \Delta p_m^{(2)}$$\end{document} terms as well as the presence of microbubbles via the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\imath \Im (\tilde{k}^2_m) p_m^{(2)}$$\end{document} terms.

• In the second equation in (30), we see the influence of microbubbles on the nonlinearity parameter as their presence leads to the additional \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\beta }_a$$\end{document} term on the right-hand side, thus effectively enlarging the nonlinearity parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _a$$\end{document} . By setting \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{\beta } = 0$$\end{document} and replacing \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\tilde{c}^2_{m}$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$c^2$$\end{document} , the system reduces to the two-harmonic expansion of Westervelt’s equation for a strongly damped medium without bubbles.

Numerical Experiments

In this section, we investigate numerically the algorithms set up in Sections 3 and 5. To discretize the resulting Helmholtz problems, we use a conforming finite element method. We note that algorithmic frameworks of the multiharmonic-FEM type used in this section are known as the Harmonic Balance Finite Element Method (HBFEM) in the context of electromagnetism; see, e.g., [3, 4, 27].

The discretization of the classical Helmholtz equation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$-\Delta u-k^2 u=f$$\end{document} with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$k>0$$\end{document} , and the issues related to the so-called pollution effect and sign-indefiniteness have been extensively studied in the literature; see, e.g., [2, 19, 29, 30], and the references contained therein. Recall that the Helmholtz problems appearing in linearized problems in (25) and (28) have the general form

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\{ \begin{aligned}&-\left( 1+\imath \omega \frac{m b}{c^2} \right) \Delta u - (k^2+\mathfrak {a}) u + \imath \mathfrak {b} u= f_\omega \quad \text { in }\, \Omega ,\quad \text {with } k= \frac{m \omega }{c},\\&(\imath \omega m \beta + \gamma ) u + \nabla u \cdot \boldsymbol{n}= 0 \quad \text { on }\, \partial \Omega , \end{aligned} \right. \end{aligned}$$\end{document}

for a known right-hand side \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$f_\omega \in L^2(\Omega ; \mathbb {C})$$\end{document} (that depends on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega $$\end{document} ), and with

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \mathfrak {a}= \mu \rho _0 n_0 \frac{m^2 \omega ^2 (\omega _0^2 - m^2 \omega ^2)}{(\omega _0^2-m^2\omega ^2)^2 + (m \delta \omega _0 \omega )^2}, \quad \mathfrak {b}= \mu \rho _0 n_0 \frac{m^3 \omega ^3 \delta \omega _0 }{(\omega _0^2-m^2\omega ^2)^2 + (m \delta \omega _0 \omega )^2}. \end{aligned}$$\end{document}

The variational form of (32) is given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} a(u, \phi ) = l(\phi ) \quad \text {for all } \phi \in H^1(\Omega ; \mathbb {C}), \end{aligned}$$\end{document}

where the sesquilinear form \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a: H^1(\Omega )\times H^1(\Omega )\rightarrow \mathbb {C}$$\end{document} is defined as

$$[eqn]$$

and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$l(\phi ) :=\int _{\Omega }f \overline{\phi }\, d x$$\end{document} . We then have

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \mathfrak {R}\left( a(u,u)\right) = \Vert \nabla u\Vert ^2_{L^2(\Omega )} - (k^2+\mathfrak {a})\Vert u\Vert ^2_{L^2(\Omega )} + \left( \gamma - \frac{b\beta (m \omega )^2}{c^2}\right) \Vert u\Vert ^2_{L^2(\partial \Omega )}. \end{aligned} \end{aligned}$$\end{document}

and

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \begin{aligned} \Im \left( a(u,u)\right) = \frac{m \omega b}{c^2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} +\mathfrak {b}\Vert u\Vert ^2_{L^2(\Omega )}+\left( \gamma \frac{b}{c^2}+\beta \right) \omega m \Vert u\Vert ^2_{L^2(\partial \Omega )} \end{aligned} \end{aligned}$$\end{document}

Therefore, we can conclude that

$$[eqn]$$

which showcases how the coercivity of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a(\cdot , \cdot )$$\end{document} is influenced by the relative behavior of the various medium and frequency-dependent parameters. Note that from the complex part of (33), we obtain the bound

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{m \omega b}{c^2} \Vert \nabla u\Vert ^2_{L^2(\Omega )} +\frac{\mathfrak {b}}{2} \Vert u\Vert ^2_{L^2(\Omega )}+\left( \gamma \frac{b}{c^2}+\beta \right) \omega m \Vert u\Vert ^2_{L^2(\partial \Omega )} \le \frac{1}{2 \mathfrak {b}}\Vert f_\omega \Vert ^2_{L^2(\Omega )}, \end{aligned}$$\end{document}

which carries over to the finite element setting.

Numerical Framework

All simulations are performed using a conforming two-dimensional finite element discretization with linear Lagrange elements, using FEniCSx v0.9.0 (see, e.g., [5]) with Gmsh [13] employed for meshing. A direct LU factorization is used to solve the resulting linear systems without iterative refinement.1

We use a circular domain of propagation \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Omega = B_{0.2}(0) \subset \mathbb {R}^2$$\end{document} with a radius of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$0.2 \, {\text {m}}$$\end{document} and place a monopole source \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h=h(x)$$\end{document} at \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_0 = (0,0)^T$$\end{document} defined as

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} h(x) = {\left\{ \begin{array}{ll} \frac{a}{4 r_\delta } \left( 1 + \cos \left( \pi \frac{\left\| x - x_0 \right\| _{l^2} }{2 r_{\delta }}\right) \right) \qquad \quad & \left\| x - x_0 \right\| _{l^2} \le 2 r_{\delta }, \\ 0 & \text {otherwise.} \end{array}\right. } \end{aligned}$$\end{document}

We present numerical results for various values of a, using \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega =\omega _0$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\delta }=0.004$$\end{document} unless stated otherwise. The physical parameters used in the simulations are chosen as typical in ultrasound contrast imaging (see, e.g., [17, 18]), and are listed in Table 1.Table 1. Overview of the parameter values used in the simulationscspeed of sound1500 m/sbdiffusivity of sound \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 \cdot 10^{-3}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\rho _0$$\end{document} mixture mass density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1000\, kg/\hbox {m}^{3}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _a$$\end{document} nonlinearity coefficient3.5 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$R_0$$\end{document} initial radius \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \,\upmu \hbox {m}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0$$\end{document} bubble number density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1 \cdot 10^{12} \, 1/\hbox {mL} $$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$P_0$$\end{document} vapor pressure100 Pa \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\kappa $$\end{document} adiabatic exponent1.4 \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\nu $$\end{document} kinematic viscosity \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.9 \cdot 10^{-6} \hbox {m}^{2}/\hbox {s}$$\end{document}

Using the values in Table 1, we compute

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \omega _0 = \sqrt{\frac{3 \kappa P_0}{\rho _0 R_0^2}} \approx 0.324 \, {\text {M}\text {Hz}} \end{aligned}$$\end{document}

and additionally \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\delta = \frac{4 \nu }{\omega _0 R_0^2}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_0 = \frac{4 \pi }{3} R_0^3$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mu = \frac{4 \pi R_0}{\rho _0}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\zeta = \frac{(\kappa + 1) \omega _0^2}{2 v_0}$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\xi = \frac{1}{6 v_0}$$\end{document} are fixed. Furthermore, we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta = \frac{1}{c}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\gamma = 1$$\end{document} in the boundary conditions for Helmholtz problems.

The Influence of Spatial Discretization

For a fixed N, we first investigate the convergence of a quantity of interest given by

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \left\| \Re (p_{h_{FEM }}^{N}) \right\| _{L^{\infty }(0,T; L^2(\Omega ))}, \end{aligned}$$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{h_{FEM }}^N$$\end{document} denotes the approximate pressure, as the spatial discretization parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }\searrow 0$$\end{document} . A common rule of thumb for finite element discretizations of the classical Helmholtz equation with linear elements is to use at least 10 elements per wavelength in each spatial direction. If the mesh is not refined accordingly, one encounters the above-mentioned pollution effect, where numerical phase errors grow with increasing k; see, e.g., [2, 19, 29].Fig. 2. Relative error (34) of the quantity of interest with respect to the mesh size \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }$$\end{document} on a log–log scale

Figure 2 illustrates the behavior of the relative difference

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\begin{aligned} \frac{ \left| \bigl \Vert \Re (p_{h_{FEM }}^{N})\bigr \Vert _{L^\infty (0,T;L^2(\Omega ))} - \bigl \Vert \Re (p^{\textrm{ref}})\bigr \Vert _{L^\infty (0,T;L^2(\Omega ))} \right| }{ \bigl \Vert \Re (p^{\textrm{ref}})\bigr \Vert _{L^\infty (0,T;L^2(\Omega ))} }, \end{aligned}$$\end{document}

as the mesh is refined, using the solution on the finest spatial mesh as reference. A clear algebraic decay is observed for decreasing mesh size  \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }$$\end{document} , indicating that the spatial discretization error is well controlled in the considered parameter regime and that the chosen mesh resolutions are sufficient to avoid significant pollution effects for the frequencies under consideration. The similarity of the convergence curves for different truncation levels N shows that the spatial discretization error dominates over the effect of the multiharmonic truncation.Table 2. Mesh sizes, degrees of freedom, and computational time for different spatial discretizations \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }$$\end{document} NodesElementsTime0.0003906259,564,7951,933,5899984.8 s0.00078125239,024478,0471244.3 s0.001562560,136120,271145.5 s0.00312515,21830,43519.8 s0.006253,8967,7915.16 s0.01251,0102,0192.72 s0.0252795572.24 s0.05851692.15 s

Additional quantitative data are provided in Table 2, which lists the corresponding mesh sizes, number of nodes and elements, as well as the computational time required for each case. This highlights the significant increase in computational cost associated with finer meshes, especially for very small \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }$$\end{document} , where both the number of degrees of freedom and the runtime grow substantially. All simulations were carried out on a standard laptop (Intel Core i7, 16GB RAM, no dedicated GPU), underscoring the feasibility of the method even without access to high-performance computing resources.

In the following simulations, we use a mesh size of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }= 0.003$$\end{document} . Unless stated otherwise in the figure captions, we set \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_{\delta } = 0.004$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 10^5$$\end{document} , and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\omega = \omega _0$$\end{document} .

Comparison of Real- and Complex-Valued Fields

We compare the five-harmonic pressure obtained with the linearized multiharmonic cut-off algorithm using the complex-valued formulation (28) and the real-valued formulation (25). In the considered numerical setting, no noticeable differences between the resulting pressure fields are observed.Fig. 3. Five-harmonic expansion \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {R}(p_{h_{FEM }}^{(5)}(x_1,x_{2,0},t_0))$$\end{document} plotted as a function of \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_1$$\end{document} for a fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$x_2=x_{2,0}$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t=t_0$$\end{document} resulting from complex fields in (28) vs. real setting in (25)

To further quantify this observation, we compare the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\Omega )$$\end{document} -norms of the harmonic components obtained from the complex- and real-valued formulations. The relative differences in these norms are below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5 \cdot 10^{-3}$$\end{document} for the dominant harmonics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m = 1,\ldots ,4$$\end{document} , and below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2 \cdot 10^{-2}$$\end{document} for higher harmonics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m \ge 5$$\end{document} . At the same time, the amplitudes of these higher harmonics are already several orders of magnitude smaller than that of the fundamental mode.

We therefore consider both formulations to be numerically equivalent in the considered parameter regime. Since the complex formulation leads to a computationally more efficient algorithm, it is used in all subsequent simulations.

The Influence of the Number of Harmonics

We next analyze how many harmonics need to be retained in the multiharmonic expansion for the considered parameter settings. We first investigate the behavior of the approximate pressure field with respect to the truncation level N for a fixed \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$h_{FEM }$$\end{document} , which motivates the choice of a suitable reference solution.Fig. 4. Relative error (34) of the quantity of interest with respect to the truncation level N on a semi-log scale (logarithmic y-axis)

Figure 4 shows the relative error of the pressure field as defined in (34) for increasing truncation levels N, using the solution with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10$$\end{document} harmonics as reference. For larger values of N, the relative error levels off at approximately machine precision, reflecting that further changes in the norm are below the accuracy of the numerical discretization and floating-point arithmetic.

To illustrate how the contributions from individual harmonics are reflected in the approximate pressure field, we additionally consider the effect of truncating the multiharmonic expansion in the time domain. For a fixed reference time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} , we compute the pointwise differences

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ \left| p_{h_{FEM }}^N(x,t_0) -p_{h_{FEM }}^{(10)}(x,t_0) \right| , $$\end{document}

where \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$p_{h_{FEM }}^N(x,t)$$\end{document} denotes the approximate pressure field obtained using the first N harmonics and the solution with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10$$\end{document} harmonics serves as reference.Fig. 5. Pointwise differences \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$|p_{h_{FEM }}^N(x,t_0) - p_{h_{FEM }}^{(10)}(x,t_0) |$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^5$$\end{document} at a fixed time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} and different truncation levels \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=2,3,5,7$$\end{document}

Figure 5 visualizes the spatial structure of the truncation error for different truncation levels N in the strongly nonlinear case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^5$$\end{document} . For small values of N, the deviations from the reference solution are large and extend over significant parts of the domain, whereas increasing the number of retained harmonics leads to a rapid reduction and spatial localization of the differences. In particular, once approximately five harmonics are included, the remaining discrepancies become small throughout the domain, which is fully consistent with the harmonic-wise decay observed in Table 3.

We now fix the truncation index to \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10$$\end{document} and examine the harmonic content of the numerical solution in more detail. To this end, we compute the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\Omega )$$\end{document} -norms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert p^N_{h_{FEM }, m}\Vert _{L^2(\Omega )}$$\end{document} of the harmonic coefficients and normalize them by the fundamental mode,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ r_m(a) := \frac{\Vert p^N_{h_{FEM }, m}\Vert _{L^2(\Omega )}}{\Vert p^N_{h_{FEM },1}\Vert _{L^2(\Omega )}}, $$\end{document}

for each driving amplitude a and harmonic index m.

From Table 3 we observe that the relative amplitudes of higher harmonics decrease rapidly with increasing m, and that this decay becomes more pronounced as the driving amplitude a is reduced, indicating a weaker nonlinear response. For \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^3$$\end{document} , the dominant contribution comes from the first two harmonics, while for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^4$$\end{document} the first three harmonics contribute noticeably. In the strongly nonlinear case \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^5$$\end{document} , appreciable contributions extend to higher harmonic indices, indicating that retaining at most four to five harmonics is appropriate to accurately represent the pressure field.Table 3. Relative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\Omega )$$\end{document} -norms of higher harmonics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_m(a)$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$m=2,\dots ,5$$\end{document} and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10$$\end{document} a \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_2(a)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_3(a)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_4(a)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_5(a)$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^3$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.23\cdot 10^{-4}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.11\cdot 10^{-8}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.94\cdot 10^{-11}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.15\cdot 10^{-14}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^4$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.23\cdot 10^{-3}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.11\cdot 10^{-6}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.94\cdot 10^{-8}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.15\cdot 10^{-10}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^5$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.23\cdot 10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.11\cdot 10^{-4}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.94\cdot 10^{-5}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.15\cdot 10^{-6}$$\end{document}

Comparison of the Pressure Field with and Without Bubbles

Finally, we compare the pressure field obtained from the multiharmonic formulation in the presence and absence of bubble dynamics in order to assess how bubble coupling affects the nonlinear response of the system. We examine the pressure field in the time domain and the harmonic coefficients.Fig. 6. Pressure field \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Re (p_{h_{FEM }}^{(5)}(x,t_0))$$\end{document} at a fixed reference time \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$t_0$$\end{document} for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a=10^5$$\end{document} , with and without bubbles in the medium

Figure 6 compares the real part of the pressure distribution within the domain at a specific time point with and without the presence of bubbles in the medium; that is, the real part of the pressure obtained using the algorithm in (30) and the same algorithm with \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0=\tilde{\beta }=0$$\end{document} . The presence of bubbles leads to an overall damping effect, resulting in a reduced pressure amplitude across the domain. This observation aligns with the behavior predicted by the system of equations in (30), as the modified speed of sound introduces attenuation. Although the source amplitude is fixed in this case, we note that the differences become more pronounced for larger source amplitudes, leading to higher overall pressure levels.

To quantify these differences, we compare the harmonic content of the pressure field in both settings. For a fixed truncation level \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10$$\end{document} , we consider the \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\Omega )$$\end{document} -norms of the pressure harmonics \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\Vert p_{h_{FEM },m}^N\Vert _{L^2(\Omega )}$$\end{document} and report their relative contributions with respect to the fundamental mode,

\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ r_m := \frac{\Vert p_{h_{FEM },m}^N\Vert _{L^2(\Omega )}}{\Vert p_{h_{FEM },1}^N\Vert _{L^2(\Omega )}}. $$\end{document}

Table 4 shows that higher pressure harmonics have systematically larger relative amplitudes in the absence of bubbles. In particular, the decay of the relative harmonic ratios is significantly slower in the bubble-free case, with comparable contributions persisting over a broader range of harmonic indices. When bubble dynamics are included, the relative amplitudes of higher pressure harmonics decrease much more rapidly, dropping below \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$10^{-6}$$\end{document} already around the fifth harmonic. In addition, the harmonic coefficients associated with the bubble variable \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$v_m$$\end{document} are several orders of magnitude smaller than the corresponding pressure harmonics for all m, indicating that nonlinear effects transferred into the bubble dynamics are strongly damped.Table 4. Relative \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$L^2(\Omega )$$\end{document} -norms \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_m$$\end{document} of the pressure harmonics for \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$a = 10^5$$\end{document} with and without bubbles in the medium ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$N=10$$\end{document} ) \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_2$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_3$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_4$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_5$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$r_6$$\end{document} with bubbles \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.23\cdot 10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$8.11\cdot 10^{-4}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.94\cdot 10^{-5}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.15\cdot 10^{-6}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.13\cdot 10^{-7}$$\end{document} without bubbles \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5.48\cdot 10^{-2}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$1.78\cdot 10^{-3}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$5.50\cdot 10^{-5}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.61\cdot 10^{-6}$$\end{document} \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.64\cdot 10^{-6}$$\end{document}

We further examine the temporal evolution of the five-harmonic expansion at different spatial points in this context.Fig. 7. Real part of the pressure from the five-harmonic expansion computed using (28) without bubbles ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0=0$$\end{document} , dashed green line) and with bubbles (blue line) over time at three spatial pointsFig. 8Real part of the pressure \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathfrak {R}({p_{h_{FEM }}^{(5)}(0,0.05,t)})$$\end{document} from the five-harmonic expansion computed using (28) without bubbles ( \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0=0$$\end{document} , dashed green line) and with bubbles (blue line, scaled to the maximum of the green curve) over time

As illustrated in Figure 7, the presence of microbubbles near the source has little effect on the pressure waves. However, as the distance from the source increases, the signal strength decreases. The microbubbles not only attenuate the waveform but also enhance its nonlinear characteristics. To better visualize these effects, we scale the pressure waveform obtained from the five-harmonic expansion with bubbles to the maximum of the waveform without microbubbles, as shown in Figure 8.

Conclusion

The numerical results demonstrate that microbubbles introduce significant attenuation and phase shifts to the wave propagation, particularly at greater distances from the source. The impact increases with microbubble density and source amplitude, and is further influenced by the driving frequency. Regarding the number of harmonics N, in the considered numerical settings using three to five harmonics within the simplified multiharmonic framework obtained from complex fields provides a good compromise between computational efficiency and accuracy. This relatively low number of harmonics makes the approach promising for use in practical applications.

Outlook

In practice, nonlinear acoustic interactions between microbubbles and ultrasound waves generate not only harmonics, which are frequency components at integer multiples of the driving frequency, but also subharmonics, which appear at fractional multiples of the driving frequency, such as \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\omega }{2}$$\end{document} , \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\frac{\omega }{3}, \, \ldots $$\end{document} ; see, e.g., [26]. Subharmonics primarily appear due to non-spherical deformations and multibubble interactions. While these effects are not included in our current model, extending the framework to incorporate them would be an interesting direction for future research.

Expanding the present theoretical and numerical framework to explore other phenomena relevant to applications of focused ultrasound waves is also of practical interest. For instance, localized heating, cavitation, and nonlocal attenuation play an important role in therapeutic ultrasound applications such as targeted drug delivery, and it would be worthwhile to investigate the potential role of multiharmonic expansions in these modeling contexts.

Furthermore, investigating inverse problems related to reconstructing spatially varying parameters, such as the bubble number density \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$n_0=n_0(x)$$\end{document} or the nonlinearity parameter \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\beta _a=\beta _a(x)$$\end{document} from measured acoustic signals is an important task in the context of contrast-enhanced ultrasound imaging, as it could help improve diagnostic accuracy in the long run.