Time-domain boundary integral equation modeling of heat transmission problems

This paper investigates the numerical modeling of a time-dependent heat transmission problem by the convolution quadrature boundary element method. It introduces the latest theoretical development into the error analysis of the numerical scheme. Semigroup theory is applied to obtain stability in spatial semidiscrete scheme. Functional calculus is employed to yield convergence in the fully discrete scheme. In comparison to the traditional Laplace domain approach, we show our approach gives better estimates.


Introduction
Since the inception of the boundary integral equation method, the thermal engineering community has been exploiting its potential in solving transient heat conduction problems [47]. The method is also a popular choice among environmental scientists in the study of pollutant transport problem [26]. Recent applications include photothermal spectroscopy [38] and diffusion in variable media [1,2]. The method enjoys sustained interest due to its remarkably simple way to handle problems on infinite domains.
Various schemes have emerged to discretize time domain boundary integral equations associated to parabolic problems. The theoretical basis for much of this work originated with [4,15] and focused on the development of Galerkin discretization of the basic integral equations for problems on the exterior of a bounded domain. While we will focus on a particular class of methods (Galerkin in space, Convolution Quadrature in time), let us mention that there is a large literature on other families of schemes applied to exterior problems for the heat equation [19,43,44,33,32].
Our context is that of a Convolution Quadrature Boundary Element Method applied to transmission problems for the heat equation. Methods combining CQ and Galerkin BEM for heat diffusion problems appeared first in [29] (a combination of the ideas of CQ first set in [27], with the then recent results on the single layer operator for the heat equation) and then reinterpreted as a Rothe-type method in [14]. From the point of view of the algorithm itself, we here combine a Costabel-Stephan formulation for transmission problems [16,36], with Galerkin semidiscretization in space and multistep or multistage CQ in time [27,28,7]. We set our main goals in the proof of stability and convergence of the method avoiding Laplace domain estimates (more on this later), using instead techniques of evolutionary equations associated to infinitesimal generators of strongly continuous analytic semigroups in certain Hilbert spaces. We obtain three kinds of results: (a) long term stability of the system after semidiscretization in space; (b) optimal order of convergence in time (with explicitly expressed behavior of the bounds with respect to the time variable) for BDFbased CQ schemes; (c) reduced (but higher than stage) order of convergence for RK-based CQ schemes.
Let us next try to clarify what kind of mathematical techniques are used for our analysis and where the novelty of this work lies. To do that, we first need to discuss the mathematical background for the field. The modern analysis for time domain boundary integral equations (TDBIE) traces its roots back to the seminal paper [6], where the bounds for boundary integral operators and layer potentials of the Helmholtz equation are derived by carrying out inversion using a Plancherel formula in anisotropic Sobolev spaces. In the regime of parabolic equations, the paper [29] converts the Laplace domain results to the time domain by bounding the inverse Laplace transform with pseudodifferential calculus, a theoretical tool not available for non-smooth domain problems. As already mentioned, [4] and [15] contain the seeds of a space-and-time coercivity analysis for the thermal boundary integral operators using anisotropic Sobolev spaces, while [29] adopts a separate strategy for the space and time variables, focusing on Sobolev regularity in space and Hölder regularity in time. Working on a different parabolic problem (Stokes flow around a moving obstacle), [5] gave an alternative proof of the bounds applicable to non-smooth boundaries, improving past results by revealing how all constants depend on time. Note that the usual approach was the study of mapping properties for the TDBIE and its inverse operator and for the associated potentials. Some kind of Laplace-domain coercivity was used to justify space or space-and-time Galerkin discretization. Instead, the work of [24] understood Galerkin semidiscretization in space as part of the problem set-up and its effects were analyzed as a continuous problem, instead of as a discretization of an existing continuous problem. The time-domain translation of this approach (using any variant of Laplace inversion or Payley-Wiener estimates) typically yields estimates that are suboptimal, due to the passage through the Laplace domain. (This effect can be seen for the TDBIE associated to the wave equation in [40].) What we do in this paper is related to the purely time-domain analysis of TDBIE initiated in [17] for wave propagation problems. This theory has seen different extensions and refinements: for instance, [37] extends the results to the Maxwell equations and [22] is the realization that a first order in time formulation makes the analysis much simpler. The goal of exploring purely time-domain techniques is multiple. First of all, in comparison with the Laplace domain approach, and when compared on the same type of bounds (the Laplace domain also provides estimates in weighted Sobolev norms that are not available with other techniques), time-domain estimates provide: (a) lower needs for the regularity of the input data to obtain the same estimates (i.e., refined mapping properties); (b) better bounds for the estimates as time grows. The time-domain analysis also emphasizes that a dynamical process is occurring through the entire discretization problem, a process that is hidden when we focus on transfer function estimates directly attached to the boundary integral operators.
The previous two paragraphs dealt with integral formulations, mapping properties, and semidiscretization in time. We now turn our attention to Convolution Quadrature. The multistep version of CQ, considered as a method to approximate causal convolutions and convolution equations, originated in [27]. A multistage version of the method was derived only a couple of years later in [28]. CQ techniques are now widely used in the realm of TDBIE, especially for wave propagation phenomena [34,8,25], but they are also useful in the context of TDBIE for parabolic problems [29,5]. The Laplace domain analysis of CQ has a black-box nature that makes it very attractive: it deals with general families of operators as long as their Laplace transforms (transfer functions) satisfy certain properties. However, as already observed in the seminal work of Lubich, the CQ process applied to TDBIE can be rewritten as the application of a background ODE solver to the associated PDE in the exterior domain. In fact, rewriting the CQ process as a time-stepping procedure expressed through ζ-transforms puts into evidence the fact that we are approximating a non-standard evolutionary PDE with non-homogeneous boundary conditions using an ODE solver. Along those lines, this paper offers two non-trivial contributions. First of all, we use functional calculus techniques and classical analysis of BDF methods [45] to show a direct-in-time analysis of BDF-CQ methods applied to the semidiscrete system of TDBIE. Second, we borrow heavily on difficult results by Alonso and Palencia [3] to offer an analysis of RK-CQ methods applied to the same problem. While we can prove estimates that improve the basic stage order (which is what you could obtain with direct Laplace domain analysis [28,9,10]), our numerical experiments will show that we are still slightly suboptimal and some additional work is needed.
The paper is organized as follows. Section 2 introduces the time domain boundary integral equation formulation (TDBIE) for the heat equation transmission problem. Section 3 proves the stability and convergence of the Galerkin-semidiscretization-in-space scheme. Sections 4 and 5 prove the convergence of BDFCQ and RKCQ in respectively. Finally, Section 6 provides several numerical experiments. An appendix presents some needed background material to ease readability.
Notation. For Banach spaces X and Y , B(X, Y ) will be used to denote the space of bounded linear operators from X to Y . We use standard function space notations: C k (I; X) for the space of k times continuously differentiable functions of a real variable in the interval I with values in the Banach space X, L 2 (O) for the space of square integrable functions on a domain O, and the Sobolev spaces If Γ is the boundary of a Lipschitz domain, H 1/2 Γ will be the trace space, H −1/2 Γ its dual, and ·, · Γ will denote the duality product of H Γ . We will use the following convention respectively. We will not have a special notation for the natural norm of H 1 ∆ (O). We will use the same notation (1.1) for the norms on Cartesian products of several copies of the same spaces. Finally, we will denote R + := [0, ∞).

