Higher-order finite element methods for the nonlinear Helmholtz equation

In this work, we analyze the finite element method with arbitrary but fixed polynomial degree for the nonlinear Helmholtz equation with impedance boundary conditions. We show well-posedness and error estimates of the finite element solution under a resolution condition between the wave number $k$, the mesh size $h$ and the polynomial degree $p$ of the form ``$k(kh)^{p+1}$ sufficiently small'' and a so-called smallness of the data assumption. For the latter, we prove that the logarithmic dependence in $h$ from the case $p=1$ in [H.~Wu, J.~Zou, \emph{SIAM J.~Numer.~Anal.} 56(3): 1338-1359, 2018] can be removed for $p\geq 2$. We show convergence of two different fixed-point iteration schemes. Numerical experiments illustrate our theoretical results and compare the robustness of the iteration schemes with respect to the size of the nonlinearity and the right-hand side data.


Introduction
In various situations, such as for high intensities, linear(ized) material laws are no longer accurate enough and nonlinear constitutive relations have to be incorporated into the models.One well-known example are Kerr-type materials [12] in electromagnetics, where the permittivity ε depends on the electric field E like ε(E) = ε 0 + ε 2 |E| 2 .In general, wave propagation in nonlinear media causes a lot of new possible phenomena such as optical bistability [9].
As simplified model, we study in the following the nonlinear Helmholtz problem where D ⊂⊂ Ω is the subdomain where the nonlinearity is "active".Detailed assumptions on the domain and the data are given below.Problems of the above form occur in nonlinear acoustics as well as in time-harmonic and suitably polarized nonlinear electromagnetics.The nonlinear Helmholtz equation has been studied analytically for instance in [6].Various numerical approaches have been suggested as well: [1] and [25] consider layered media and study a finite volume approach or approximate it as the steady state to a Schrödinger equation, respectively.[26] focuses on different iteration schemes for the nonlinearity and uses a pseudospectral method in space.A multiscale finite element method is proposed and analyzed for the heterogeneous nonlinear Helmholtz equation in [16].The present work is inspired by [24], where the linear (i.e., p = 1) finite element method (FEM) is studied and a priori error estimates are shown.The main findings of [24] are that, under a smallness of the data assumption between k, ε, f and g as well as the resolution condition k 3 h 2 sufficiently small, a unique finite element solution exists and the discretization error is of the order kh + k 3 h 2 .We emphasize that both the resolution condition as well as the a priori error estimate are similar to the well-studied linear case in the so-called pre-asymptotic regime, cf.[5].In a similar spirit as [24], the recent work [11] provides a finite element error analysis of the nonlinear Helmholtz equation with perfectly matched layer at the boundary and Newton's method as iteration scheme.
In fact, a major ingredient in the numerical analysis of [24,11] is the study of an auxiliary linearized Helmholtz problem which is solved in each iteration step.The finite element error analysis of the linear Helmholtz equation is much more mature than of its nonlinear counterpart.Seminal results are the asymptotic hp-FEM analysis of [18,19] and the pre-asymptotic error analysis for arbitrary, but fixed polynomial degree of [5].These results have been obtained for the constant coefficient case, but recently much progress has been made for the heterogeneous Helmholtz equation as well.In the asymptotic regime, arbitrary but fixed polynomial degree is treated in [3] and the hp-FEM in [10,14,2].Pre-asymptotic estimates for the absolute error can be found in [7,20], whereas [13] studies the relative error.By the difference between "arbitrary but fixed polynomial degree" and "hp-FEM" results, we mean that in the first case, constants may (implicitly) depend on the polynomial degree.By now higher-order-and hp-FEM approximations are the state of the art -in comparison to linear FEM -for the Helmholtz equation as they allow a relaxed resolution condition of k(kh) 2p sufficiently small (pre-asymptotic, fixed polynomial degree) or kh/p sufficiently small and p ln k (asymptotic, hp version).
Our main contribution is the rigorous a priori error analysis of higher-order finite element methods for the nonlinear Helmholtz problem.We essentially show that under a smallness of the data assumption similar to [24,11], a resolution condition of k(kh) p+1 is sufficient for existence and uniqueness of a finite element solution and the discretization error is of the order h + (kh) p + k(kh) 2p .We also rely on the numerical analysis of a linearized Helmholtz equation, where we prove a solution splitting (into an analytic and a less oscillatory part) in the spirit of [18,19] as well as pre-asymptotic stability and error estimates.The linearized Helmholtz equation has a non-constant, discontinuous refractive index, so that we cannot directly apply recent results for the heterogeneous Helmholtz equation [15,2,20].If we assume smoothness of the domain where the nonlinearity is active, one might transfer the recent results from [7].However, here we use a perturbation argument that the deviation from the constant-coefficient case is sufficiently small due to the smallness of the data assumption.Moreover, we show a discrete stability result in the L ∞ -norm which causes the tighter resolution condition k(kh) p+1  1 compared to the linear case.Along this analysis, we treat two different iteration schemes: the frozen nonlinearity scheme considered in [24] and a scheme suggested in [26].For the latter, we provide the first proof of linear convergence by interpreting it as a fixed-point iteration.This result fills a theoretical gap in [26] and may be of own interest.
The paper is organized as follows.We present the setting and study the iteration schemes in Section 2. Section 3 describes the finite element discretization and the main error estimates, whose proofs are then presented in Section 4. Finally, we illustrate our theoretical results with numerical experiments in Section 5, where we also compare the two iteration schemes numerically in detail.

