{'fraction_non_alphanumeric': 0.047248534055029316, 'fraction_numerical': 0.017685310479627125, 'mean_word_length': 4.0027265889432115, 'pattern_counts': {'":': 0, '<': 20, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A time-domain preconditioner for the Helmholtz equation

Time-harmonic solutions to the wave equation can be computed in the frequency or in the time domain. In the frequency domain, one solves a discretized Helmholtz equation, while in the time domain, the periodic solutions to a discretized wave equation are sought, e.g. by simulating for a long time with a time-harmonic forcing term. Disadvantages of the time-domain method are that the solutions are affected by temporal discretization errors and that the spatial discretization cannot be freely chosen, since it is inherited from the time-domain scheme. In this work we address these issues. Given an indefinite linear system satisfying certain properties, a matrix recurrence relation is constructed, such that in the limit the exact discrete solution is obtained. By iterating a large, finite number of times, an approximate solution is obtained, similarly as in a time-domain method for the Helmholtz equation. To improve the convergence, the process is used as a preconditioner for GMRES, and the time-harmonic forcing term is multiplied by a smooth window function. The construction is applied to a compact-stencil finite-difference discretization of the Helmholtz equation, for which previously no time-domain solver was available. Advantages of the resulting solver are the relative simplicity, small memory requirement and reasonable computation times.


Introduction
Time-harmonic solutions to the wave equation can be computed in the frequency or in the time domain. In the frequency domain, a discrete version of the Helmholtz equation is solved. In the time domain, the periodic solutions to a discrete wave equation with time-harmonic forcing term are sought. Frequency domain methods involve less degrees of freedom. However, the indefinite linear systems resulting from discretizing the Helmholtz equation are often difficult to solve. Time-domain methods can be attractive because they require relatively little memory and are easy to implement if a time-domain solver is available. The introduction of [1] contains a recent overview of time-harmonic wave equation solvers.
Time-domain methods are based on the correspondence between Helmholtz and wave equations. In this paper we will use as example the damped wave equation where u = u(t, x) is the wave field, f = f (t, x) is the forcing term, ∆ is the Laplacian, c(x) is the spatially dependent wavespeed, R is a spatially dependent damping coefficient, and x is in some domain Ω with Dirichlet and/or Neumann boundary conditions. Let u(x, t), f (x, t) and U (x), F (x) be related by (2) u(x, t) = e iωt U (x), and f (x, t) = e iωt F (x).
Then u satisfies a wave equation with forcing term f if and only if U satisfies the following Helmholtz equation with forcing term F (3) − ∆U − ω 2 c 2 U + iωRU = F. This equation is supplemented with Dirichlet and/or Neumann boundary conditions that carry over from those for (1).
We briefly review some time-domain approaches. The most basic time-domain method is derived from the limiting-amplitude principle. This principle states that solutions u(x, t) to (1) with zero initial conditions and forcing term (4) f (x, t) = e iωt F (x) satisfy (5) u(x, t) = e iωt U (x) + O(1), t → ∞ under certain conditions on the problem, see [4] and references therein. Thus, if u(x, t) is a (numerical) solution to this initial boundary-value problem, and T is some large time, measured in periods, then an approximate solution to the Helmholtz equation is given by (6) e −i2πT u(x, 2πω −1 T ).
We will call this the limiting-amplitude approximate solution for time T . A more advanced method is the exact controllability method [5]. In this method the periodicity of the solutions is enforced using optimization. The starting value for the optimization procedure is typically some partially converged limiting-amplitude solution. Recently more insights in and some improvements to this method were obtained [11], and its parallel implementation was studied [10]. In [1] an optimization approach called WaveHoltz was introduced. This method again uses a form of optimization but with a different optimization functional. An important feature of the methods just described is that they approximate periodic solutions of a given discrete wave equation. This has two consequences that are in general not desirable. First, the results will be negatively affected by both spatial and temporal discretization errors, while solutions to discrete Helmholtz equations only have spatial discretization errors. Secondly, the method is limited to situations in which a time-domain scheme is available, and not (directly) applicable if one only has a discretization of (3) available, or perhaps only a linear system with similar properties.
In this paper, we will address both of these shortcomings by developing a new time-domain solver for discrete Helmholtz equations. The method takes as a starting point a linear system (7) HU = F, where H is a complex N × N matrix, such that (8) Re H is symmetric, and (9) Im H is symmetric positive semidefinite.
In equation (7) F is a vector in C N and U is the unknown, also in C N . The applications we have in mind are discretized Helmholtz equations and other time-harmonic wave equations. For a certain finite difference discretization of the Helmholtz equation the method will be worked out in detail and numerical results will be given. We next describe the steps involved in the construction of the new time-domain solver. First an N × N system of second order ODE's (10) ∂u 2 ∂t 2 + Au + B ∂u ∂t = f and a frequency parameter ω are constructed. We will look for time-harmonic solutions with frequency ω of the system (10). Note that ω is in general not the physical frequency parameter used to derive (7). It is a computational parameter, chosen together with the matrices A and B, and it depends on H in a way to be specified. Time-harmonic functions u = e iωt U and f = e iωt F , with U, F ∈ C N , satisfy (10) if and only if (11) (−ω 2 I + iωB + A)U = F. Therefore we will choose A, B and ω such that (12) H = −ω 2 I + iωB + A.
The system (10) plays the role of a semi-discrete wave equation. Secondly, this system is time-discretized. This is done in such a way that time-harmonic solutions of the discrete-time system are exactly those of the continuous-time system. I.e. if u n , f n are related to U, F by u n = e iωn∆t U and f n = e iωn∆t F then u n , f n are solutions to the discrete-time system if and only if U, F satisfy (11). For this purpose we present two modified leapfrog methods. This is related to ideas from the papers [12,22], see also the optimized time-stepping method in [3].
Having a time-discretization of (10) at hand, the third step is to define a map from a right-hand side F in (7) to an approximation for the solution U . For this we follow the idea of equation (6) with one modification, which is the inclusion of a smooth window function in the time-harmonic forcing term.
It is known that the convergence of limiting-amplitude approximate solutions can be slow in some cases, e.g. in case of resonant wave cavities. Therefore we will not use the approximate solution operator directly, but use it as a preconditioner in an iterative solution method such as GMRES or BiCGSTAB (Krylov accelleration). This is the fourth and last step of our construction. The new approximate solution operator will be called a time-domain preconditioner. The term time-domain preconditioner was used before in [26]. In that work a related problem was solved, but the resulting method was substantially different. The above distinction between frequency-and time-domain methods appears no longer satisfactory for this method: It is based on time-domain methodology, but solves a frequency-domain discrete wave equation.
Conditions for stability of the new time-discretization and the convergence for large T of the approximate solutions to the exact solution of (7) will be established theoretically, including estimates for the convergence if Im H = 0. Numerical examples confirm the converge of the approximate solutions and show the usefulness of Krylov accelleration in cases where a resonant low-velocity zone is present. Results depend weakly on the large time parameter T used in the preconditioner, and a parameter for the window function, i.e. there is a large set of suitable parameter choices.
A few situations that do not fit in the classical time-domain setting, but can be handled with the new solver are as follows. First one can use it with discretizations that have been designed specifically for the Helmholtz equation. For example finite-difference discretizations that minimize dispersion errors such as those from [2] for the 2-D case and [18,21,24] for the 3-D case. Our examples are about this application. One can also imagine situations where the physical time-domain model is complicated to simulate, but reduces to a relatively simple Helmholtz equation in the frequency domain. A third use case is in the context of a multigrid method. In multigrid methods the coarsest level system still has to be solved by another (non-multigrid) method. It is based on the linear system one started with and on the choice of multigrid method. In Helmholtz equations the convergence can be very sensititve to the choice of the coarse level system and it is recommended to use certain prescribed discretizations [20,18].
However, the usefulness of the method is not restricted to these cases, and some of the ideas could also be incorporated into other time-domain methods.
The contents of the remainder of the paper is as follows. In section 2, the construction of the new solver is described. In section 3 examples of this construction are given in case H results from certain finite-difference discretizations. Section 4 contains the theoretical results. After that, section 5 contains the numerical examples. So far it was assumed that Im H is diagonal, in order for the time stepping method to be explicit. Section 6 describes a variant of the method that involves an explicit time stepping method even if Im H is nondiagonal. We conclude the main text with a discussion section. Appendices contain some remarks on the compact-stencil finite-difference discretization described in section 3 and the derivation of a Fourier transform used in section 4.

Method
In this section we describe in detail the construction of a time-domain preconditioner for a matrix H satisfying (8) and (9). We recall from the introduction that there are three main steps: (i) the definition of a suitable second order system of ODE's of the form (10); (ii) the definition of a suitable time-integration method for this system of ODE's; (iii) the definition of a linear map that produces approximate solutions based on the limiting-amplitude principle.
Step (ii), the time-integration method, will be discussed first, since the properties of the time-integration method affect the choice of the system of ODE's (10). Then steps (i) and (iii) and the application of the method as a preconditioner are discussed.
2.1. Frequency-adapted time discretizations of (10). The leapfrog or basic Verlet method is a standard method to integrate equations of the form (10) in case that B = 0. It is obtained, basically, by discretizing the second order time derivative using standard second order finite differences. To allow for nonzero B, the damping term has to be discretized as well. A standard way to do this is with central differences [6,17,16]. This yields the equation from which u n+1 can be solved. Note that u n denotes the discrete approximation to u(n∆t), and that u(t) ∈ C N . This leads to an explicit method only if B is diagonal. In section 6 we will discuss a variant that results in an explicit method in case B is non-diagonal.
We will formally define the time-integrator resulting from (13).
Central differences damped leapfrog will be defined as the time integrator given by (15) u n+1 = I cd (u n , u n−1 , f n ) := I + 1 2 L −1 (2 − K)u n − (I − 1 2 L)u n−1 + g n . The stability of these methods is discussed in section 4. According to Theorem 1, I cd is stable if (16) K and L are positive semidefinite, and (17) 4I − K is positive definite.
The latter condition leads to a CFL bound, that will be discussed below.
As mentioned in the introduction, we look for discretizations such that the time-harmonic solutions of frequency ω of the discrete-time system are exactly those of the continuous time system (10). Due to discretization errors this is not the case for I cd . In the next proposition we will show that the time-harmonic solutions to (10) satisfy a recursion of the same form as (13), but with different choices of A, B. From these recursions modified schemes can derived that have the desired property. Proposition 1. Let u n and f n be related to U, F ∈ C N by (18) u n = e iωn∆t U, f n = e iωn∆t F and let and β = ω∆t sin(ω∆t) .
Then U, F satisfy (11) if and only if u n , f n satisfy Proof. To prove the first claim,Ã,B andc will be constructed such that if and only if (11). Inserting u n = U e inω∆t into (22), results in Using the definitions of α and β to rewrite the left-hand side, this is equivalent to Multiplying by α results in This is equivalent to (11) ifc = α −1 andÃ andB are defined as in (21).
Based on the proposition we define the following time integrators. The time-harmonic solutions with frequency ω of these integrators correspond exactly to time-harmonic solutions of (10).
Definition 2. Frequency adapted central differences damped leapfrog will be defined as the time integrator given by (26) u n+1 = I acd (u n , u n−1 , f n ) := I + 1 2 L −1 (2 − K)u n − (I − 1 2 L)u n−1 + g n , where K, L and g n are given by The time integrators I cd and I acd are of the same form with different choices for K and L. Therefore for I acd the stability conditions are again (16) and (17), but now with K, L as in (27).

2.2.
Choice of the semi-discrete system and the parameters ω and ∆t. We now look for a second order system of ODE's of the form (10), and a parameter ω such that the time-harmonic solutions to (10) satisfy HU = F . Recall that ω is the computational frequency parameter, that is chosen in the construction of the algorithm, and not the physical frequency used in the underlying Helmholtz equation. We assume the time-integrator I acd is used and will discuss the choice of the parameter ∆t as well. An explicit and stable scheme is obtained when Im H is diagonal.
Section 6 treats the case that Im H is non-diagonal. Because of (12), we set where ω is still to be determined.
We assume the integrator I acd is used, and discuss under which conditions this leads to stable and explicit scheme. Clearly, the resulting discrete time system is explicit if Im H is diagonal, and generally not otherwise, cf. (29) and definition 2. The stability conditions are given in (16) and (17). Because Im H is positive semidefinite, L is also positive semidefinite. To ensure that K is positive semidefinite, we set (30) ω = −λ min (Re H).
(or to a lower bound for −λ min (Re H) if this value is not exactly known). The other stability condition (17) is a form of the well-known CFL condition. It implies that the eigenvalues of K must be less than 4. For central differences damped leapfrog integration (not frequency adapted) it implies the condition For the frequency adapted variant it implies, by (27) and (19) the condition sin( ω∆t ).
The choices and requirements (29), (30) and (32) define a suitable choice of parameters of the systems of ODEs and the I acd integrator in case Im H is diagonal.

Time-domain approximate solution operators.
In this section the time-domain approximate solution operator will be defined. We also show that a complex approximate solution can be computed by solving a real time-domain wave equation. This appears to be a standard trick in the field. The main novelty is that the formula for the time-harmonic forcing term (4) is modified so that the forcing is turned on gradually.
We first specify the window functions that will be used in the forcing term.
Definition 3. A function χ : R → R is called a C k admissible window function if the following requirements are satisfied.
The time-domain approximate solution operator associated with the integrator I acd will be denoted by S acd T and is defined as follows.
Definition 4. Let A, B, ω and ∆t satisfy the requirements of subsection 2.2. Let χ be an admissible C ∞ window function and let T be a positive real constant, such that n steps := 2πω −1 T /∆t is an integer. For F ∈ C N , let The time-domain approximate solution operator for H associated with the integrator I acd is the linear map S acd T : C N → C N defined by (34) S acd T F = e −i2πT u nsteps , where u n , n = 0, 1, . . . , n steps is given by In section 4 the convergence of S acd T F to H −1 F will be established under the assumption that requirements for stability of the the time integrators discussed in subsection 2.2 are satisfied. Section 5 contains numerical examples for S acd T . There is still a lot of freedom to choose the window function χ. In the numerical examples we will introduce therefore an additional parameter ρ, 0 < ρ ≤ 1, such that the window TimeDomainPreconditioner(F, K, L, T, n P ) 1 u −1 = u 0 = 0 2 for j = 0, . . . , T n P − 1 Table 1. Algorithm for a time-domain preconditioner. Computed is U = S acd T F , with U, F ∈ C N . The number of time steps per period is denoted by n P .
function is 1 on [0, (1 − ρ)] and positive but strictly less than one on (1 − ρ, 1). On the interval (1 − ρ, 1) a sine square is used, which leads to the following form of χ The parameter ρ will be called the window parameter. Based on the convergence analysis in subsections 4.3 and 4.4, ρ should not only be strictly larger than 0 but also strictly smaller than 1. We next show that it is sufficient to solve a real time-domain wave problem to compute the limiting-amplitude approximate solution. The argument is given in the continuous case, but is applicable equally well in the discrete case.
Let F be a complex right-hand side for the Helmholtz equation and u(t, x) be the solution to (1) with right-hand side f (t, x) = F (x)e iωt . Assuming the limiting-amplitude principle holds, cf. (5), an approximate solution to the Helmholtz equation is given by for some large t.
The field Re u(t, x) can be determined by solving the real wave equation with real forcing term, i.e. with forcing term Approximations to Re U (x) and Im U (x) can hence be obtained from Re u(t, x) by 2π ω with T a large integer. Therefore real time-domain simulation is sufficient. Table 1 contains an algorithm for computing S acd T U using time stepping with real fields. The number of time steps per period is denoted n P . It is assumed n P is a multiple of 4 to facilitate the computation of Re U and Im U as in (40). If n P is not a multiple of 4 one could alternatively take u j for different choices of j, with relative phase not equal to 0 or 180 degrees, and extract Re U (x) and Im U (x) by solving a 2 × 2 linear system.

2.4.
Time-domain preconditioned GMRES. The approximate Helmholtz solver can be used as a preconditioner for iterative methods like GMRES of BiCGSTAB. Without preconditioning, the system to be solved is HU = F . Applying left-preconditioning means that instead the system The vector V is solved from this system and the solution to the original problem is then given by U = P V . We propose a solution method where GMRES is applied to the leftpreconditioned system.

Examples
The examples we consider are finite-difference discretizations of the damped Helmholtz equation HereR is the spatially dependent damping in units of damping per cycle. We start with standard second order differences and then apply the method to the discretization from [18]. The method is not limited to finite-difference discretizations. Finite-element discretizations can also lead to linear systems HU = F with H satisfying (8) and (9).
We will see that in our method, the discrete system is a discretization of a modified PDE, see (48) below, and that the CFL bound for this modified PDE is independent of the velocity.
3.1. Second order finite differences. It is instructive to start with a simple second order finite-difference discretization. For readability we describe the two-dimensional case. The degrees of freedom will be denoted by where h denotes the grid spacing, and Dirichlet boundary conditions are assumed, i.e. U (i,j) = 0 if (i, j) / ∈ D. In this case Im H is diagonal, and ω and ∆t are chosen based on the values We find that ω = k max while ∆t is chosen from (32) using that λ max (A) = 8h −2 +k 2 max −k 2 min . The scheme in the computational time domain, can be written as This follows from (20), by multiplying this equation with α, and using that A = Re H + ω 2 , B = ω −1 Im H and that in this case ω = k max . The scheme (47) differs substantially from the standard second order FDTD scheme. In fact, it is a discretization of the PDE The highest order part of this PDE contains a wave operator ∂ 2 ∂t 2 − ∆ with constant velocity, even if the velocity in the Helmholtz equation we started with is non-constant (this velocity is proportional to 1/k(x)). Thus we have identified another difference with alternative timedomain methods like those of [5,11,1], that use a discretization of the standard wave equation.
On the other hand, this difference can also be largely removed. If before applying our method the operator H would be rescaled, multiplying from the left and the right by a diagonal matrix with entries c (i,j) on the diagonal, then the resulting computational time domain scheme would be close to a standard FDTD scheme.

3.2.
An optimized finite-difference method. Next we consider a discretization with a 27 point cubic stencil (in 3-D) that was described in [18]. There exists different variants of such compact stencil methods, see among others [2,21,24]. For each gridpoint, the number of neighbors that the gridpoint interacts with is relatively small (when compared e.g. to higher order finite elements). Even so, the dispersion errors of these schemes are also relatively small, as shown in [18], so that these methods can be used with relatively coarse meshes 1 .
The equation discretized in [18] was the real part of (43). We will discuss here the situation with constant k, for variable k we refer to appendix A. The stencil is the 27 point cube that we will denote with {−1, 0, 1} 3 . For constant k, due to symmetry there are four different matrix coefficients. The values of these coefficients are h −2 f s ( hk 2π ), for s = 0, 1, 2, 3 respectively, where the dimensionless function f s describes the coefficient as a function of the dimensionless quantity hk 2π . To be precise, In [18] the f s are precomputed functions, in a class of piecewise polynomial functions, obtained by optimization to minimize phase errors. For the 2-D case, in [2] explicit expressions for the optimal functions f s , s = 0, 1, 2 were derived. The imaginary part was discretized as above, i.e. it was diagonal with entries (k (i,j,l) ) 2R(i,j,l) π . Because Im H is diagonal, the integrator I acd is used. The discretization in [18] has in addition the property that it was second order and such that (50) Re H + k 2 I is symmetric positive definite.
Thus ω = k in case of constant k. The upperbound for Re H that determines the maximum for ∆t is given by This bound is typically somewhat less than the bound associated with second order finite differences with the same parameters. It is straightforward to follow the recipe of subsection 2.2. This scheme can again be considered a discretization of (48) rather than of (1).

Analysis
Here theoretical results will be obtained concerning the stability of the time integrators defined in subsection 2.1, the solutions of these time-integration schemes and and the convergence of the approximate solution operators defined in subsection 2.3.
The consequences of the stability analysis for the choice of ∆t were already discussed in subsection 2.2.

Stability.
To study the conditions under which (15) and (26) are stable, we study the growth of the solutions to the homogeneous recursion (52) (1 + 1 2 L)u n+1 + (−2 + K)u n + (I − 1 2 L)u n−1 = 0. Throughout we will assume that K and L are symmetric matrices and that L is positive semidefinite. The basic idea is to show that the following energy functions stays bounded Here ·, · denotes the standard inner product. If the energy function is coercive, than the solution also stays bounded. The following theorem, states the stability conditions that follow from such an analysis. The conditions under (ii) in the theorem were already mentioned above, see equations (16), (17). This result is not new, see [9], page 8, for an earlier discussion of the same energy function. Proof. To prove the result, a bound is derived for Using that K is symmetric in combination with basis rules for standard products results in By using the recursion equation in two places results one obtains the estimate where the last inequality holds because by assumption L is positive semidefinite. The theorem follows from this result.

Eigenvalue analysis.
There are two cases for the eigenvalue analysis of the solutions to (52), depending on whether K and L have joint eigenvectors. The standard case in which they do is worked out because it is used subsection 4.4. We also analyse the other case. In the presence of damping (Im H = 0, hence L = 0) the intuition is that undamped solutions must correspond to vectors in the kernel of L, and hence to eigenvalues of K, in which case it is clear that an undamped solution can exist. The purpose of our analysis was to confirm that there are no other undamped solutions, because such solutions could affect convergence (even though for the convergence proof the result is not needed).
So let v be a joint eigenvector of K and L with eigenvalues k and . Solutions to equation (52) of the form u n = Λ n v exist with Λ given by If k < 0 or k > 4, these solutions are exponentially growing. Thus, for stability, K and 4I − K need to be positive semidefinite. In addition we require 4I − K to be strictly positive definite, to avoid a non-decaying solution of the form (−1) n v while k = 4, > 0. Next we assume that K and L do, in general, not have joint eigenvalues. We study the presence of linearly growing solutions and of solutions that neither grow nor decay using the first order form of the recursion (52), that is given by The properties of the solutions are directly related to the eigenvalues, (generalized) eigenvectors and the Jordan blocks of Ξ. We next list some properties of Ξ that relate to eigenvalues ξ of Ξ with |ξ| = 1. These correspond to undamped solutions. Claim (iii) follows from the energy growth equations. Solution to second order recursion would be of the form u n = ξ n w + nξ n−1 v, and the dominant term in the energy in the limit n → ∞ would be n(1 − ξ)v, (4I − K)n(1 − ξ)v . This term would be unbounded as n → ∞ which is impossible, hence any Jordan block of size > 1 must be associated with an eigenvalue ξ = 1.
In the situation of claim (iv), solutions are of the form u n = ξ n v with v ∈ C N . Energy must be conserved hence u n+1 − u n−1 , L(u n+1 − u n−1 ) = 0 hence (ξ − ξ −1 )u ∈ ker L, hence u ∈ ker L. But then it follows that Ku = (−ξ −1 + 2 − ξ)u, so that u is a simultaneous eigenvector of K, and L.
The first part of part (v) follows from the bound on the energy. Jordan blocks of size > 2 would lead to solutions of the form The dominant contribution to the term u n − u n−1 , (4I − K)(u n − u n−1 ) would be 2(n − 1)v, (4I − K)2(n − 1)v , which would be unbounded, which is not possible. The second part of (v) follows from the form of Ξ. If eigenvector is We first give our conventions regarding the Fourier transform, and discuss the Fourier transform of functions defined on ∆tZ. Let g(n) be a function Z. Considering n as a position variable, we have the following Fourier transform / inverse Fourier transform pair, denoting the frequency by ν Here ν in the periodic interval [−π, π] which we will also denote by T 2π . We will also denote the Fourier transform of a function g by Fg.
When considering discrete approximations of continuous quantities, it can be convenient to write u(n∆t) instead of u n and f (n∆t) instead of f n . In other words, u and f are considered as functions u(t), f (t) with t ∈ ∆tZ, and not as two-sided sequences. The forward and inverse Fourier transforms between functions on ∆tZ and functions on T 2π/∆t are defined by For functions of t ∈ ∆tZ the Dirac delta function will be defined including a factor ∆t −1 , i.e.
For function g(t), h(t), t ∈ ∆Z, we have where * denotes convolution on T 2π/∆t . In this subsection we will treat u and f as functions on ∆tZ. Associated with the integrator I acd is then a second order difference operator C acd , defined as where K, L are as defined in (27). Solutions to the time-integration satisfy The left-hand side can be see as the convolution of u with a kernel Γ acd that is a matrix valued function on ∆tZ. (The explicit expression for Γ acd is easily derived from (66).) In the Fourier domain the convolution becomes a multiplication and we have From proposition 1 it follows that In addition we know that Γ is C ∞ , because Γ is nonzero only for t ∈ {−∆t, 0, ∆t}. By Φ acd (t), t ∈ ∆tZ we will denote the causal Green's function associated with (67). It is a function ∆tZ → C N ×N defined as the solution to where δ(t) is as defined in (64). Assuming that the conditions in (16), (17) hold, Φ acd grows at most polynomially if t → ∞ (by Theorem 1), and its Fourier transform Φ acd is well defined as a matrix valued distribution on the circle T 2π/∆t .
We will first show that the approximate solution operator S acd T can be expressed as a weighted mean of Φ(τ ) around τ = ω. For this purpose we first define a helper function ψ 1/T , depending on χ and various parameters. Let θ T (t), t ∈ ∆tZ be given by then ψ 1/T is defined by We have the following result.
Theorem 3. Assume ψ and Φ acd are as just defined, then In words, s acd T equals the convolution product of ψ 1/T and Φ acd , evaluated at τ = ω. The convolution product is taken on the torus T 2π/∆t .
Proof. LetT = 2πω −1 T . By linearity, the time-domain approximate solution S T (F ) satisfies In this formula we easily recognize a Fourier transform The equality (73) follows because the Fourier transform of the product is 1 2π times the convolution product of the Fourier transforms, see (65).
We continue by studying the family of functions ψ 1/T (τ ). Lemma 1. LetT = 2πω −1 T . Suppose χ is a C ∞ admissible window function and ψ 1/T is defined as in (72). The function ψ 1/T (τ ) satisfies Proof. The Fourier transform of the function s → χ(s/T ), s ∈ R is given by σ →T χ(σT ), σ ∈ R. When s is restricted to ∆tZ then the Fourier transform is a function on the torus T 2π/∆t given by

This shows (76).
Because χ is C ∞ and compactly supported, its Fourier transform χ(σ), defined for σ ∈ R, satisfies for any m ≥ 0 and N ≥ 0. Furthermore, 1 2π χ(σ) has integral 1. As a consequence ψ 1/ is an approximation to the identity (this is actually proven in the upcoming proof). This fact and Theorem 3 imply the convergence of the time-domain approximate solutions to the correct solutions.  (16) and (17) then Proof. The Fourier transform Φ acd satisfies (81) Γ acd (τ ) Φ acd (τ ) = I.
Note that Γ acd (τ ) is a C ∞ function and Φ acd is a distribution and that this equality is in the sense of distributions. By the assumption (79) and the property (69), the matrix valued function Γ acd (τ ) is non-singular on a neighborhood of τ = ω. It follows that Φ acd (τ ) is continuous on a neighborhood of τ = ω and that We define a periodic cutoff θ per . This will be supported in [(−π − )/∆t, (π + )/∆t] for some < π and such that Next we insert the sum, and distinguish between the cases j = 0 and j = 0. The case j = 0 is further split into the set near τ = 0 and away from τ = 0 using a cutoff function θ 1 that is supported on a neighborhood of τ = 0 where Φ(ω − τ ) is smooth and is equal to one on a smaller neighborhood of 0. We thus have three terms In the first term on the right hand side, the norms θ per (τ ) χ(T (τ + j2π/∆t)) C k decay as T −N and j −N for any N , because of (78). Since Φ acd is a distribution of finite order, this term goes to zero asT → ∞. The second term goes to zero similarly. For the third term, we may assume that θ per = 1 on the support of θ 1 (τ ). We may expand where g 1 (τ ) is bounded by a constant M (actually it is bounded by the first derivative of Φ acd (ω − τ ) on the support of θ 1 ). We write the third term and then apply a change of variables as follows The first term on the right hand side converges to Φ acd (ω), since (the latter equality by the inverse Fourier transform), while the second term goes to zero, since it can be bounded by where the latter integral is finite by (78). This concludes the proof.

4.4.
Speed of convergence. The speed of convergence is the topic of this and the next subsection. We primarily use a spectral decomposition and Theorem 3 and restrict to the case that Im H = 0. In this case Re H has an orthogonal eigenvalue decomposition, so that an essentially scalar analysis can be used.
After applying a change of basis, H may be assumed diagonal, given by Matrices such as K, Φ acd (t) and Φ acd (τ ) are also diagonal. The diagonal components can be obtained by applying our approach to the scalar problem in which H = λ, where we must take into account that ω and ∆t are chosen using using the set of eigenvalues {λ 1 , . . . , λ N } of H. So, for example there is the requirement for ω. In the analysis of the scalar problem we will obtain estimates uniformly in λ. With these considerations it suffices to analyse the scalar problem.
For the quantities of this scalar problem we use the following notations. We write k(λ) for the matrix K with the above scalar problem, and φ acd (t, λ) for the Φ acd (t). For example, according to definition equations (29) and (27) the function k(λ) is given simply by We will study the decay of the error (93) error = s acd T (λ) − λ −1 , as a function of λ and T , in the limit T → ∞. We wil take into account the possibility that |λ| is small, obtaining appropriate uniform estimates.
The first step is an explicit construction of φ acd (t, λ), which amounts to solving the recursion given by (59), (58) for the case that L = 0 and K is scalar. Proposition 2. Let and U (n) be the unit step function of (173). In case 0 < k < 4 Proof. The second order recursion in this case becomes This can be written as a first order recursion It follows that Assuming 0 < k(λ) < 4, the matrix Ξ has two eigenvalues, given by, cf. (57) The decomposition of Ξ is Associated with Λ + is the following frequency parameter This implies that Λ n ± = e ±in∆tΩ(λ) . The Fourier transform φ acd (τ, λ) can now be obtained from the Fourier transform of U (n) given in Proposition 5 in Appendix B, and the rule that modulation of a function corresponds to translation of its Fourier transform. This leads to the following Proposition 3. If 0 < k(λ) < 4 and Ω(λ) is as in (102) then (103) φ acd (τ, λ) = 1 2iα sin(Ω∆t) U (∆t(τ − Ω(λ))) − U (∆t(τ + Ω(λ))) .
If k(λ) = 0 then We next state and prove a result concerning the convergence speed. Multiple bounds are provided, each uniform in λ over a different subsets of the domain containing λ. This is related to special behavior that occurs for certain values of λ. A first special case is when λ is near −ω 2 which corresponds to the case that k(λ) is near 0. In this case the factor 1 sin(Ω∆t) becomes large, or, in the limiting case that k(λ) = 0, the formula for the φ acd is of a special form involving U instead of U (note that U is more singular). A second special case is when λ is near zero. In this case ω is near Ω(λ). The closer ω is to the singularity at Ω(λ), the narrower ψ 1/T must be to have a good approximation, so the larger T should be. In each case we assume k(λ) < 4 uniformly, which implies that The assumption that χ is C ∞ is replaced by different, more specific assumptions on χ, so that the effect of the precise properties of χ becomes clearer.
Theorem 5. Instead of the assumption that χ is C ∞ , assume that χ is a C 0 admissible weight function that satisfies the following additional assumptions: First there is an integer N 1 ≥ 3 such that and secondly there is an integer N 2 ≥ 1 such that (107) χ(σ)σ N 2 dσ < ∞ and the moments of χ of order 1, . . . , N 2 − 1 are zero.
In the remainder of this subsection we will prove this result. In subsection 4.5 we will discuss the result.
Proof of Theorem 5. By rescaling problem quantities, we may assume that ∆t = 1. In terms of original problem quantities this amounts to a change of variables (111) ν = ∆tτ.
In this notation, case (I) is the case that 0 ≤ ν 1 < constant < ν 2 , and case (II) is the case that 0 < constant < ν 1 < constant < π. For brevity we will also omit the subscript from φ acd . It is given by The quantity s acd T (λ) satisfies, where we omit the superscript, We first consider case (I). To be shown is that In this case, to handle the factor 1 ν 1 around ν 1 = 0, we write φ(ν) as which is valid for k(λ) = 0 and for k(λ) = 0. The integral over s can be substituted on both sides of the equality sign in (115) and it is sufficient to obtain an estimate for the integrand. Letν = sν 1 , s ∈ [−1, 1]. It is sufficient to show that T → ∞ uniformly inν. We insert the expression for U from Proposition 5, and separately consider the contributions from the delta function and the cotangent function. Let (118) then the left hand side of integral (117) equals πI 1 − i 2 I 2 , and it is sufficient to show that (119) and Using the rules for differentiation and integration of delta functions and Lemma 1, I 1 can be written as (120) The first part of (119) now follows from Lemma 2.
Using Lemma 1 and a periodic cutoff function as in the proof of Theorem 4, the integral I 2 is rewritten as We split the integration domain in three regions for j = 0 and in four regions for j = 0. In both cases, region 2 is of width ∼ 1/T aroundν, the associated cutoff function θ 2 is equal to 1 one a region of this width, and vanishes outside a large region of this width. Regions 1 and 3 are below and above region 2 respectively. In case of j = 0, an region of width ∼ 1 around ν = ν 2 , staying ∼ 1 away from −ν is split off from region 3, with a cutoff function θ 4 . Will absorb θ per in these cutoff functions. For region and cutoff function k and summation index j there are then the contributions (122) I 2,j,k =T 2π R χ(T (ν 2 − ν))θ k (ν + j2π) cot ( 1 2 (ν +ν)) dν.
where the integral stays ∼ 1/T away from the singularity of cot and, for j = 0 it stays ∼ 1 away from ν = ν 2 where χ is large. By these properties and Lemma 2 equation (123) is true for region 3. In a similar way it is true for region 1.
We next consider the contribution for case (I) around ν = ν 2 . To show this we need assumption (107) and extend the reasoning around (86) to higher order. Let g(ν) = cot ( 1 2 (ν +ν))θ 4 (ν). This function can be expanded where h N 2 is bounded by a constant times the N 2 -th derivative of g. Then the relevant integral is transformed with h N 2 bounded by a bound for d N 2 g dν N 2 (ν). By assumption (107) we find that This concludes the proof of (108). We next consider case (II). In this case we assume T 1/|λ|. It is sufficient to show that We note that |ν 2 − ν 1 | ∼ |λ|, while ν 2 + ν 1 ∼ 1. For the + sign, in fact O(T −N (II) ) can be obtained. We will treat the − sign, the + sign can be done in the same way with minor modifications. Using the expression for U in Proposition 5 the integral on the left can be written as a sum πI 3 − i 2 I 4 where I 3 and I 4 are defined by (130) and it is sufficient to show that By Lemma 2 the first part of (131) follows. The integral I 4 can be written as We split the integration domain in three regions for j = 0 and in four regions for j = 0. In both cases, region 2 is of width ∼ 1/T aroundν, the associated cutoff function θ 2 is equal to 1 one a region of this width, and vanishes outside a large region of this width. Regions 1 and 3 are below and above region 2 respectively. In case of j = 0, a region of with ∼ |λ| around ν = ν 2 , staying ∼ |λ| away from ν 1 is split off from region 1 or region 3, with a cutoff function θ 4 . Will absorb θ per in these cutoff functions. For region and cutoff function k and summation index j there are then the contributions (134) I 4,j,k =T 2π R χ(T (ν 2 − ν))θ k (ν + j2π) cot( 1 2 (ν − ν 1 )) dν.
and I 4 = I 4,0,4 + j∈Z 3 k=1 I 4,j,k . We claim that for regions k = 1, 2, 3 For region 2 we observe that the action of cot on a test function supported on [−w, w] can be bounded by The estimate (135) follows from this fact and the assumptions on χ. For region 3 there is a contribution of the form (137) The support of the integrand stays ∼ 1/T away from the singularity of the cotangent function and, in case j = 0 it stays ∼ |λ| away from the point ν = ν 2 which is where χ is large. This implies (135) for region 3. Region 1 is treated similarly. The treatment in equations (126) and (127) is used with minor modifications to show that (138) This concludes the proof of the theorem.
4.5. Discussion. We briefly discuss Theorem 5 and its assumptions. By elementary properties of Fourier transforms assumption (106) is directly related to the smoothness of χ. For example, for a piecewise C ∞ function, with finite jumps in the k-th derivatives, such an estimate holds with N 1 = k + 1. The k-th moment of the Fourier transform χ is proportional to χ (k) (0). The second requirement of (107) can hence be satisfied by assuming χ = 1 on neighborhood of 0. The first part of (107) is again related to the smoothness of χ.
It is clear that the inversion problem becomes more difficult if the smallest absolute eigenvalue of H is closer to zero. Define the notion of "gap" by (139) gap = min{|λ 1 |, . . . , |λ N |}.
When the matrix S acd T is applied as a preconditioner, the convergence factor is of interest. In the context of subsection 4.4 this is given by The theorem leads to the following result. Suppose one wants to obtain a convergence factor ρ c < R c , with R c some fixed constant < 1. There is a constant C such that one can set which is clearly less attractive than the bound (142).

Numerical examples
In this section we study some examples. The examples all involve an optimized finitedifference discretization described in subsection 3.2 in two and three dimensions.
To study the 2-D method, a simple implementation in Julia was made. In the 3-D case the time stepping was done in C, the computation of coefficients and the GMRES iteration were done in Julia. In the C implementation, for each grid point and timestep, 14 real coefficients had to be read, of which one also had to be written back to. The C implementation used AVX extensions but was otherwise straightforward. In most computations double precision numbers were used. The exception was the time-stepping performed in C. We found that the time-domain preconditioner could also be run in single precision, since the GMRES solver would automatically take care of rounding errors in subsequent iterations.  figure 1 for non-constant velocity models and some solutions in these models.
The first objective was to the study the convergence of the approximate solution S acd T F to H −1 F . For this we took models of size 320 × 320 and 640 × 640 (excluding damping layers), and let the right-hand side be a point source. For the constant model the parameter k = ω ph c equalled k = 335.1 for size 320 × 320 and k = 670.1 for size 640 × 640. The "taper" parameter ρ was varied, we considered values in {0, 0.25, 0.5, 0.75}. The value 0 was not exactly zero, in this case the time-harmonic right-hand side grew from 0 to 1 over one period. The convergence (relative difference between approximate and true solutions) is given in Figure 2. Note the differing axes in the Figure 2(b). The main conclusions from these figures are that the convergence in the resonant model is relatively poor, and that the choice of ρ = 0 (no tapering) leads to poor convergence. The detailed amount of tapering is less important. The lower value (ρ = 0.25) appears to perform best unless extremely small relative errors are required.
GMRES convergence for the three velocity models of size 640 gridpoints (excluding damping layers) are given in Figure 3. In each case the number of periods of the preconditioner was varied. In Figures 3(a-c) the taper parameter was 0.25 and both the error indicator from GMRES and the true error are plotted. In Figure 3(d) the GMRES error was plotted and the taper parameter was varied. On the x axis, the number of iterations times the number of periods per iteration was displayed. This is roughly, but not quite proportional to the cost, the main difference being that one more preconditioner application is needed for the righthand side in the preconditioned system. In velocity models 1 and 3, there was no speedup compared to the direct application of the time-domain solver. Neither was the method much slower. However, in the resonant velocity model (velocity model 2), the performance using GMRES was substantially better than when the time-domain solver was used directly.
In the examples about 50 to 100 periods for the preconditioner worked best. This is approximately the size of the example in wavelenghts.
It is confirmed here that, in case of resonances, long time simulation (i.e. applying the limiting amplitude principle directly) leads to poor convergence. For our time-domain preconditioner, a substantial speedup can be obtained by using Krylov accelleration and a smooth window function. Exact controllability is another way of obtaining a speedup, cf.   We also studied the performance of the algorithm in three dimensions. In this case, the velocity model was the SEG/EAGE Salt Model 2 . In this case our goal was to obtain a first estimate of the computational cost of the method. A minimum of 6 gridpoints per wavelength was used. The Julia/C code was run on a 2019 MacBook Pro with a 2.6 GHz 6-core Intel i7 processor and 16 GB main memory. A GMRES error reduction with a factor 10 −5 was required. The results are in Table 5.2.

5.3.
Comparison with diagonally preconditioned GMRES. In Table 5.3 a comparison between the time-domain preconditioner and diagonally preconditioned GMRES is made. In both cases the number of matrix applications, the computation time and the error compared to the true solution is displayed. The computations were done in Julia.     Table 3. Comparison of time-domain preconditioner with a diagonal preconditioner Diagonally preconditioned GMRES is consistently slower. This is because per matrix application, the diagonally preconditioned GMRES method involves more work, since the vectors spanning the Krylov subspace must be manipulated. In specific scenarios there are more advantages of the time-domain preconditioner. When GPUs are used it is an advantage that the Krylov subspace is manipulated infrequently, which means that the Krylov subspace doesn't have to be stored in the scarce device memory. In inner-outer iterative methods the time-domain preconditioner could be used as an inner method, with the advantage that it is linear.

Explicit scheme for non-diagonal damping
If the matrix B in (13) is non-diagonal, the associated time-stepping method becomes implicit, which is not practical for a hyperbolic time-dependent system. Here we introduce a time-integrator that is explicit even in the case of non-diagonal damping by following  6.1. Definition of the method. The starting point of our frequency adapted backward differences damped leapfrog scheme is a straightforward modification of (13) in which backward differences are used to discretize the damping term (144) 1 ∆t 2 (u n+1 − 2u n + u n−1 ) + B ∆t (u n − u n−1 ) + Au n = f n .
The associated time integrator will be called backward differences damped leapfrog.
Definition 5. Let K, L and g n be as in (14). Backward differences damped leapfrog will be defined as the time integrator given by By choosing K, L differently, the resulting time integrator can be made exact for timeharmonic signals of frequency ω. This is shown in the following Proposition, which is an equivalent of Proposition 1. The associated time integration method will be denoted by I abd .
Proposition 4. Let u n and f n be related to U, F ∈ C N by (146) u n = e iωn∆t U, f n = e iωn∆t F and let α and β be as defined in (19). Then U, F satisfy (11) if and only if u n , f n satisfy Proof. To prove this claim,Â,B andĉ will be constructed such that if and only if (11). Inserting u n = U e inω∆t into (149), results in Using the definitions of α and β and multiplying by α results in the equivalent equation This is equivalent to (11) ifĉ = α −1 andÂ andB are defined as in (148).
Definition 6. Frequency adapted backward differences damped leapfrog will be defined as the time integrator given by where K, L and g n are given by We proceed by discussing the choice of A, B, ω and ∆t. Obviously we still have (29). The requirements for ω and ∆t should follow from stability conditions for I abd , cf. subsection 2.2. These stability conditions are (16) and (154) 4I − K − 2L is positive definite, as shown in subsection 6.2 below. The method for choosing ω and ∆t is somewhat more complicated than in subsection 2.2, because the stability requirement leads to two conditions that both involve ω and ∆t. It is convenient to use ω∆t and ∆t as parameters instead of ω and ∆t. The parameter ω∆t should be between 0 and π. Given ω∆t the following expression for ω 2 can be derived from the condition that K is positive semidefinite and equation (153) (155) ω 2 = −λ min (Re H) + β ω∆t (1 − cos(ω∆t))λ max (Im H) (instead of λ min (Re H) and λ max (Im H) lower and upper bounds can be used respectively). From the condition that 4I − K − 2L is positive definite we then get the following scalar condition (156) (ω∆t) 2 α λ max (Re H) ω 2 + 1 + β α (1 + cos(ω∆t)) λ max (Im H) ω 2 < 4.
From here a value of ω∆t can be obtained that satisfies the conditions by solving a simple quadratic equation. If a larger value of ω∆t is desired, one can look numerically for a value as large as possible for which (156) is still satisfied.
The time-domain approximate solution operator is defined similarly as in Definition 4.
Definition 7. Let χ be an admissible C ∞ window function and let T be a positive real constant, such that n steps := 2πω −1 T /∆t is an integer. For F ∈ C N , let

Concluding remarks
In this paper we constructed a time-domain preconditioner for indefinite linear systems that come from discretizing time-harmonic wave equations, and studied its properties and behavior analytically and with numerical examples. It should be emphasized that the method does not compute in the physical time domain. We call it a time-domain method, because of the similarity with classical time-domain methods for time-harmonic waves. For the practical application it would be useful to have further examples, with different matrices H and a comparison with alternatives. We will make a few brief remarks in this direction.
First it is clear that the method requires relatively little memory, compared to alternatives that use LU (or LDL T ) decompositions, such as domain-decomposition and direct methods, see [19,15] and references in [8]. In the 3-D implementation here, with GMRES with restart m = 10, most memory was used for the approximately m + 2 complex double precision vectors needed for GMRES. This was in part because the timestepping was done in single precision, and using real fields. In principle, memory use could be somewhat reduced by setting m = 5 or using a different iterative method like BiCGSTAB.
For finite element discretizations the cost in general will be different. These methods often have stricter CFL bounds, there may be more nonzero matrix elements, and the use of unstructured meshes may also influence efficiency. On the other hand, the size of the grid cells can be adapted to the local velocity so that the number of grid cells can be smaller.
It is difficult to compare performance of different methods, as methods are run on different computer systems, with different examples and implementations are optimized to different degrees. We refer to [1,7,13,15,19,23,25] for some alternative methods. We believe that computation times of the method as outlined are modest, and there is potential for further improvements, by using GPUs or by reorganizing the code to better make use of cached data. It remains challenging to get the most out of modern computer hardware in time-domain finite-difference and finite-element simulations. We hope that recent developments in this area, cf. [14], will lead to further improvements.