Model problem and TDBIE Formulation
We are concerned with a transmission problem for the heat equation in free space in presence of a single homogeneous inclusion. Both the inclusion and the free space medium possess homogeneous and isotropic thermal transmission properties, characterized by two positive constants: κ as the thermal conductivity and ρ as the density scaled by heat capacity.
Let Ω − ⊂ R d (d = 2, 3) be a bounded Lipschitz domain with boundary Γ. The normal vector field ν : Γ → R d is defined almost everywhere on the boundary, pointing from the interior Ω − to the exterior domain Ω + := R d \ Ω − . We can thus define two trace operators where upper dots denote differentiation in time. We next give a crash course (based on [40]) on the few ingredients that are needed to have a rigorous setting for the weak definition of the heat boundary integral operators applied in the sense of distributions. Let F : C + := {s ∈ C : Re s > 0} → X be an analytic function such that F(s) X ≤ C F (Re s)|s| µ ∀s ∈ C + , where C F : (0, ∞) → (0, ∞) is non-increasing and is allowed to blow-up as a rational function at the origin, i.e., there exists a constant C > 0 and ≥ 0 such that C F (σ) ≤ Cσ − when σ → 0. It is then possible to prove [40,Chapter 3] that there exists an X-valued causal tempered distribution f whose Laplace transform is F, i.e., L{f } = F. The precise set of distributions whose Laplace transforms satisfy the above conditions is described in [40,Chapter 3] and denoted TD(X). These Laplace transforms (symbols, or transfer functions) include the ones that appear in parabolic problems (see [5]) where now F is well defined and analytic in C := C \ (−∞, 0] and satisfies where D F : (0, ∞) → (0, ∞) has the same properties as C F above. Since and therefore admits a unique solution. To see that, note that the associated bilinear form is coercive as This solution of (2.4) can be expressed using two s−dependent bounded operators acting on the data: Once again by definition, the following limit relations hold: ∂ ± ν S(s) = ∓ 1 2 I + K T (s), γ ± D(s) = ± 1 2 I + K(s).
Theorem 2.1. For s ∈ C + , denote σ := Re s 1/2 > 0 and σ := min{1, σ}. There exists a constant C only depending on the boundary Γ for all s ∈ C such that Proof. Since we can write the heat equation in Laplace domain as ∆u − (s 1/2 ) 2 u = 0 for s 1/2 ∈ C + , i.e., s ∈ C , s in the estimates in [24, Table 1] can be replaced by s 1/2 .
Applying the results of [40, Chapter 3], we can define the operator-valued distributions in the time domain through the inverse Laplace transform, using an inverse diffusivity parameter m > 0 When m = 1, the subscript will be omitted. As is well known, convolutions in time correspond to multiplications in the Laplace domain. For instance, if the Laplace transform of λ ∈ TD(H −1/2 Γ ) is Λ = L{λ}, then L(S m * λ) = S(s/m)Λ(s) and S m * λ ∈ TD(H 1 (R d )). The convolution operator λ → S m * λ is the heat single layer potential. More details about the distributional convolution can be found in [40,Section 3.2], with a more general theory given in [46]. The next theorem is the Green's representation theorem for the heat equation, which is a consequence of the analogous result in the Laplace domain.
Even though we will not need them for our numerical scheme, we now write down the explicit expression for the heat layer potentials. The time domain fundamental solution of the heat equation is The single layer potential is then given by while the integral form of the double layer potential is These integral operators are well defined for smooth enough densities λ and φ and they coincide with the distributional defintions. Precise mapping properties in anisotropic spacetime Sobolev spaces are given in the fundamental work of Martin Costabel [15]. The transmission problem in the sense of distributions is a weak version of (2.1). The data are now β 0 ∈ TD(H 1/2 Γ ) and β 1 ∈ TD(H −1/2 Γ ) and we look for u ∈ TD(H 1 The upper dot is now the distributional differentiation with respect to the time variable. The bracket in the right-hand sides of the equations in (2.6) clarifies where the equations hold. For instance, when we say ρu = κ∆u in L 2 (Ω − ), we mean that both sides of the equation are equal as L 2 (Ω − )-valued distributions. A rigorous understanding of such an equation requires the elementary but careful use of steady-state operators, like distributional differentiation in the space variables or the restriction of a function to a subdomain.
An integral system. We finally derive a system of time domain boundary integral equations (TDBIE) that is equivalent to the transmission problem (2.6). This follows exactly the same pattern as the work of Costabel and Stephan for steady-state (or time-harmonic) problems [16], recently extended to transmission problems for the wave equation [36]. Since the ideas are exactly the same as in those references, we will just sketch the process. We first choose the interior trace and normal derivative of u from (2.6) as unknowns We then define two scalar fields each of them defined on both sides of the boundary, and related to the solution of (2.6) by where the symbol χ O is used to the denote the characteristic function of the set O. In theory, this doubles the number of unknowns of the problem, even if we know that u − and u + vanish identically in Ω + and Ω − respectively. Later on, it will be clear that the doubling of unknowns is a natural byproduct of semidiscretization in space. The solution of (2.6) can be reconstructed as If we now substitute the representation formula (2.8) in (2.9c)-(2.9d), it follows that λ and φ satisfy We summarize the relations between the boundary integral equations and the partial differential equation in the following proposition. Its proof follows from elementary arguments using the jump relations of potentials and the definitions of the associate boundary integral operators.
is a solution of (2.10) and let is then a solution of (2.9) and u : is a solution of (2.10).

Galerkin semidiscretization in space
In this section we address the semidiscretization in space of the system of TDBIE (2.10), using a completely general Galerkin scheme, and the postprocessing of the boundary unknowns using the potential expressions (2.8).

The semidiscrete problem
We start by choosing a pair of finite dimensional subspaces (Note that following [24] we will only need X h and Y h to be closed.) Their respective polar sets are The semidiscrete method looks first for λ h ∈ TD(X h ), φ h ∈ TD(Y h ) satisfying a weakly-tested version of equations (2.10) The expression (3.1) is a compacted form of the Galerkin equations: when we write that the residual of the equations is in X • h × Y • h , we are equivalently requiring the residual to vanish when tested by elements of X h × Y h . Once the boundary unknowns have been computed, the potential representation yields two fields defined on both sides of Γ and satisfying the corresponding wave equations. If we subtract (3.1) by (2.10), we obtain the system satisfied by the error of unknown densities on the boundary The error corresponding to the posprocessed fields is easily derived by subtracting (2.8) from An exotic transmission problem. A transmission problem will encompass the solution of the semidiscrete system (3.1)-(3.2) and the associated error system (3.3)-(3.4). The problem looks for w − , w + ∈ TD(H 1 ∆ (R d \ Γ)) such that 2) is equivalent to the above system with λ = 0, φ = 0. On the other hand, if we set β 0 = 0, β 1 = 0, then the spatial semidiscrete error (e h − , e h + ) of (3.4) is the solution to (3.5). The main results for this section require some additional functional language and will be given in Section 3.3, after we have embedded a stronger version of the distributional system (3.5) in a framework of evolutionary problems on a Hilbert space.