Nonlinear Helmholtz equation in the continuous setting
In this section, we formulate our model problem and discuss the solution of the nonlinear problem via iteration schemes in the continuous setting.
Throughout this article, all our functions are complex-valued unless otherwise mentioned.For any (sub)domain S, (•, •) S denotes the complex L 2 -scalar product (with complex conjugate in the second argument).We use standard notation on Sobolev spaces H s (S) and their norms.Further, we use the following (semi)norms • 0,S := • L 2 (S) , | • | 1,S := ∇ • 0,S , and • 2,S := • H 2 (S) .As usual in the Helmholtz context, we also employ the following k-weighted norm We will omit the subdomain S in the notation of norms and scalar products if it equals the full computational domain Ω and no confusion can arise.Last, we use the notation a b to indicate that there exists a generic constant C, independent of h and k but possibly dependent on the polynomial degree p, such that a ≤ Cb.

Model problem
Let Ω ⊂ R d , d ∈ {2, 3}, be a bounded star-shaped domain with analytic boundary Γ = ∂Ω and outer normal ν.In this work, we are interested in approximating the (weak) solution u ∈ H 1 (Ω) of the following nonlinear Helmholtz problem for all v ∈ H 1 (D).Here, k is the wave number, ε the Kerr coefficient, and χ D denotes the characteristic function of D. We make the following assumptions on the data throughout the whole article.Assumption 2.1.We assume that f ∈ L 2 (Ω), g ∈ H 1/2 (Γ) and ε ∈ R >0 .Suppose that k 1 in the sense there exists a constant k 0 > 0 such that k ≥ k 0 and, subsequently, all constants in our estimates may depend on k 0 .Finally, we assume that D ⊂⊂ Ω is a non-empty compactly embedded subdomain with Lipschitz boundary.
In the following, we abbreviate C data := f 0 + g H 1/2 (Γ) .[24] shows that there exists θ 0 such that if there exists a unique solution u ∈ H 1 (Ω) to (2.1).Further, u satisfies the following a priori estimates The (sufficient) condition k d−2 εC 2 data < θ 0 for these results to hold is called a smallness of the data assumption and it comes from a Banach fixed-point argument, see also Section 2.2 below.While the exact condition itself does not seem to be sharp in numerical experiments, it is well known that some condition on k, ε, f and g is required for uniqueness.Note that the power of k in the smallness assumption depends on the stability constant for the linear Helmholtz equation with ε = 0. Throughout this article, we assume this stability constant to be O(1), which is well established for star-shaped domains Ω, see, e.g., [17,4].Remark 2.2.For some results, the assumption g ∈ L 2 (Γ) instead of g ∈ H 1/2 (Γ) would be sufficient.For simplicity, we omit to track this and work under Assumption 2.1 and with the constant C data .

Iteration schemes
We present and discuss two iteration schemes for the nonlinear problem in the continuous setting.Both schemes will subsequently be combined with the spatial discretization in Section 3.1 to obtain a discrete solution in practice.
[24] considers the following fixed-point iteration based on a frozen nonlinearity approach.Given the previous iterate u (l−1) , the next iterate u (l) ∈ H 1 (Ω) is defined as where The iteration starts from some u (0) ∈ H 1 (Ω) ∩ L ∞ (D) with sufficiently small energy norm.For simplicity, we consider u (0) ≡ 0 throughout.Under the smallness of data assumption (2.2) with suitably chosen θ 0 , the sequence {u (l) } l∈N forms a strict contraction in the sense that This is the main idea in the proof of existence and uniqueness of a solution u to (2.1) in [24].Precisely, the sequence {u (l) } converges to u strongly in H 1 (Ω).[26] proposes a different iteration scheme -the motivation stems from Newton's method, but it is neither Newton's method itself nor any variant thereof.Given the previous iterate u (l−1) , the next iterate u (l) ∈ H 1 (Ω) is defined as (2.6) Again, we let the iteration start from u (0) ≡ 0 for simplicity.[26] observes numerically that the scheme converges with a linear rate and that it has the advantage of allowing larger values of k, ε and the data than (2.4).To the best of our knowledge, these observations have not been rigorously confirmed in theory.In the rest of this section, we will hence address the convergence of (2.6) and relate it to the convergence of the frozen nonlinearity approach.

