Homotopy analysis of the Lippmann-Schwinger equation for seismic wavefield modeling in strongly scattering media


 We present an application of the homotopy analysis method for solving the integral equations of the Lippmann–Schwinger type, which occurs frequently in acoustic and seismic scattering theory. In this method, a series solution is created which is guaranteed to converge independent of the scattering potential. This series solution differs from the conventional Born series because it contains two auxiliary parameters ϵ and h and an operator H that can be selected freely in order to control the convergence properties of the scattering series. The ϵ-parameter which controls the degree of dissipation in the reference medium (that makes the wavefield updates localized in space) is known from the so-called convergent Born series theory; but its use in conjunction with the homotopy analysis method represents a novel feature of this work. By using H = I (where I is the identity operator) and varying the convergence control parameters h and ϵ, we obtain a family of scattering series which reduces to the conventional Born series when h = −1 and ϵ = 0. By using H = γ where γ is a particular pre-conditioner and varying the convergence control parameters h and ϵ, we obtain another family of scattering series which reduces to the so-called convergent Born series when h = −1 and ϵ ≥ ϵc where ϵc is a critical dissipation parameter depending on the largest value of the scattering potential. This means that we have developed a kind of unified scattering series theory that includes the conventional and convergent Born series as special cases. By performing a series of 12 numerical experiments with a strongly scattering medium, we illustrate the effects of varying the (ϵ, h, H)-parameters on the convergence properties of the new homotopy scattering series. By using (ϵ, h, H) = (0.5, −0.8, I) we obtain a new scattering series that converges significantly faster than the convergent Born series. The use of a non-zero dissipation parameter ϵ seems to improve on the convergence properties of any scattering series, but one can now relax on the requirement ϵ ≥ ϵc from the convergent Born series theory, provided that a suitable value of the convergence control parameter h and operator H is used.