Functional framework
The handling of the double transmission problem (3.5) (with two fields defined on both sides of the interface) can be carried out with theory of differential equations associated to the infinitesimal generator of an analytic semigroup. Consider first the following spaces To separate components of the elements of these spaces we will write w = (w − , w + ). Given a constant c = 0 (we will need c ∈ {ρ, ρ −1 , κ}) we will write cw := (cw − , w + ) for the associated multiplication operator acting on the first component of the vector. We will also consider the following bilinear forms and four copies of the boundary spaces, equipped with their product duality pairing The two-sided trace and normal derivative operators are remixed to transmission operators where the matrices Note that the first three components of γ D and the last three components of γ N appear in the transmission conditions of (3.5) and This integration by parts formula can be understood as a different way of rephrasing [36, formula (4.7)]. We finally consider the operator and the spaces which are respective polar spaces. In the coming paragraphs we will study the following problem: we are given data This is a strong form (restricted to the interval [0, ∞) and with strong derivatives, instead of distributional ones) of (3.5) when we choose χ D = (β 0 , 0, −φ, 0) and χ N = (0, −κλ, 0, β 1 ). The last ingredient for our framework consists of two spaces In some future arguments we will find it advantageous to collect the transmission conditions (3.9b)-(3.9c) in a single expression, using Proposition 3.1. With the above notation: admits a unique solution and The constant C depends only on Γ, ρ, and κ.
The operator A is maximal dissipative and self-adjoint.
(c) The operator A is the generator of a contractive analytic semigroup in H.
Proof. To prove (a), consider the coercive variational problem: The coercivity constant of the bilinear form in (3.12b) depends only on the constants κ and ρ if we use the standard (3.12) with smooth functions that are compactly supported in R d \ Γ, we can prove that ρw = κ∆w. Therefore, by (3.6), it follows that and therefore, by (3.6), If we define the identity (3.13) and a simple computation show that (3.14) In the sequel we will also use Hölder spaces C θ (R + ; X) for θ ∈ (0, 1), where X is a Hilbert space, and the seminorms |f | t,θ,X := sup The problem (3.9) admits a unique solution satisfyinġ Moreover, there exist constants independent of h such that for all t Proof. The proof is based on the decomposition of the solution of (3.9) into a sum w = w χ + w 0 , where w χ takes care of the data (using Proposition 3.1), while w 0 will be handled using a Cauchy problem (Theorem A.1). Let w χ (t) be the solution of (3.11) with g = 0, ξ D = χ D (t) and ξ N = χ N (t), and note thatẇ χ = wχ ∈ C θ (R + ; D), since we have applied a time-independent operator to the transmission data. Let now f := w χ −ẇ χ ∈ C θ (R + ; H), which satisfies f (0) = 0. Note that for all t ≥ 0 with a constant C depending exclusively on the parameters and geometry. Let finally w 0 : By Theorem A.1, w 0 and therefore w have the required regularity. Using the bound for w 0 in Theorem A.1 and (3.16a)-(3.16b), we can easily prove (3.15a). Using the bound for Aw 0 in Theorem A.1 and (3.16a)-(3.16c), we can prove that Proving (3.15b) from the above estimate is the result of a simple computation. Finally, in view of (3.14), the estimate follows and therefore (3.15c) is a simple consequence of (3.16a), (3.15a), and (3.15b).