Linearized Helmholtz equation as auxiliary problem
Both iteration schemes (2.4) and (2.6) solve a linear(ized) Helmholtz problem in each step, which we can formulate as follows.Let Φ ∈ H 1 (Ω) ∩ L ∞ (D) be given.For (2.4), find w ∈ H 1 (Ω) such that for all v ∈ H 1 (Ω).For (2.4), find w ∈ H 1 (Ω) such that for all v ∈ H 1 (Ω), where f = f − k 2 εχ D |Φ| 2 Φ.Note that (2.7) and (2.8) are closely related because A(Φ; w, v) = B lin ( √ 2Φ; w, v) and the right-hand side is (slightly) different.Consequently, we can deduce many results for (2.8) from their counterparts for (2.7).We emphasize that (2.7) is of Helmholtz-type, but with a variable, i.e., x-dependent, refractive index n := 1 + χ D ε|Φ| 2 induced by Φ.Moreover, our assumptions on ε and D imply that n is discontinuous over ∂D, i.e., the interface between the nonlinear material and the linear "background".[24] shows that there exists a constant θ 1 such that if kε Φ 2 L ∞ (D) ≤ θ 1 , a unique solution w ∈ H 1 (Ω) to (2.7) exists and it satisfies the a priori estimates (2.9) In fact, these results on w are the crucial ingredient to establish existence, uniqueness and a priori estimates (cf.(2.3)) for the solution u to the nonlinear problem (2.1) in [24].By exploiting the correspondence between (2.7) and (2.8), we directly obtain that, if kε Φ 2 L ∞ (D) ≤ θ 1 /2, a unique solution w ∈ H 1 (Ω) to (2.8) exists.From (2.9) and (2.8), we deduce the a priori estimates as well as These estimates will play a central role in the convergence proof for (2.6) in Section 2.4 below.
In the analysis of the finite element method for the linearized Helmholtz problem in Section 3.2, we need a splitting of the continuous solution w into an H 2 -regular and an analytic part.The recent results on solution splittings obtained by [14,15] are not applicable since n may exhibit a discontinuity across ∂D in our case, see above.Under a smallness of the data assumption, we show that the well-known solution splitting of the standard Helmholtz equation with ε = 0 already yields the desired result for (2.7) as well.We note once more that the argument then transfers to (2.8) by the obvious modifications in the scaling of Φ and the form of f .Proposition 2.3.There is a constant θ 2 such that, if kε Φ 2 L ∞ (D) ≤ θ 2 , the solution w of (2.7) can be split as w = w H 2 + w A with w H 2 ∈ H 2 (Ω) and w A analytic.Further, there exist k-and Φ-independent constants C, γ > 0 such that for any m ∈ N 0 The proof is presented in the appendix.