INTRODUCTION
There exist a range of different numerical methods for seismic wavefield modeling (Carcione, 2002), including differential equation methods (e.g., Robertsson et al., 2012) and integral equation methods (Pike and Sabatier, 2000;Oristaglio and Blok, 2012;Jakobsen, 2012;Jakobsen and Wu, 2016;Malovichko et al., 2018). The majority of researchers in the seismic community use differential equation methods (Carcione et al., 2002), but the integral equation approach has actually several advantages compared with the differential equation approach: (1) it is naturally target oriented (Haffinger et al., 2016;Huang et al., 2019), (2) it gives the sensitivity matrix directly in terms of Green's functions (Jakobsen and Ursin, 2015) which is convenient for uncertainty estimation (Eikrem et al., 2018) and (3) it is compatible with the use of domain decomposition and renormalization methods from modern physics Wu, 2016, 2018). However, the integral equation approach can be less efficient than the differential equation approach, depending on how it is implemented (Malovichka et al., 2018;Jakobsen and Wu, 2018;Jakobsen et al. 2019). An integral equation solution based on matrix inversion can be very accurate, but very memory-depending and costly (Jakobsen and Wu, 2018). Efficient implementations of the integral equations approach are typically based on the use of iterative methods and /or scattering series solutions Wu, 2016, Malovichko et al., 2018;Jakobsen et al., 2019, Huang et al., 2019b.
Many geophysicists are familiar with the scattering series of Born that one can easily obtain from the Lippmann-Schwinger equation via simple iteration (Jakobsen and Wu, 2016). However, the Born series represents an example of a so-called naive perturbation expansion which is only guaranteed to converge in the special case of relatively small contrasts (Kirkinis, 2008;Jakobsen and Wu, 2016). To obtain a convergent scattering scattering series in the presence of strong contrasts, it may be required to use a non-perturbative method for strongly nonlinear systems. Previously, researchers have developed convergent scattering series solutions of the Lippmann-Schwinger equation by using renormalization procedures (Abubakar and Habashy, 2003;Wu et al., 2007;Osnabrugge et al., 2016;Jakobsen et al., 2016;Jakobsen et al.. 2019). In this study, however, we have employed the so-called homotopy analysis method, which is based on concepts and ideas form topology (Liao, 2003;Hetmaniok et al., 2014).
The homotopy analysis method (HAM) used in this study was developed by Liao (1998Liao ( , 2003Liao ( , 2004Liao ( , 2009Liao ( , 2012Liao ( , 2014. However, the development of related globally convergent homotopy methods for solving nonlinear equations started around 1976 (see Watson, 1989). Historically, there have been several attempts to apply homotopy methods to model and invert geophysical data (see Watson, 1989), but the paper of Huang and Greenhalgh (2019a) appears to represent the first geophysical application of the modern homotopy analysis method developed by Liao (2003), which differs from the Convergent scattering series via homotopy analysis 3 one discussed by Watson (1989). In any case, the homotopy methods allows one to solve operator equations of any kind by using ideas and concepts of topology, which is a branch of pure and applied mathematics dealing with quantities that are preserved during continuous deformations. A homotopy describes a continuous transformation between two states and has been compared with the concept of scale-invariance in renormalization group theory (Palit and Datta, 2016;Jakobsen et al., 2019a;Pfeffer, 2019). The homotopy analysis method have been used to solve a range of different nonlinear problems, ranging from heat conduction problems (Abbasbandy, 2006) to problems within theoretical physics (Pfeffer, 2019). Most applications of the homotopy analysis method is based on a differential operator formulation, but there have also been successful attempts to solve integral equations of the Fredholm and Volterra types using the homotopy analysis method (see Hetmanio et al., 2014).
Although the homotopy analysis method may potentially be very useful for practical nonlinear inversion of seismic waveform data (see Han et al., 2005;Fu and Han, 2006), we shall focus on the forward problem. This is partially because there is still an important need for more work on the nonlinear direct scattering problem (Jakobsen et al., 2019a,b) and the corresponding nonlinear inverse scattering problem is much more difficult to solve due to its ill-posed nature. It will be demonstrated that the homotopy analysis method can be used to construct a scattering series solution of the Lippmann-Schwinger equation in the context of seismic wavefield modeling. Although such convergent scattering series have been developed on the basis of renormalization methods in the past (Abubakar and Habashy 1983;Osnabrugge et al., 2016;Jakobsen et al., 2019a), we think it is interesting to study convergence properties of the direct scattering series solution from different perspectives, since this may give us new ideas and insights that may be useful for future studies of nonlinear inverse scattering as well as direct scattering problems.
In section 2, we present fundamental equations and establish our notation. In section 3, we present a general method for obtaining convergent series solutions of nonlinear operator equations that does not depend on any parameter being small. In section 4, we derive a convergent scattering series solution

THE LIPPMANN-SCHWINGER EQUATION AND CONVENTIONAL BORN SERIES
The scalar wave equation in the frequency domain (the inhomogeneous Helmholtz equation) can be written as (Morse and Feshbach, 1953;Osnabrugge et al., 2016;Huang et al., 2019b,c) where k(x) is the wave number at position x. Following Osnabrugge et al. (2016), we now decompose the actual medium with wavenumber k(x) into an arbitrary homogeneous dissipative reference medium with complex wave number k d given by k 2 d = k 2 0 + i (where is an arbitrary small positive number) and a corresponding complex scattering potential V (x) (with compensating gain, rather than dissipation). It follows that where The second term on the right-hand side of equation (2) represents the so-called equivalent sources. By treating the contrast-sources just like real sources, one can derive the Lippmann-Schwinger equation (Jakobsen and Ursin, 2015) where G (0) (x − x ) is the Green's function for the homogeneous reference medium, that satisfies Note that the introduction of a non-zero imaginary part to the squared wave number k 2 0 in the homogeneous reference medium makes the energy associated with Green's function finite and the wave fields more localized (Osnabrugge et al., 2016;Jakobsen et al., 2019a;Huang et al., 2019c). Although most workers tend to set to zero, the use of a non-zero parameter improves the convergence properties of any scattering series (Abubakar and Habashy, 2003;Osbabrugge et al., 2016;Huang et al., 2019c;Jakobsen et al., 2019a).
In symbolic operator notation, the Lippmann-Schwinger equation (4) can be written as where the scattering potential operator V is local (but see Jakobsen and Wu, 2017) and can be represented by a diagonal matrix in the real-space representation (Jakobsen and Ursin, 2015). The above equation has the following exact formal solution: where I is the identity operator.
The solution (7) is valid independently of the contrast volume, but it involves the inversion of a huge operator or matrix (in the coordinate representation), which can be very costly in the case of a realistic model. In principle, one could try to solve the Lippmann-Schwinger equation by iteration. The well-known Born series can be regarded as the simplest possible iterative solution of the Lippmann-Schwinger equation and can be presented as where ψ 0 = ψ (0) and The Born series is very popular due to its simplicit. However, the Born series represents an example of a naive perturbation expansion (Kirkinis, 2008) which is only guaranteed to converge if the contrast is relatively small, in the sense that the largest eigenvalue of the operator G (0) V must be smaller than unity (Weinberg, 1983, Newton, 2002Osnabrugge et al., 2016).

THE HOMOTOPY ANALYSIS METHOD
The homotopy analysis method can be used to solve operator equations of the form (Liao, 2003) where N denotes a nonlinear operator and ψ is the unknown function (or state vector). The first step is to define the homotopy operator H by (Liao, 2003) where λ ∈ [0, 1] is the so-called embedding parameter, h = 0 is the so-called convergence control parameter, H is a convergence control operator (see section 4), ψ 0 is our initial guess of the solution to equation (10) and L is an auxiliary linear operator that can be selected arbitrarily as long as L [0] = 0.
By setting H [Φ, λ] = 0 we get the so-called zero-order deformation equation (Liao, 2003) ( where ψ is the solution of equation (10) we are looking for. A gradual change in the embedding parameter λ from 0 to 1 therefore means a continuous transition of Φ(λ) from the initial guess ψ 0 to the exact solution ψ of the original equation (10).
If we now expand the auxiliary field Φ(λ) into a Maclaurin series with respect to the embedding parameter λ then we obtain (Liao, 2003) By introducing the definition (Liao, 2003) the above equation (13) can be expressed as (Liao, 2003) If the above series (15) is convergent for λ = 1 then the solution we are looking for is given by (Liao, It is of course not obvious that the series (15) is convergent for λ = 1, but by adjusting the auxiliary parameter h and the auxiliary operator H we can make sure that this series is indeed convergent (Liao, 2003).
In order to determine the different ψ m terms , we now differentiate the left and right side of the 0th-order deformation equation (12) m times with respect to the auxiliary parameter λ, divide the result by m! and set λ = 0. In this way suggested by Liao (2003), we obtain the so-called mth-order deformation equation (m > 0): where and The different R m parameters will depend on the nature of the non-linear operator N . In the next section, we shall evaluate the R m -parameters for the nonlinear operator corresponding with the Lippmann- By using the above definitions of the linear and nonlinear operators L and N in conjunction with the mth-order deformation equation 17, we obtain By using the definition of the nonlinear operator N given in equation (20) in conjunction with the expression for the R m parameters in equation (18), we get The above equation implies that or By using the above expression for R m in conjunction with the recursive formula (21), we obtain and for m ≥ 2: where Equations (25) and (A.1) for the first and higher-order terms in the homotopy analysis scattering series represents the main results of this paper. This homotopy analysis method iterative solution of the Lippmann-Schwinger equation differs from the conventional Born series via the convergence control parameter h and the operator H that can be selected arbitrarily to ensure that the series is convergent.
The HAM series converges if the spectral radius σ of M is smaller than unity; which can occur even if the spectral radius of the operator G 0 V is larger than unity; that is, when the conventional scattering series of Born diverges.

COMPARISON WITH EXISTING ANALYTICAL RESULTS
If we use our freedom to set ψ 0 = ψ(0), h = −1 and H = I then the homotopy series in equations (25-26) reduces to the conventional Born series (9). As discussed earlier, the conventional Born series have a rather small range of convergence, since the largest eigenvalue of the operator G (0) V must be smaller than unity.  (23) that If we now set h = −1 and H = γ then the above equation becomes where Also, it follows from equation that control parameters and h we can accelerate the convergence of the HAM series. This point will be illustrated in the next section dealing with numerical experiments based on a strongly scattering seismic model.By introducing an imaginary part to the wavevector of the background medium, we make the total energy represented by the background Greens function finite and localized. The imaginary term in the background medium is compensated exactly by an imaginary term in the scattering potential.
Therefore, the final solution remains the same as the solution without any dissipation.

NUMERICAL RESULTS AND DISCUSSION
We In each experiment, we used one of the combinations of the ( , h, H) -parameters given in Table 1 and generated a scattering series solution of the Lippmann-Schwinger equation by using the recursive formula (25-27).
We quantified the convergence properties of the different scattering series by calculating the normalized overall error δ k as a function of the number of iterations k, where and ψ (r) is a reference solution obtained by solving equation (7) via matrix inversion (see Figure   2), which is exact apart from very small numerical discretization errors (Jakobsen, 2012). In each numerical experiment, we iterated until the normalized overall error became smaller than 10 −3 (in the case of convergence) or larger than 10 (in the case of divergence). However, this stopping criteria is of course flexible and dependent on the desired accuracy.
If the scattering series diverges then the resulting wavefield ( Figure 3) will of course look very different from the reference wavefield ( Figure 2). If the scattering series converges in the sense that the overall normalized error becomes smaller than 10 −3 than the resulting wavefield ( Figure 4) will necessarily be very similar to the true wavefield ( Figure 2). Since the resulting wavefield is independent of the auxiliary parameters in the case of convergence, we focus on the behaviour of δ k rather than the wavefield itself.
In numerical experiments 1-6 ( Figure 5) we assumed H = I and varied the dissipation parameter and the convergence control parameter h. As discussed in the previous section, when = 0 and h = −1 the numerical results correspond with the conventional Born series, whereas the use of differentand h-values represents different modifications of the conventional Born series. Clearly, one can see from the blue curve in Figure 5 that the conventional Born series corresponding with = 0 and h = −1 diverges for this strongly scattering medium. When = 0 and h = −0.95 corresponding with the green curve in Figure 5, the scattering series still diverges. When = 0 and h = −0.9 corresponding with the red curve in Figure 5, the scattering series is starting to converge, but extremely slowly. When = 0 and h = −0.8 corresponding with the cyan curve in Figure 5, the scattering series converges faster. When = 0 and h = −0.1 corresponding to the black curve in Figure 5, the scattering series is still convergent, but the convergence rate is much smaller than when using h = −0.8. When = 0.5 and h = −0.8 corresponding to the black curve in Figure 5 than the scattering series converges faster than for all the other experiments 1-5. Therefore, it appears that the use of a non-zero -value in conjunction with an optimal h-value helps to accelerate an already convergent scattering series.
In numerical experiments 7-12 (  Figure 6, the scattering series is as expected divergent. When = 0.5 c and h = −0.5 corresponding with the red curve in Figure 6, the scattering series is still divergent. However, when = 0.5 c and h = −0.25 corresponding with the cyan curve in Figure 6, the scattering series become convergent again. When = 0.5 c and h = −0.125 corresponding with the curve in Figure 6, the scattering series is still convergent, but the convergence rate is smaller than when using h = −0.25. When using = 0.25 c corresponding with the black curve in Figure 6   The auxiliary parameters and h may be referred to as local and global convergence control pa-rameters, respectively. This is because a non-zero value leads to dissipation in the reference medium (and gain in the scattering potential), which makes the wavefield update more localized in space; and different h-values are associated with different degrees of global wavefield scaling. Having both local and global convergence control parameters in addition to the auxiliary convergence control operator H makes this homotopy analysis method very general and flexible.

CONCLUDING REMARKS
We have used the homotopy analysis method (HAM) to derive a general scattering series solution of the Lippmann-Schwinger equation which is guaranteed to converge independent of the scattering potential, provided that one select the dissipation parameter as well as the convergence control pa- this approach by introducing subseries. We hope to develop convergent inverse scattering series using the HAM scattering series. However, it is not obvious that the inverse scattering series will converge even though the forward scattering series is convergent. Finally, we note that the work reported here may be useful in future applications of the homotopy analysis method within the context of nonlinear inverse scattering to solve the so-called regularized normal equations (Jegen et al., 2001).

ACKNOWLEDGMENTS
We thank the Research Council of Norway for funding related to Petromaks II project 267769/E3 (Bayesian inversion of 4D seismic waveform data for quantitative integration with production data) and the National IOR Centre of Norway. Ru-Shan Wu would like to thank the sponsors of the WTOPI consortium for their continuous support.   Table 1.  Table 1.  Table 1.

APPENDIX A: DESCRIPTION AND IMPLEMENTATION OF THE HAM SERIES
Similar to the conventional Born series, every iteration is associated with multiple scattering processes of different orders. However, we have reorganized the different terms in the conventional Born series so that the spectral radius of M is smaller than unity. This implies that each new term is smaller than the previous one, so that the scattering series does not diverge when the number of iterations becomes large. Mathematically, this is done by introducing an integral operator with spectral radius smaller than unity via the use of control parameter h and the convergence control operator H. The series converges if the spectral radius of the operator M is less than unity. We have also introduced an element of dissipation in the reference medium, which ensures that the energy associated with Greens function is finite and localized. It should be emphasized that the dissipation parameter can be selected arbitrary.
This is because the dissipation is compensated exactly by a corresponding gain term in the scattering potential, suggesting that the final results are independent of this dissipation in the reference medium.
The dissipation aspect of our HAM algorithm is similar to the convergent Born series of Osnabrugge et al. (2016). However, our convergent scattering series is more general than the convergent Born series, since the convergent control parameters can be selected rather arbitrary, as long as the spectral radius of the M-operator is smaller than unity. Some details for implementation of the new convergent scattering series using based on the HAM is provided in Algorithm 1. In addition, Table 1 shows the norm M and spectral radius σ (M ) of the operator M with numerical experiments.
The HAM algorithm is represented by equations (25)-(27). However, the formulation in the main text is based on the real-space coordinate representation of the relevant integral operators. As discussed by Osnabrugge et al., (2016), the operation of Green's function with contrast-source terms has a convolution structure that can be implemented more efficiently by using the wave vector representation; that is, by using the Fast Fourier Transform (FFT) algorithm in this context. This is because convolution in real-space is equivalent to multiplication in the Fourier space, and the computational cost of the FFT-operation is much smaller than that of matrix multiplication and inversion.The memory requirements scales like N 2 and N when using the position and wave vector representations, respectively.
The computational cost should theoretically scale like N 3 and N logN when using the position and wave vector representations, respectively. The iterative FFT algorithm is implemented as where fft and ifft are the forward and inverse fast Fourier transform operators, k is the Fourier transformed coordinates. Figures A1 and A2 show the frequency domain wavefield using exact solutions, the efficient and memory-saving FFT implementation with 100 iterations, respectively.
Algorithm 1 Pseudo code for the scattering series Initialisation: frequency, maximum iteration number N max , the parameter true model and background model;  Table 1.