Main results on the semidiscrete problem
The first step towards the analysis of the two problems that are hidden in (3.5) is the reconciliation of the solution of the classical differential equation (3.9) with a distributional form, where we look for w ∈ TD(D) such thaṫ for given data η ∈ TD(H Γ ). is equivalent to the variational problem for all s ∈ C + , which can be easily proved using the techniques of the proof of Proposition 3.1. We will prove that (3.19) is uniquely solvable and that its solution can be bounded as for some ν ≥ 0 and non-increasing C : (0, ∞) → (0, ∞) that is allowed to grow rationally at the origin. These statements imply that W = L{w} where w ∈ TD(D) (note that the needed bounds for the Laplacian of W (s) follow from equation (3.18)) and w satisfies (3.17), which is the inverse Laplace transform of (3.18).
In order to deal with (3.19) and (3.20), we proceed as follows. For fixed s ∈ C + , we consider the coercive transmission problem  [6] and an 'extension' to non-smooth boundaries can be found as Lemma 2.7.1 in [40]), it follows that We then consider the coercive variational problem What is left for the proof is very simple indeed. First of all, it is clear that W (s) := W D (s) + W 0 (s) is the solution to (3.19). Second, it is simple to see that which yields a bound for |||W 0 (s)||| |s| in terms of H N (s) −1/2,Γ and |||W D (s)||| |s| . Finally (3.22) can be used to prove (3.20).
The following process mimics the one in [22] and in [40,Chapter 7]. It involves two aspects: (a) an extension of zero of the data to negative values of the time variable; (b) a hypothesis on polynomial growth of the data. The reason to deal with (a) lies in the fact that the distributional equations (3.17) are for causal distributions of the real variable, not for distributions defined in the positive real axis. This is due to the fact that the heat potentials and operators have memory terms that involve the entire history of the process including values at time t = 0. The extension by zero to negative time will be done through the operator The reason why (b) is important is the fact that the equation (3.17) can be shown to have a unique solution in the space TD(D), which imposes some restrictions on the growth of the solution (and hence the data) at infinity.
Proof. The hypotheses imply that Eχ ∈ TD(H Γ ). The bounds of Proposition 3.2 imply that w is polynomially bounded as a D-valued function and therefore v := Ew ∈ TD(D). Finally, since w(0) = 0, it follows that Eẇ =v, which finishes the proof, since E commutes with any operator that does not affect the time variable.
We are almost ready to state and prove the two main results concerning the semidiscrete system: semidiscrete stability and an error estimate. To shorten up some of the expressions, to come, we introduce the bounded jump operator We also consider the function spaces tagged in the parameter θ ∈ (0, 1) Theorem 3.5. If β ∈ B 1+θ , ψ h = (φ h , λ h ) is the solution to the semidiscrete system of TDBIE (3.1) and u h = (u h − , u h + ) is given by the potential representation (3.2), then u h ∈ U θ , ψ h ∈ B θ , and we can bound is a collection of cummulative seminorms in B 1+θ .
Proof. If β = (β 0 , β 1 ), we define χ D := (β 0 , 0, 0, 0) and χ N = (0, 0, 0, β 1 ). Let u h be the solution to (3.17) with data η = (χ D , χ N ) and ψ h := J u h . We can identify u h and ψ h with the solution of (3.1) and (3.2). We also note that The bounds in the statement of the theorem follow from the fact that Proposition 3.4 identifies u h | R + with the solution of (3.9), with (χ D , χ N )| R + as data, and we can thus use the estimates of Proposition 3.2.
Proof. First of all, note that the errors ψ h − ψ can be defined as the solution of the semidiscrete TDBIE (3.3) and the associated potential errors e h := u h − u are given by a potential representation (3.4) using ψ h − ψ as input densities. If we define χ D = (0, 0, −φ, 0) and χ N = (0, −κλ, 0, 0), we can see that e h is the solution to (3.17) with data η = (χ D , χ N ) and that ψ h − ψ = J e h . Using the same arguments as in the proof of Theorem 3.5, we can prove that e h (t) 1,R d \Γ + J e h (t) H Γ ≤ C max{1, t}Acc(ψ, t, θ). Note now that we can decompose ψ h −ψ = (ψ h −Π h ψ)−(ψ −Π h ψ) and consider ψ −Π h ψ as the exact solution in the argument and ψ h − Π h ψ as its Galerkin approximation. In other words, if we input Π h ψ as exact solution of the TDBIE, then Π h ψ is also the solution of the semidiscrete TDBIE and the associated error is zero. Therefore, we can rewrite (3.26) as e h (t) 1,R d \Γ + J e h (t) H Γ ≤ C max{1, t}Acc(ψ − Π h ψ, t, θ), which proves (3.25).