Convergence of scheme (2.6)
We now show that scheme (2.6) satisfies a contraction property and therefore, the iteration sequence converges to the (unique) solution u of (2.1).
for some q < min{ 1 6 , C1θ1 2 } with θ 1 introduced in Section 2.3, we have for iteration scheme and the sequence {u (l) } l converges linearly to the solution u of (2.1).
Proof.First step: A priori estimates for u (l) .We show by induction that for all l ≥ 1, u (l) is well-defined and satisfies The case l = 1 directly follows from u (0) = 0 and (2.10)- (2.11).Let the statement be satisfied for l.Since kε u 2 by the assumptions, the discussion in Section 2.3 yields that u (l+1) is indeed well-defined.Moreover, we deduce from (2.10) that which recursively yields with q < 1 that Employing (2.11), we furthermore obtain which finishes the first step.
Second step: Contraction property.Direct calculation shows that u (l+1) − u (l) solves The estimates from the first step yield where we used (2.12) in the last step.The assumption q < 1/6 finishes the proof.
The proposition explains the linear convergence observed in practice [26].However, the required smallness of the data assumption is more restrictive than for the frozen nonlinearity scheme.Precisely, following the proofs of [24], we see that the contraction property (2.5 , which is more relaxed in comparison to (2.12).Hence, Proposition 2.4 does not explain the better "robustness" of the scheme with respect to the data observed in [26] as well as in our experiments in Section 5.

Nonlinear Helmholtz equation in the discrete setting
In this section, we turn to the finite element approximation of (2.1).We introduce the discretization using finite elements with higher-order polynomials in Section 3.1.We then present the results of a priori error analysis, where we first consider the linearized problems in Section 3.2 and then the nonlinear problem in Section 3.3.All proofs are collected in Section 4 and the appendix.

Finite element discretization and notation
Since we assume Γ to be analytic, we will consider curved elements in order to have a conforming discretization.We follow the typical procedure as outlined in, e.g., [18,Section 5].We assume that there exists a polyhedral/polygonal domain Ω and a bi-Lipschitz mapping ξ : Ω → Ω.Let T h denote an admissible, shape regular simplicial mesh of Ω.We assume that the restrictions ξ| T are analytic for all T ∈ T h .We then set T h = {ξ( T ) : T ∈ T h } as our mesh on Ω with mesh size h := max T ∈T h diam T .Note that for any T = ξ( T ) ∈ T h , there exists an affine, bijective mapping A T : T → T from the reference element T (the unit simplex).Consequently, we have a mapping F T : T → T via F T = R T • A T with R T = ξ| T .We assume F T , R T and A T to satisfy the smoothness and scaling assumptions of [18,Assumption 5.2].
For such a so-called quasi-uniform regular simplicial mesh T h , we denote the finite element space of piecewise (mapped) polynomials of degree p by V h,p , i.e., where P p denotes the polynomials of degree p.We now seek the discrete solution for all v h ∈ V h,p .This yields a nonlinear system that we can solve via the discrete versions of the iteration schemes (2.4) or (2.6).As usual, these discrete versions are obtained by a Galerkin procedure, i.e., ansatz as well as test functions come from the space V h,p .As already done in the continuous case, we start the iterations with u (0) h,p ≡ 0 for simplicity.We collect further finite element-related notation that will turn out useful in the error analysis.Let P h : H 1 (Ω) → V h,p be the elliptic projection as defined by [24] via This projection is well-defined and satisfies Further, following [5], we introduce discrete H j (Ω)-norms on V h,p .We define the discrete operator denote its eigenvalues, which are all positive, and let ϕ j,h for j = 1, . . .dim V h,p be the corresponding discrete eigenfunctions.For any real number j, the operator A j h is defined via The discrete norms on V h,p are then defined for any integer j via 2. for any integer 0 ≤ j ≤ p + 1,

FEM error analysis for the auxiliary problem
We can now analyze the linear auxiliary problem in the discrete setting.According to the discussion in Section 2.3, we focus on problem (3.8) below using B lin , but note that everything carries over to the discrete version of (2.8) by the relation of A and B lin .Let Φ ∈ L ∞ (D)∩H 1 (Ω) be given.Define w h ∈ V h,p as the solution of with the constants introduced in Section 2, there exists a constant C 0 > 0 such that, if k(kh) 2p ≤ C 0 , the finite element solution w h to the auxiliary problem (3.8) exists, is unique and satisfies where the constant in may depend on p. Further, w h is stable in the following sense In the proof of the above lemma, we will also establish the following L 2 -error estimate Lemma 3.1 essentially transfers [5] to a case where the coefficient n in the Helmholtz problem is no longer constant.In contrast to the approach in [20, Sec.2.4] and [7], we treat n as a sufficiently small perturbation from the constant coefficient case.Therefore, some of our assumptions are different, in particular we can allow for lower regularity in n.Further, since we do not have a coefficient in the gradient part, i.e., A = 1, we can use Robin boundary conditions everywhere, cf. the discussion in [20, Rem.2.62].Concerning the occurrence of kh as first term in (3.9), we refer to the discussion after Theorem 3.5.For convenience, we include the proof of Lemma 3.1 (along the lines of [5]) in the appendix.
Besides the stability of w h in the energy norm, we have the following L ∞ -estimate, which is important for the nonlinear case.
with the constants introduced in Section 2, there exists a constant C 1 > 0 such that, if k(kh) p+1 ≤ C 1 , the solution w h of (3.8) satisfies where p = 1 for p = 1 and p = 0 for p ≥ 2.
The proof follows the lines of [24,11] using the interior L ∞ estimates of [21].The latter only introduce an | ln h|-dependence in the case p = 1.Note that [11] recently showed that even in the linear case, the | ln h| can be removed if d = 2. Remark 3.3 (Why the resolution conditions in Lemmas 3.1 and 3.2 differ).Comparing to the linear Helmholtz case, the resolution condition in Lemma 3.1 is the known one from [5], while the condition in Lemma 3.2 corresponds to an older, sub-optimal condition in [27].To get the improved result similar to Lemma 3.1, [5] uses suitable discrete H j -norms, cf.(3.6), and negative Sobolev norms in the estimation of L 2 -scalar products (between a discrete and a projection error).In the L ∞ -norm estimate, however, we do not get such L 2 -scalar products and therefore cannot use this technique.This is the main reason for the tighter resolution condition in Lemma 3.2.Note that both conditions agree for p = 1, in particular the condition in Lemma 3.2 agrees with the result for p = 1 in [24].

Finite element method for the nonlinear problem
We are now prepared to analyze the higher-order finite element method for the iteration schemes (2.4) and (2.6) applied to the nonlinear Helmholtz equation.We emphasize once more that the analysis in [24,11], which partly inspires our proofs, is limited to p = 1 and either iteration scheme (2.4) or Newton's method.We start with the convergence of the discrete schemes and, thereby, existence and uniqueness of the solution u h,p to (3.1).
, where C1 , C2 associated with schemes (2.4) and (2.6), respectively, are some constants.If σ j < 1, the associated sequence {u h,p ≡ 0 converges to the unique solution u h,p of (3.1) with rate σ j , i.e., Further, u h,p satisfies the stability estimates We hence obtain the continuous as well as the discrete solution as limit of a sequence of solutions to linearized Helmholtz problems.The proposition gives us the convergence of the two iteration schemes also in the discrete setting as well as existence and uniqueness of the solution to (3.1).Of course, this unique solution exists as soon as σ j < 1 for j = 1 or j = 2.In other words, we can choose the less restrictive condition when we consider properties of the discrete solution to the nonlinear problem.Using the error estimates in the linear case (cf.Lemma 3.1), we can conclude our main result on the finite element error.Note that by combining the previous theorem and (3.13) we deduce an error estimate for u − u (l) h,p and any of the two schemes (2.4) or (2.6) by the triangle inequality.Theorem 3.5 bounds the error between the exact and the discrete solution of the nonlinear Helmholtz equation, where the dependence on k, h, and p is the same as in the linear case.We provide a pre-asymptotic error bound with the so-called pollution term k(kh) 2p under a resolution condition discussed further in Remark 3.3.Note that the first term h occurs in Theorem 3.5 since we do not assume more than H 2 (Ω)-regularity of u.Nevertheless, this term can be interpreted as a higher order term in the sense that it does not dictate the rate of convergence due to the resolution condition, cf.[20, Rem.2.40].As in [5] for the linear Helmholtz equation, our result assumes a fixed polynomial degree in the sense that the involved constants depend on p.We strongly believe that the famous hp-error analysis in the asymptotic regime [18,19] can be transferred to the nonlinear Helmholtz equation for sufficiently small data as well.An hp-version of the result in Lemma 3.1 in the asymptotic regime is already available, see [14], but to the best of our knowledge, nothing is known about an hp-version of the L ∞ -estimate in Lemma 3.2.One can circumvent the application of Lemma 3.2 as in [16], but the price to pay is a worse k-dependence in the smallness of the data assumption.

Proofs of the results in Section 3
In this section, we prove our main results Lemma 3.2, Proposition 3.4, and Theorem 3.5.

Proof of Lemma 3.2
Proof.We set wh = P h w with P h defined in (3.2) and η h = wh − w h .Further, we introduce η ∈ H 1 (Ω) as the solution of (∇η, ∇v) for all v ∈ H 1 (Ω).
By the triangle inequality, it holds First step: Estimate of T 1 : By elliptic regularity theory and (3.11) we deduce We observe that η h satisfies with L h as defined in (3.4).Hence, η h is the finite element approximation of η.Standard finite element theory then yields where we used (3.11).Next, we re-write the equation for η and observe that it solves Hence, we obtain by (2.9) together with (4.2) and (3.11) that where we used the resolution condition k(kh) p+1 ≤ C 1 in the last step.

Proofs from Section 3.3
Proof of Proposition 3.4.First step: Iteration scheme (2.4):As before, let for simplicity u (0) h,p ≡ 0. From Lemma 3.1 we obtain that a unique solution u (1) h,p ∈ V h,p to the discrete version of (2.4) exists and moreover, that it satisfies due to (3.10) and Lemma 3.2 The latter implies that kε u min{θ 1 , θ 2 } so that the assumptions of Lemma 3.1 are satisfied.Inductively, we conclude that the whole iteration sequence (2.4) exists.Each u (j) h,p is the unique solution of a linearized Helmholtz problem and satisfies the stability estimates h,p and observe that v As we have kε u h } l∈N0 forms a strictly contracting sequence or, in other words, the iterations {u (l) h,p } l∈N0 form a Cauchy sequence.Therefore, they converge to some u h,p and it is easy to verify that u h,p is a solution to (3.1).Further, u h,p satisfies the stability estimates (3.14) as they are satisfied in all iterations (see above).Following the arguments in the estimate for v (l) h above, we can also conclude the uniqueness of u h,p .Finally, (3.13) for j = 1 follows as in the Banach fixed-point theorem.
Second step: Iteration scheme (2.6):We transfer the proof of Proposition 2.4 to the discrete setting.We show by induction that for all l ≥ 1, u (l) h,p is well-defined and satisfies where the constants in may depend on σ 2 , but not on l.
The case l = 1 directly follows from u h,p = 0 and (3.10) and Lemma 3.2.Let the statement be satisfied for l.Since kε u With the assumption σ 2 < 1 we obtain recursively Employing Lemma 3.2, we furthermore obtain As in the first step, we define v (l) h,p , which solves The a priori estimates for u (l) h,p from above yield with (3.10) h } l∈N0 forms a strictly contracting sequence or, in other words, the iterations {u (l) h,p } l∈N0 form a Cauchy sequence.Therefore, they converge to some u h,p , which one can easily identify as the unique solution to (3.1) from the first step.The stability estimates (3.14) are then already known, and, finally, (3.13) for j = 2 again follows as in the Banach fixed-point theorem.
Proof of Theorem 3.5.We proceed as in [24].Let {u (l) } l∈N0 and {u (l) h,p } l∈N0 be iteration sequences for u and u h,p , respectively.We bound the error u (l) − u (l) h,p and at the very end, let l → ∞.Since we eventually consider the limit, we simply work with the iteration sequences defined according to (2.4).The adaption to (2.6) is straightforward and yields a similar result.We define a sequence {ũ We split the error as h,p ) and estimate both terms separately.For the first term, we obtain directly from Lemma 3.1 that h,p .We observe that η From Lemmas 3.1 and 3.2, we hence obtain data is sufficiently small, we consequently have By induction together with (4.3) and η Finally, employing the triangle inequality, (4.3) and Letting j → ∞ finishes the proof.