Multistep CQ time discretization
In this section we introduce and analyze Convolution Quadrature schemes, based on BDF time-integrators, applied to the semidiscrete TDBIE (3.1) and to the potential postprocessing (3.2). All convolution operators -in the left and right hand sides of (3.1) and in the retarded potentials in (3.2) -will be treated with BDF-CQ.

The algorithm
Consider a constant k > 0 and a sequence of discrete time steps t n := nk for n ≥ 0 and let be the characteristic polynomial of the BDF(q) method. If we write the ζ-transform of a sequence of samples of a causal X-valued function v, Note that the approximation of the derivative is only good when v is a smooth causal function. In terms of the ζ-transform of data and of the fully discrete unknowns and then postprocess their output to build This short-hand exposition of the BDF-CQ method can be easily derived by taking the Laplace transform of (3.1) and (3.2) and substituting the Laplace transformed variable s by the discrete symbol k −1 δ(ζ). For readers who are not acquainted with CQ techniques, we explain in Appendix A.2 the meaning of formulas (4.3) and (4.4).

The analysis
We can think of (4.3)-(4.4) as a 'frequency-domain' system of BIE followed by potential postprocessing associated to a transmission problem with diffusion parameters ρκ −1 k −1 δ(ζ) and k −1 δ(ζ). Therefore, if we write u h,k n := (u h,k −,n , u h,k +,n ), χ(t) := ((β 0 (t), 0, 0, 0), (0, 0, 0, β 1 (t)), and recall (4.2), it follows that is the backward derivative associated to the BDF scheme. On the other hand, the semidiscrete solution u h satisfies very similar equations at the discrete timeṡ Therefore, the error e n := u h,k n − u h (t n ) satisfies the equations is the error associated to the finite difference approximation of the time derivative. Note that Theorem 4.1. If u h is smooth enough, then for all n ≥ 0 The proof of (4.10) is purely algebraic, based only on the fact that α 0 I − kA can be inverted. The rational functions in (4.11) are bounded at infinity and have all their poles at −α 0 < 0, which allows us to apply Theorem A.2 and show that However, if u h is smooth enough (as a function from R to H = L 2 (R d \ Γ) 2 ), Taylor expansion yields and that f n := ∂ k e n ∈ D h satisfies the recurrence ∂ k f n = Af n + ∂ k θ n . Using a simple argument on Taylor expansions and Theorem A.2, it follows that This can be used to give a bound for A e n = ∂ k e n − θ n and (3.14) then proves that Using this bound, (4.8), and the boundedness of J , (4.9b) follows.

Multistage CQ time discretization
In this section we introduce and analyze some Runge-Kutta based Convolution Quadrature (RKCQ) schemes for the full discretization of the semidiscrete system of TDBIE (3.1) and the potential postprocessing (3.2). Some background material and references on RKCQ and the needed Dunford calculus can be found in Appendices A.3 and A.4.

The algorithm and some observations
We consider an implicit s-stage RK method with Butcher tableau c Q b T c 1 a 11 a 12 · · · a 1s c 2 a 21 a 22 · · · a 2s . . .
c s a s1 a s2 · · · a ss b 1 b 2 · · · b s and stability function r(z) := 1 + zb T (I − zQ) −1 1, . . , 1) T ∈ R s and I is the s × s identity matrix. We will assume the following hypotheses on the RK method: (a) The method has (classical) order p and stage order q ≤ p − 1. We exclude methods where q = p for simplicity. (For instance, the one-stage backward Euler formula, which was covered as the BDF(1) method in the previous section, is not included in this exposition.) (b) The method is A-stable, i.e., the matrix I − zQ is invertible for Re z ≤ 0 and the stability function satisfies |r(z)| ≤ 1 for those values of z.
Assuming the usual simplifying hypothesis for RK schemes Q1 = c, stiff accuracy implies that c s = 1, that is, the last stage of the method is the step. Stiff accuracy also implies that (cf. [10, Lemma 2]) The matrix Q is invertible. This hypothesis and A-stability imply that the spectrum of Q is contained in C + .
Examples of the Runge-Kutta methods satisfying all the hypotheses above are provided by the family of s-stage Radau IIA methods with order p = 2s − 1 and stage order q = s. Given a function of the time variable, we will write . . .
to denote the s-vectors with the samples at the stages in the time interval [t n , t n+1 ]. We sample the boundary data and collect the vectors of time samples in formal ζ-series The unknowns for the fully discrete method can be collected in The pairs (λ h,k n , φ h,k n ) ∈ X s h × Y s h are computed using (4.3) where the symbol δ is now the matrix-valued RK differentiation operator (these two matrices can be easily seen to be equal using the stiff accuracy hypothesis) and the testing condition has to be modified, imposing that the residual is in ( The discrete potentials at the different stages can be computed using (4.4), with the new definition of δ(ζ). In all the expressions for the RK-CQ fully discrete equations, analytic functions are evaluated at k −1 δ(ζ) via Dunford calculus (see Appendix A.3). This is meaningful since the spectrum of δ(ζ) lies in C + for ζ small enough, which is due to the fact that δ(ζ) is a small (rank-one) perturbation of Q −1 and the spectrum of Q is in C + (see hypotheses (b) and (d) above).
Before we embark ourselves in the error analysis of the fully discrete method, which involves using quite non-trivial results from [3], we are going to make some important remarks that will be pertinent to the analysis.
(1) Using classical results on interpolation spaces on (bounded and unbounded) Lipschitz domains and the identification of Sobolev spaces with or without Dirichlet condition for low order indices (see [30,Theorem 3.33, Theorem B.9, Theorem 3.40]) it follows that for all µ < 1/2 and therefore (2) Because of the hypotheses on the RK method, the rational function R(z) := r(−z) satisfies the conditions of Theorem A.2 (it is bounded at infinity and has all its poles in the negative real part complex half-plane). We also know that the operator −A is self-adjoint and non-negative (Proposition 3.1(c)). Therefore, where we have used A-stability of the RK scheme. Consequently sup n r(kA) n H→H ≤ 1 ∀n, ∀k > 0.
The entries of the matrix-valued rational function z → (I+zQ) −1 are rational functions with poles in {z : Re z < 0} and converging to zero as |z| → ∞. Using Theorem A.2 in a similar way to how we proved (5.2) above, it follows that Here we have used traditional tensor product notation for Kronecker products of s × s matrices by operators.

Error estimates
The following proposition basically says that using RKCQ in the semidiscrete system of integral equations and then in the potential representation is equivalent to applying RK to the transmission problem associated to the heat equation satisfied by the potential fields.
Proposition 5.1. If U h,k n := (U h,k −,n , U h,k +,n ) ∈ D s is the vector of the internal stages of the RKCQ method in [t n , t n+1 ], then Proof. If we collect data and potential fields in it then follows that Looking at the time instances of the discrete-in-time equations enconded in the ζ-transformed equations (5.7) and (5.6b), the result follows.
Proposition 5.2 (L 2 error estimate for the steps). Let u h,k n = e T s U h,k n be the n-th step approximation provided by the RKCQ method. We have the bounds: if p = q + 1, then whereas if q ≤ p + 2 and ε ∈ (0, 1/4], Here, Proof. If P : H Γ → M ⊥ is the orthogonal projection onto M ⊥ (the orthogonality is with respect to the H Γ inner product), then the semidiscrete equations (see (3.9)) are equivalent tou (Note that the last condition is equivalent to Bu(t) − χ(t) ∈ M .) If we apply the RK method to equations (5.9), and we recall that the method is stiffly accurate, we obtain that the computation of the internal stages is given by These equations are clearly equivalent to equations (5.5), which have been shown to be equivalent to the equations satisfied by the fields obtained in the RKCQ method. In summary, we are dealing here with the direct application of the RK method to equations (5.9). This result is now a consequence of one of the main theorems of [3]. Unfortunately the reader will be now teleported from the middle of this proof to the core of a highly technical article. We will just give a translation guide to help with the application of the results of that paper. In our context, and taking advantage of our limitation to methods with p ≥ q + 1, the key theorem of [3] is Theorem 2. This result is stated for problems with homogeneous boundary conditions, but Section 4 of [3] explains how to handle the non-homogeneous boundary conditions and why the result still holds. The following translation table 1 4 − ε} 1 0 can be used to navigate [3,Theorem 2] and relate it to our particular problem. Note that when p ≥ q + 2, it is important to have D ⊂ V ⊂ [H, D(A)] 1/4−ε with bounded embeddings. This is a hypothesis needed in the application of the theorems and allows us to eliminate the estimates in terms of interpolated norms and write them in the natural Sobolev norm of V .
We note that the error for the case p = q + 1 can be written with the H = L 2 (R d \ Γ) 2 norm in the right-hand-side. We will keep the V = H 1 (R d \ Γ) 2 overestimate in this case for simplicity. Proposition 5.3 (L 2 error for the internal stages). With c(u h , t n ) defined as in Proposition 5.2, we have the bounds Proof. Let E n := U h,k n − u h (t n + c k) ∈ D s h = D(A) s and note that Using Taylor expansions and the fact that a method with stage order q satisfies Qc −1 = c for 1 ≤ ≤ q (powers of a vector are taken componentwise), we can easily prove that (5.12) Therefore, by (5.4), we can bound and the result follows by applying (5.12) and Proposition 5.2.
Proposition 5.4 (H 1 error and estimates for boundary unknowns). With the definition of c(u h , t n ) given in Proposition 5.2, we have the estimates Proof. From (5.11), we have and therefore by Proposition 5.3 and (5.12). Therefore, by (3.13) The result is therefore clear, since the errors in the boundary quantities can be derived from jump relations applied to E n and the H 1 (R d \ Γ) 2 norm can be bounded by the sum of the V seminorm and the H norm.

Numerical experiments
In order to confirm our theoretical findings, we conduct some numerical experiments in R 2 . For this we combine a frequency-domain Galerkin-BEM code with the fast CQ-algorithm from [11]. Note that a faster implementation can be achieved using the Fast and Oblivious CQ method [41]. For the discretization in space we use piecewise polynomial spaces; for X h we use discontinuous polynomials of degree p and for Y h we use globally continuous piecewise polynomials of degree p + 1. We denote these pairings as P p − P p+1 .
Testing on a manufactured solution. We start our experiments with the case where the exact solution can be computed analytically. This allows us to test the predicted convergence rates from Sections 4 and 5.
Since we are mainly interested in the performance of the CQ-schemes, i.e., the discretization in time we try to use a spatial discretization of higher order than the time discretization, while keeping the ratio k/h of timestep size and mesh-width constant. This means, whenever we cut the timestep size in half, we perform a uniform refinement of the spatial grid. See Table 1 for the degrees used. We note that for smooth solutions, we expect convergence rates in space of order O(k p+1 ) for the quantity λ − λ h L 2 (Γ) + φ − φ h H 1/2 (Γ) when using the space P p − P p+1 (see e.g. [39]).
The domain Ω − is the quadrialteral with vertices (0, 0), (1, 0), (0.8, 0.8), (0.2, 1). The thermal transmission constants are chosen to be m := ρ −1 κ = 0.8, κ := 1.2. We can then prescribe a solution by picking a source point outside the polygon x sc = (1.5, 1.6) and defining: For our computations, we solve the system up to the fixed end-time T = 4. When using a Runge-Kutta method, we expect a reduction of order phenomenon, which depends on the norm under consideration, see Proposition 5.4. Since they are easier to compute, we would like to consider the L 2 -errors for the Dirichlet and Neumann traces. For the Dirichlet-trace φ, the L 2 (Γ)-norm is weaker than the H 1/2 Γ -norm in Proposition 5.4. We therefore expect a slightly higher rate of convergence. Namely, if we use the multiplicative trace estimate (see e.g. [13, est. (1.6.2)]) we get: Considering the error quantity we therefore expect a rate of p e,φ := min q + 3/4, 3 4 p + q 4 (up to arbitrary ε > 0) for the convergence in time.
For the normal derivative, the L 2 norm is stronger than what is covered by our theory. Therefore, we also include an estimation for the true H where Π L 2 denotes the L 2 -projection onto the boundary element space X h . Since the exact solution is smooth and the space-discretization error is taken to be of higher order than the time discretization, this should give a good estimate for the H −1/2 -error. We predict rates of p e,λ := q. For multistep methods there is no reduction of order phenomenon, and we expect to see the full convergence order. We collect the expected rates in Table 1.
Compared to the experiments, our estimates in Proposition 5.4 do not appear to be sharp. It should be possible to improve the convergence estimate for A(u h (t n ) − u h,k (t n )) 0,R d \Γ to O(k q+1/4−ε ). Proving this estimate is more involved and part of future work. Using such sharper estimates, it should be possible to show an improved rate of p e,φ := min q + 1, 3 4 p + q 4 + 1 16 p e,λ := q + 1/4, for the L 2 -error of φ and the H −1/2 Γ -error of λ respectively. We include these rates as conjectured in Table 1.
In Figure 1 we observe that the Runge-Kutta based methods slightly outperform the expected rates p e,φ and p e,λ . Instead we get very good correspondence to the conjectured rates p e,φ and p e,λ . For multistep methods we see the full convergence rate, as was predicted in Section 4.
A simulation. We finally illustrate the use of our method for a simulation without known analytic solution. We choose the domain Ω − as a horseshoe shaped polygon, with a high conductivity parameter κ := 100. The density ρ was set to 1. We then placed point-sources of the form  on points x sc j uniformly distributed on a circle. Making the ansatz for the solution u tot = u src + u, we can compute u using our numerical method and recover u tot by postprocessing. See Figure 2a for the geometric setting and initial condition. (Note: in order to avoid the singularity at t = 0, we shifted the functions by a small time t lag := 0.001).
We then solved the evolution problem up to the final time T = 1. We used the BDF(4) method with a step-size of k := 1/2048. We used P 3 − P 4 elements with h ≈ 1/64. The results can be seen in Figure 2.