Numerical experiments
In this section, we investigate the higher-order FEM for the nonlinear Helmholtz problem numerically.We first focus on the discretization error in dependence on the mesh size h, the polynomial degree p and the wave number k, where we aim to illustrate the theoretical findings of Theorem 3.5.In the second part of the numerical experiments, we examine the convergence of the nonlinear iteration, where we aim to analyze the dependence of the contraction factor on k, h, p, ε and the data, cf.(3.13).We also compare the two different iteration schemes presented in Section 2.2.For all experiments, we set Ω = B 1 (0) ⊂ R 2 and D = B 0.5 (0).All simulations1 were obtained with NGSolve [22,23].

Convergence of the discretization error
As data, we choose g ≡ 0 and with x 0 = (−0.55,0).Since an analytical solution is not known, we use the finite element solution on a mesh with h = 2 −7 and polynomial degree p = 3 as reference solution for the following error plots.We always depict relative errors in the • 1,k norm.Unless otherwise mentioned, we solve the nonlinear system using the frozen nonlinearity iteration (2.4) until either the (relative) residual is smaller than 5 • 10 −7 or the maximum of 20 iterations is reached.
The results for two different values of ε are depicted in Figure 5.1.Firstly, we observe that the (asymptotic) error behavior is not influenced by ε as expected by Theorem 3.5.Moreover, we confirm the expected convergence rates h p for this smooth right-hand side.Similar to and as expected from the linear case, we further observe that the plateau of error stagnation is larger   for growing wave numbers, but this pollution effect can be reduced by increasing the polynomial degree.We emphasize that we could always achieve a relative residual of at least 5 • 10 −7 in less than 20 iterations in the convergence regime.On the other hand, for very coarse meshes (when the condition k(kh) 2p 1 is not met), the fixed-point iteration may not converge.We stress that this does not contradict our theory.
To compare the error convergence for different p from another perspective, we investigate how many degrees of freedom are required to obtain a relative energy error below a certain tolerance, say, 0.06.As Figure 5.1 suggests that the behavior for the two different ε values is quite similar, we fix ε = 0.01 in the following.Table 5.1 summarizes our findings, where −− indicates that we could not achieve the desired energy error with the considered meshes.We make two important observations.First, for fixed wave number, the required number of degrees of freedom for the targeted accuracy decreases when we increase the polynomial degree.This is explained by the better convergence rate and the shorter stagnation phase.Second, for the wave number k = 8, 16, 32, we can obtain the desired accuracy with (at most) 7, 746 degrees of freedom.For this, we need to slightly increase the polynomial degree, which is especially visible if the results for k = 16 and k = 32 are compared.Note that we expect to reduce the required number of degrees of freedom for k = 64 if we considered even higher order spaces with p ≥ 4.These observations agree very well with the hp-FEM convergence analysis of the linear case where the polynomial degree should be adapted like p log k, see [18,19].