Conclusions
We have presented a collection of fully discrete methods for transmission problems associated to the heat equation in free space. The problem is reformulated as a system of time domain boundary integral equations associated to the heat kernel, thus reducing the computation work to the interface between the materials. The system is discretized with Galerkin BEM in the space variable and Convolution Quadrature (of multistep or multistage type) in time. Part of our work has consisted in dealing with the error analysis for the fully discrete method directly in the time domain, thus avoiding typically suboptimal estimates based on Laplace transforms. All the results have been presented for the case of a single inclusion, but the extension to multiple inclusions is straightforward. The case where the inclusion has piecewise constant material properties is more complicated and will be the aim of future work.

A.2 BDF-CQ
Let δ be the backward differentiation symbol introduced in (4.1). When F : C → B(X 1 , X 2 ) is an operator-valued analytic function, is given by the Taylor expansion of the left hand side about ζ = 0 since δ(ζ) takes values in the domain of F for small |ζ|. (This is due to the A(θ)-stability of the BDF formulas of order less than or equal to six.) We thus obtain a sequence of operator-valued coefficients ω F n (k) ∈ B(X 1 , X 2 ).
Then given G(ζ) := ∞ n=0 g n ζ n : C → X 1 , when we multiply F( 1 k δ(ζ))G(ζ) = ∞ n=0 n m=0 ω F n−m (k)g m ζ n , (A. 4) we are just computing the discrete causal convolution of the sequence of operators {ω F n (k)} to the discrete sequence {g n }. Computational strategies for efficient implementation of multistep CQ (in the context we find it in (4.3)-(4.4)) can be found in the literature. The lecture notes [21] contain a simple introduction to the topic.

A.3 Rudiments of functional calculus
We here introduce some minimun requirements on (Dunford-Riesz) functional calculus needed for our work. First of all, here is a result concerning bounds for rational functions of nonnegative self-adjoint operators. The theorem is a consequence of more general results related to functions of operators, which can be found in general introductions to functional calculus (see, for instance, [20, Corollary 7.1.6, Theorem 2.2.3]).
Theorem A.2. Let B : D(B) ⊂ X → X be a self-adjoint non-negative operator in a Hilbert space X. Let P, Q ∈ P(C) be two polynomials such that the rational function R := P/Q is bounded at infinity (deg P ≤ deg Q) and has all its poles in {s ∈ C : Re s < 0}. Then R(B) : X → X is bounded and We will also need the evaluation of analytic functions on matrices, using Dunford calculus. Let F : O ⊂ C → X be an analytic function defined on a simply connected open set of C. Let M be a matrix whose spectrum is contained in O. We then define

A.4 RK-CQ
Let F : C + → B(X 1 , X 2 ) be an operator-valued analytic function and consider the expansion (A.3), where now δ(ζ) ∈ R s×s is defined by (5.1). Since for ζ small δ(ζ) is a small perturbation of Q −1 and Q has its spectrum contained in C + , then F(k −1 δ(ζ)) can be defined using functional calculus as in Section A.3. The expansion is then a simple Taylor expansion about the origin for an analytic function with values in B(X 1 , X 2 ) s×s . The coefficients of this expansion can be given by the Cauchy integrals where C is any simple positively oriented closed contour in C + surrounding the spectrum of Q −1 . Discrete causal convolutions in the form (A.4) can be made now for sequences {g n } in X s 1 . A simple practical introduction to Dunford calculus related to RK-CQ methods can be found in [21,Section 6.3]. Practical computational strategies involve diagonalizing δ(ζ), cf. [7,12,21].