Convergence of the nonlinear iteration
We now turn to the behavior of the iteration schemes and we aim to shed light on their dependence on k, h, p, ε and the data.In our investigations, we will study (i) the number of iterations required to reach a (relative) residual of 5 • 10 −7 or (ii) the contraction factor in the l-th iteration step.First, we fix k = 8, ε = 0.1, g ≡ 0 and f ≡ 50 and investigate the h and p-dependence.We give the number of iterations till convergence (as explained above) for the iteration scheme (2.4) in Table 5.2.In the second and third column we see that the iteration scheme does not seem to be affected by p.This is in good agreement with (3.13) in Proposition 3.4, where the contraction factor does not contain p.In the other columns we compare the dependence on h for p = 1 and p = 2.Each time we refine the coarse mesh three times.The number of iterations only increases by a few steps.Since [11] shows that the | ln h| can be removed in two space dimensions, this is in good alignment with the theory.
As next step, we investigate the dependence of (2.4) on the wave number k.We fix ε, f , and g as above, set h = 2 −6 and p = 3.Table 5.3 clearly shows a decrease of the required iterations (until the residual is below 5 • 10 −7 ) with growing wave number.This indicates that the k-independent contraction factor for d = 2 in Proposition 3.4 seems to be sub-optimal.Our results would suggest that the contraction factor may even decrease like k −1 .However, we also note that in all the experiments on the h, p and k-dependence of the nonlinear iteration so far, the contraction factors σ (l) varied rather considerably over the iteration number l in each experiment.For this reason, we gave these results in terms of the required iterations.The results of Tables 5.2 and 5.3 are qualitatively the same also for the iteration (2.6) and therefore omitted here.
We are mainly interested in how the convergence of the iteration depends on the L 2 -norm of f and on the size of ε.We fix k = 16, h = 2 −5 and p = 2 and let the schemes (2.4) and (2.6) iterate until either the (relative) residual is below 5 • 10 −7 or the maximum number of 50 iterations is reached.Recall that we choose f as a constant function on the whole domain for this experiment.Figure 5.2 shows the contraction factors σ (l) for the iteration schemes (2.4) and (2.6) for different values of f .We first observe that there is an initial phase with varying σ (l) before an almost constant contraction factor is reached.This constant limit regime numerically verifies that both iteration schemes are of linear order like a fixed-point scheme.Additionally, we see that the initial phase seems to be longer for the frozen nonlinearity scheme (2.4).Comparing the behavior across different values of f , we clearly see that the contraction factor grows with larger f in accordance with Proposition 3.4.
To make this better visible, we plot the average value of σ across all iterations versus the value of f in Figure 5.3 left.For small values of f , both iteration schemes perform similarly, but for larger f , (2.6) has better contraction and convergence properties and, in that sense, is more robust.In particular, for f ≡ 150, the frozen nonlinearity (2.4) seems to be not converging (after 50 iterations, the relative residual is still 0.15), while the scheme (2.6) converges.In a similar spirit, we depict the average contraction factors over different values of ε in Figure 5.3 right.Again, we see that the frozen nonlinearity is less robust than (2.6).The dependence on ε seems in general to be more severe that the one on f .While the frozen nonlinearity converges only for ε = 0.1, 0.2.0.4 in our example, the scheme (2.6) converges for all considered values up to ε = 2.For ε = 1.6 and ε = 2, 50 iterations did not suffice to reach the residual of 5 • 10 −7 , but the contraction factors clearly suggest that (2.6) should be able to reach that tolerance if we allowed more iterations.In fact, we obtain a final residual of 5.4 • 10 −7 for ε = 1.6 and of 2.3 • 10 −6 for ε = 2.
Summarizing, we could numerically confirm the better "robustness" of scheme (2.6) over the frozen nonlinearity with respect to ε and the (right-hand side) data, cf.[26].Unfortunately, we could not reflect this in our theory, where the assumption on the data are more restrictive for the iteration (2.6).Further, our experiments indicate that the k-dependence of the contraction factor may be relaxed from the theoretical prediction in Proposition 3.4 and [24].

Conclusion
In this contribution, we studied the finite element method with arbitrary but fixed polynomial degree for the nonlinear Helmholtz equation.By employing an error analysis for a linearized Helmholtz problem with small perturbation of the wave speed, we showed well-posedness and a priori error estimates under a smallness of the data assumption and the resolution condition k(kh) p 1.In the treatment of the nonlinearity, we considered two different iteration schemes which can both be interpreted as fixed-point iterations.Our numerical experiments confirm the theoretical estimates and, moreover, indicate that the results on the hp-FEM can be transferred from the linear case as well.Additionally, we compared the two iteration schemes concerning the performance and robustness with respect to the smallness of data assumption.The contribution leaves some interesting future research questions, namely on a true hp-FEM analysis -our constants may depend on the polynomial degree presently -and on the different robustness of the two iteration schemes.

A. Proof of Proposition 2.3
Proof of Proposition 2.3.Let L Ω , L Γ , H Ω , H Γ be the high-and low-frequency filters on Ω and Γ, respectively, as introduced in [19,Sec. 4.1.1].For convenience, we briefly repeat the construction.Let F be the Fourier transform on R d and set Note that the sign for ik in the boundary conditions is flipped in comparison to [19], but the estimates remain valid.
By the linearity of (2.7), it suffices to prove Proposition 2.3 with only volume data f or boundary data g separately.We show how the assertion follows from the well-known results in case g = 0, the other case can be proven similar.We set w I By the definition of the solution operators, we deduce that the remainder r : By the estimates for w I,II H 2 and w I,II A from [19], we obtain where q can be chosen arbitrarily small by adjusting η above and C is a k-and Φ-independent constant.Here, we implicitly used that Ω is star-shaped and, hence, the stability constant of the Helmholtz equation with ε = 0 is of the order one, cf.[17,4].Clearly, we see the existence of a constant θ 2 such that, if kε Φ 2 L ∞ (D) ≤ θ 2 , we have f 0 < q f 0 for some q < 1. Iterating this argument, we can write w as sum of series (one series of analytic functions, one series of H 2 -functions) that can be bounded with the help of the geometric series.

B. Proof of Lemma 3.1
Proof of Lemma 3.1.We proceed similar to [5].Let P h be defined via (3.2) and write w − w h = w − P h w + w h − P h w = ρ + η h .We have by (3.3) and the solution splitting for w ρ where we used the approximation properties of V h,p and Proposition 2.3 in the last step.Hence, we only have to consider η h = w h − P h w in the following.Observe that η h satisfies (∇η h , ∇v h ) − (k 2 (1 + χ D ε|Φ| 2 )η h , v h ) + ik(η h , v h ) Γ = (k 2 (1 + χ D ε|Φ| 2 )(P h w − w), v h ) for all v h ∈ V h,p .First step: Insert v h = η h and consider the imaginary part.We obtain with the standard L 2 -projection Π h that k η h 2 Γ = {(k 2 (1 + χ D ε|Φ| 2 )ρ, η h )} ≤ k 2 Π h ρ 1−p,h η h p−1,h + kθ 2 ρ 0 η h 0 As in [5, p. 792], we have Π h ρ 1−p,h h p ρ 1,k .Further, (3.3) gives ρ 0 h ρ 1,k .Hence, Second step: Recall the definition of A h in (3.5).Consequently, we have for any v h ∈ V h,p .For given 1 ≤ m ≤ p, set v h = A m−1 h η h to obtain From trace inequalities (cf.[5]) and (3.7) we obtain Altogether, we have for 1 ≤ m ≤ p that and by Young's inequality we obtain As in [5], it holds that k Π h ρ m−1,h h 1−m ρ 1,k , so that we deduce Using the splitting according to Proposition 2.3 for z and the properties of P h , we have Inserting these estimates as well as ρ 0 h ρ 1,k into the one for η h , we get where we used kh 1 Hence, if k(kh) 2p ≤ C 0 sufficiently small and θ 3 1 sufficiently small, we have η h 0 (h + (kh) p ) ρ 1,k (h 2 + h(kh) p + (kh) 2p )C data by the stability of w and P h .Combining with (B.1) for m = 1, we also obtain the bound for η h 1,k , which finishes the proof of (3.9).The triangle inequality and the stability of w then also give us (3.10) as well as existence and uniqueness of the discrete solution.

Figure 5 .
Figure 5.1.:Relative error in the energy norm versus mesh size for the first experiment.Top row ε = 0.05, bottom row ε = 0.01.In each row from left to right p = 1, 2, 3.

Figure 5 . 3 .
Figure 5.3.: Average contraction factors over values of f (left) and of ε (right) for the two iteration schemes.

Table 5 .
2.: Number of "frozen nonlinearity" iterations and final residuals for different values of h and p.

Table 5 .
3.: Number of "frozen nonlinearity" iterations and final residuals for different values of k.