A wave kernel that cancels at a curved boundary

An image point gives an exact Dirichlet kernel for a flat wall. For a curved wall, the useful replacement is a second coordinate chart: its radial curves travel to the wall and then reflect. The incident and reflected terms have the same singular distribution on the wall. Their smooth amplitudes can therefore be matched there, before any remainder is estimated.

The scale matters. A source at distance $\varepsilon$ from the wall is examined for times and path lengths smaller than $4\varepsilon$. Dividing space and time by $\varepsilon$ turns the source into a fixed point and the wall into a smooth perturbation of a plane. This makes the construction uniform as the source approaches the boundary. It does not include tangent reflection or arbitrarily many reflections.

1. The equation and the two pieces of input

Let $X$ be a smooth manifold with smooth boundary, and work in a coordinate neighborhood of a compact subset $K\subset\partial X$. An open bounded smooth domain is a particular case. The spatial operator has positive scalar principal symbol and, in a local frame of a complex bundle of rank $r$, is $$P=-\partial_j\bigl(g^{jk}(q)\partial_k\bigr)I_r +b^j(q)\partial_j+c(q). \tag{R1}$$ Repeated spatial indices are summed. The matrix $(g^{jk})$ is real, symmetric, smooth and positive definite up to the boundary; $b^j,c$ are arbitrary smooth complex matrices, acting on the left. The metric on vectors is $g=(g_{jk})=(g^{jk})^{-1}$. Coefficients are extended smoothly across the boundary, retaining positive definiteness after reducing the neighborhood. Every conclusion is local over $K$; no global uniform ellipticity or global trivialization is assumed.

We use coordinate distributions and the source normalization $$\rho(q)=(\det g^{jk}(q))^{-1/2},\qquad \delta_\mu(q,y)=\rho(y)^{-1}\delta_y(q),\qquad d\mu=\rho\,dq. \tag{R2}$$ Thus a kernel acting with $d\mu(y)$ has the identity source $\delta_\mu I_r$. In a new coordinate chart the density and the coordinate Dirac mass must both be transformed. The wave operator is $\partial_t^2+P$, with negative spatial Laplacian in the flat case.

The construction uses the following basic facts and earlier lessons.

• Calculus and distributions. smooth finite-dimensional calculus, the inverse and implicit function theorems, Taylor’s formula with remainder, smooth cutoffs, compactness, and local smooth dependence for the geodesic ordinary differential equation. Distributional differentiation, smooth multiplication, coordinate change and Fourier transformation have their usual meanings. The required linear matrix transport is proved below.
• Interior geometry and transport. Sections 4–6 of Building a local inverse from radial singularities provide geodesic coordinates and jointly smooth matrix amplitudes $U_\nu(q,y)$, with $U_0(y,y)=I_r$, obtained from the negative-principal operator (R1). Only the local geometry and these transport amplitudes are used. No boundary result or reflection coefficient is imported.
• Causal wave kernels. Causal kernels, initial data, and short-time geometry supplies the distributions $E_\nu$, $\nu\in\mathbb N_0$, on $\mathbb R_t\times\mathbb R^n_z$, with the normalization $$\begin{split} (\partial_t^2-\Delta_z)E_0&=\delta_{(0,0)},\\ (\partial_t^2-\Delta_z)E_\nu&=\nu E_{\nu-1},\qquad -2\partial_{z_j}E_\nu=z_jE_{\nu-1}\quad(\nu>0),\\ E_\nu(at,az)&=a^{2\nu+1-n}E_\nu(t,z),\qquad a>0. \end{split} \tag{R3}$$ They are rotationally invariant, supported in $t\geq|z|$, and smooth in $t\geq0$ with values in spatial distributions. Their one-sided initial derivatives vanish through order $2\nu$, and the next is $\nu!\delta_0$. Away from the cone vertex, the only singularity is a conormal distribution to $t=|z|$. For $\nu>(n-1)/2$, their causal cone formula is a constant multiple of $\boldsymbol1_{t>0}(t^2-|z|^2)_+^{\nu-(n-1)/2}$; the wave unit proves integer regularity strictly below that exponent, including at the vertex. The Hölder refinement is proved in Section 7 below. The same prerequisite supplies the finite incident identity obtained from (R3) and the preceding interior amplitudes.

The reflected proof applies to every finite rank and retains the order of matrix multiplication. It is a Dirichlet theorem for a scalar principal symbol; it does not close general complementing boundary systems or operators with several distinct principal propagation cones.

2. A uniform coordinate neighborhood at boundary distance one

Choose local boundary coordinates $p$ and follow the inward unit normal geodesic from each boundary point. This gives boundary-normal coordinates $(p,a)$, with $$X=\{a\geq0\},\qquad g_{jn}=g_{nj}=0\ (j<n),\qquad g_{nn}=1. \tag{R4}$$ Here and below an index smaller than $n$ is tangential; the normal coordinate itself is written $q_n$ when ambiguity is possible. To justify the construction, at distance zero the differential of the normal-geodesic map consists of a tangential basis and the inward normal, so it is invertible. Smooth geodesic dependence and the inverse function theorem give a collar. Orthogonality in (R4) follows by differentiating the inner product of the normal velocity and a tangential variation field: metric compatibility and the geodesic equation make this derivative zero, and its initial value is zero. The normal velocity has constant unit length. Finitely many such neighborhoods cover $K$, and reducing their thickness gives uniform bounds. In this collar the coordinate $q_n$ is the distance to the boundary: every curve reaching the boundary has length at least the change of $q_n$, since (R4) makes its speed at least $|\dot q_n|$, and a normal segment achieves equality.

Write the center as $y=(p,\varepsilon)$, $\varepsilon>0$, and set $$q=(p,0)+\varepsilon x,\qquad t_{\rm phys}=\varepsilon t, \qquad e=(0,\ldots,0,1), \qquad G_{p,\varepsilon}(x)=(g_{jk}((p,0)+\varepsilon x)). \tag{R5}$$ The capital $G$ here denotes the covariant metric of the scaled space. On every fixed bounded set of $x$, this metric is smooth jointly in $(x,p,\varepsilon)$ down to $\varepsilon=0$. At zero it is constant, has normal coefficient one and no tangential-normal cross terms. In particular it need not be tangentially Euclidean in the original coordinates.

Let $$A_y(S)=\exp^{G_{p,\varepsilon}}_e S,\qquad B_y=\{S:|S|_{G_y(e)}<4\}. \tag{R6}$$ For sufficiently small $\varepsilon$, the map is defined and is a diffeomorphism on a neighborhood of the closure of $B_y$, uniformly over compact subsets of a boundary chart. Here is a useful quantitative version of that assertion. The Christoffel symbols of $G_y$ are $O(\varepsilon)$, with all fixed derivative bounds uniform; on the slightly larger ball $|S|_{G_y(e)}\leq 4+\eta$, the geodesic equation and its differentiated integral equations give $$A_y(S)=e+S+O(\varepsilon),\qquad d_SA_y=I+O(\varepsilon), \tag{R7}$$ uniformly with each fixed number of derivatives. Choose a common Euclidean velocity ball containing these ellipsoids. On it, reduce $\varepsilon$ until the derivative error has norm below $1/2$. Integration along a line segment proves $|A_y(S)-A_y(T)|\geq|S-T|/2$. This proves injectivity as well as local invertibility; it is stronger than applying the inverse function theorem separately at each point. Compactness supplies the same choices for the finitely many charts. Ordinary exponential coordinates around $y$ follow by undoing (R5). In particular $(y,s)\mapsto(y,\exp_y s)$ is a local diffeomorphism near the zero section and the diagonal.

3. Which rays hit, and how they return

Let $F_y(S)$ be the last coordinate of $A_y(S)$. By (R7), $$F_y(S)=1+S_n+O(\varepsilon),\qquad d_SF_y=dS_n+O(\varepsilon). \tag{R8}$$ The level surface $\Sigma_y=\{F_y=0\}\cap B_y$ divides the velocity ball into incident endpoints in $X$ and endpoints beyond the wall. Denote the latter region, including the hitting surface, by $B_y^-$.

Every $S\in B_y^-$ has one first hitting time $\lambda_y(S)\in(0,1]$, defined by $$F_y(\lambda_y(S)S)=0. \tag{R9}$$ The word “first” can be verified uniformly. A zero with $|S|_{G_y(e)}<4+\eta$ has $\lambda S_n=-1+O(\varepsilon)$. Consequently $S_n\leq-c<0$, $\lambda\geq c>0$, and $\partial_\lambda F_y(\lambda S)=S_n+O(\varepsilon)\leq-c$, for a uniform $c$. The entire relevant segment has negative normal velocity, since its velocity differs from $S$ by $O(\varepsilon)$. Hence it crosses at most once. The implicit function theorem proves smoothness of $\lambda_y$, including up to $\Sigma_y$. At $\varepsilon=0$, $\lambda_y=-1/S_n$. The bounds also show that an endpoint with $F_y>0$ has not secretly crossed and returned before time one: if its normal velocity is ever capable of reaching the wall in the allowed time, that velocity has the same strict negative sign along the segment.

Let $z=A_y(\lambda S)$ and let $T$ be the geodesic velocity at the hit when the geodesic is parameterized on the same unit-time scale as $S$. Reflect it in the metric normal: $$T^r=T-2\langle T,N_z\rangle_{G_y}N_z,\qquad R_y(S)=\exp_z^{G_y}\big((1-\lambda)T^r\big). \tag{R10}$$ In (R4) this changes only the sign of the last tangent component. The outgoing velocity has uniformly positive normal component, and its variation on the remaining bounded interval is $O(\varepsilon)$; it therefore has no second hit. Both the hitting and reflected constructions are smooth up to $\Sigma_y$, and $$R_{p,0}(S)=(S',-1-S_n). \tag{R11}$$ The same derivative and line-segment argument as in (R7), applied to a slightly larger region $S_n<-1+c_1$, proves that $R_y$ is a diffeomorphism onto its image $\Omega_y$. The extension to this region is obtained by allowing the implicit hitting parameter a little beyond one; the actual reflected segment uses only $\lambda\leq1$. Thus a collar of the reflecting surface causes no differentiability problem. The image contains a neighborhood of $e$: the vector $-2e$ gives the normal ray down to the wall and back to the center, with total length two, exactly in boundary-normal coordinates. Define the reflected length by $$r_y^r(x)=|R_y^{-1}(x)|_{G_y(e)},\qquad x\in\Omega_y. \tag{R12}$$ It is smooth and bounded away from zero there, locally uniformly with parameters. The incident distance $r_y(x)=|A_y^{-1}(x)|_{G_y(e)}$ is smooth off $x=e$. At the wall the two lengths agree, since the reflected segment has length zero. There is no assertion of a single smooth chart obtained by joining $A_y$ and $R_y$ across the switch. Each branch is used on its own domain.

4. The radial identity survives a reflection

For a smooth family of geodesic segments, the variation of length equals the final unit momentum paired with final displacement minus the initial unit momentum paired with initial displacement. This follows by differentiating the length integral, using metric compatibility, and integrating the velocity-variation term by parts; the interior term is the covariant acceleration and vanishes for a geodesic. Apply the formula to both segments of the broken path. The two terms at the moving hit point $z$ add to $$\langle T/|T|-T^r/|T^r|,\dot z\rangle_{G_y}=0, \tag{R13}$$ because the vector difference is normal and $\dot z$ is tangent to the boundary. Thus the derivative of the total length is just its final unit momentum paired with the final endpoint variation. This proves orthogonality of each reflected wavefront to its final rays; it also proves the full identity needed for transport, not just orthogonality at one endpoint.

Indeed put $H_0=G_y(e)$ and let $\widehat G_y=R_y^*G_y$. Along a radial change of $S$, the endpoint moves along the outgoing geodesic with velocity corresponding to the entire path length. Since the squared total length is $S^tH_0S$, (R13) gives $$\widehat G_y(S)(S,W)=H_0(S,W)\quad\hbox{for every }W, \qquad \widehat G_y(S)S=H_0S. \tag{R14}$$ Writing $a_y=\widehat G_y^{-1}$, this is $a_y(S)H_0S=S$. The identical calculation for a segment without reflection proves the incident identity. An orthonormalizing linear change is optional; retaining $H_0$ makes the density normalization transparent.

Pull the actual operator $P_y$ into the reflected chart and write it uniquely in divergence form relative to $dS$: $$Q_y=-\partial_{S_j}(a_y^{jk}\partial_{S_k})I_r +\beta_y^j\partial_{S_j}+\gamma_y, \qquad h_y(S)=\beta_y^j(S)(H_0S)_j. \tag{R15}$$ The pullback of a differential operator acting on functions includes its coordinate first-order terms. They are part of $\beta_y$, even when the original operator is a Laplace operator. In particular $$Q_{p,0}=-g^{jk}(p,0)\partial_{S_j}\partial_{S_k}I_r, \qquad h_{p,0}=0,\qquad\gamma_{p,0}=0. \tag{R16}$$ The minus sign in (R16) is fixed by (R1) and (R3).

For a radial smooth function $f(S^tH_0S)$, (R14) gives $a_y\nabla f=2Sf'$. Taking the divergence shows that the principal divergence operator on such a function is its constant-metric counterpart. The same assertion holds for the cone distributions $E_\nu(t,|S|_{H_0})$. On this reflected domain $|S|$ never vanishes, so one can use $(t,|S|,\omega)$ coordinates: multiplication and differentiation of the one-dimensional distribution in its radial argument obey the ordinary chain rule. Equivalently, regularize that one-dimensional distribution, apply the smooth identity and pass to distributions. No unverified multiplication at the cone is required.

For $\nu>0$, the product rule and (R3) now give $$(\partial_t^2+Q_y)(vE_\nu) =(Q_yv)E_\nu+ \big(\nu v+S\cdot\nabla_Sv-\tfrac12h_yv\big)E_{\nu-1}. \tag{R17}$$ The coefficient acts on the left of the matrix amplitude. For $\nu=0$, the same product rule has no point source on this domain, and its radial derivative term vanishes precisely when $2S\cdot\nabla v=h_yv$. These calculations explain which operator sign and which multiplication order belong in the reflected transport equations.

5. Transport starts at the wall

Let smooth matrix data $d_\nu$ be prescribed on $\Sigma_y$. Set $v_{-1}=0$. There is a unique smooth family on $B_y^-$, including its boundary, satisfying $$(S\cdot\nabla_S+\nu-\tfrac12h_y)v_\nu=-Q_yv_{\nu-1}, \qquad v_\nu|_{\Sigma_y}=d_\nu. \tag{R18}$$ The assertion is local on the bounded velocity region under discussion, and smoothness includes every parameter in a compact set.

Here is an explicit solution and uniqueness proof. Fix $S$, put $\lambda=\lambda_y(S)$, and let $M(r,s)$ solve $$\partial_rM(r,s)=\frac{h_y(rS)}{2r}M(r,s),\qquad M(s,s)=I_r, \qquad \lambda\leq r,s\leq1. \tag{R19}$$ The coefficient is smooth, since $\lambda\geq c>0$. Iterating its integral equation gives an ordered integral series whose term of order $k$ has norm at most $C^k/k!$. Each fixed number of derivatives in $(S,y,r,s)$ contributes at most a polynomial factor in $k$, so the differentiated series converges uniformly as well. The reverse equation $\partial_rM(s,r)=-M(s,r)h_y(rS)/(2r)$ proves invertibility. This argument does not commute any matrices.

If $f_\nu=-Q_yv_{\nu-1}$, variation of constants gives $$v_\nu(S)=\lambda^\nu M(1,\lambda)d_\nu(\lambda S) +\int_\lambda^1r^{\nu-1}M(1,r)f_\nu(rS)\,dr. \tag{R20}$$ To verify the formula, on a ray multiply the equation by $r^\nu M(\lambda,r)$; its derivative is $r^{\nu-1}M(\lambda,r)f_\nu(rS)$. Integration gives (R20). The boundary data and the equation determine the solution uniquely on that ray, and every point lies on exactly one such outgoing ray. Smoothness of the endpoint $\lambda(S,y)$, the integrand and the ordered series proves smoothness of the solution. Induction in $\nu$ proves the entire statement. The factor $\lambda^\nu$ is necessary even in the scalar case. Unlike transport from a center, transport here has a nonzero initial radius and allows an independent datum for every $\nu$.

6. Scaling the incident wave and matching all coefficients

Under (R5), the scaled spatial operator is $$P_y=-\partial_{x_j}(g^{jk}((p,0)+\varepsilon x)\partial_{x_k})I_r +\varepsilon b^j((p,0)+\varepsilon x)\partial_{x_j} +\varepsilon^2c((p,0)+\varepsilon x). \tag{R21}$$ In fact $\partial_{q_j}=\varepsilon^{-1}\partial_{x_j}$ and $\partial_{t_{\rm phys}}=\varepsilon^{-1}\partial_t$; multiplying the original equation by $\varepsilon^2$ proves every factor in (R21). Let $$A_\nu(x,y)=\varepsilon^{2\nu}U_\nu((p,0)+\varepsilon x,y), \qquad I_N(t,x,y)=\sum_{\nu=0}^N A_\nu(x,y)E_\nu(t,r_y(x)). \tag{R22}$$ Then the finite incident identity is $$(\partial_t^2+P_y)I_N =\rho(y)^{-1}\delta_{(0,e)}I_r +(P_yA_N)E_N(t,r_y(x)). \tag{R23}$$ For a direct scaling check, multiply the physical kernel by $\varepsilon^{n-1}$ before substituting (R5). Its $\nu$-th term acquires the factor $\varepsilon^{n-1}\varepsilon^{2\nu+1-n}=\varepsilon^{2\nu}$. The operator contributes $\varepsilon^2$, and the spacetime delta contributes $\varepsilon^{-n-1}$. Their product is one. Thus the source is at $(0,e)$, not at the unscaled coordinate value $(0,y)$. In (R23), $P_yA_N=\varepsilon^{2\nu}P_y[U_N((p,0)+\varepsilon\,\cdot,y)]$ with $\nu=N$; the operator differentiates the composed function of $x$.

Use (R18) with data $$d_\nu(S,y)=A_\nu(A_y(S),y),\quad S\in\Sigma_y, \qquad V_\nu(x,y)=v_\nu(R_y^{-1}(x),y). \tag{R24}$$ All these functions are smooth down to $\varepsilon=0$. At zero, the incident data are $d_0=I_r$ and $d_\nu=0$ for $\nu>0$, because $U_0(y,y)=I_r$ and the positive-index factor is $\varepsilon^{2\nu}$. Equations (R16), (R18) and uniqueness give $$V_0(x,p,0)=I_r,\qquad V_\nu(x,p,0)=0\quad(\nu>0). \tag{R25}$$ The latter conclusion uses $Q_{p,0}I_r=0$ at the first induction step. Smooth extension of $V_\nu$ through $\varepsilon=0$ does not mean it contains the factor $\varepsilon^{2\nu}$; the differentiated transport equations will show a curvature obstruction to that stronger statement.

7. A finite Dirichlet parametrix, with its exact error

On the common incident and reflected coordinate region put $$K_N(t,x,y)=\sum_{\nu=0}^N \{A_\nu(x,y)E_\nu(t,r_y(x)) -V_\nu(x,y)E_\nu(t,r_y^r(x))\}. \tag{R26}$$ There is a sufficiently small collar over $K$, chosen independently of $N$, such that this is defined for $r_y^r(x)<4$ and $t<4$. The coefficients for every fixed finite $N$ are smooth jointly in $(x,p,\varepsilon)$ down to zero. The distribution has zero Dirichlet trace and satisfies $$\begin{split} (\partial_t^2+P_y)K_N &=\rho(y)^{-1}\delta_{(0,e)}I_r+\mathcal R_N,\\ \mathcal R_N &=(P_yA_N)E_N(t,r_y(x)) -(P_yV_N)E_N(t,r_y^r(x)). \end{split} \tag{R27}$$ For every nonnegative regularity index satisfying $$\mu<N-\frac{n-1}{2}, \tag{R28}$$ the error is locally $C^\mu$, uniformly on compact subsets of the scaled variables and the closed parameter collar. For integral $\mu$, this means continuously differentiable through order $\mu$; for nonintegral $\mu$, it means the corresponding local Hölder class. If no nonnegative $\mu$ satisfies (R28), the identity still holds as a distributional finite parametrix, without a positive regularity claim for its error.

Proof of cancellation. At $x_n=0$, equations (R12) and (R24) give equal lengths and equal coefficients term by term. The traces exist as distributions in $(t,x',y)$. One can check this without a general restriction theorem. The lengths at the boundary are uniformly positive. For such a smooth positive length $r$, (R3) is locally a one-dimensional distribution in $t-r$, smoothly depending on all other variables: in the expression involving $t^2-r^2$, the derivative with respect to $t$ is $2t\ne0$ near its forward cone. Pullback is consequently defined by a change of the time variable. Outside that cone the distribution is smooth. Restricting the smooth length and amplitude to $x_n=0$ gives the claimed trace, and the two traces coincide. The apparent cone singularity therefore does not obstruct the cancellation.

The incident chart is available on this whole reflected region. Work first with the slightly larger normal ball used in (R7). The incident Gauss identity bounds the length of any curve escaping that ball below by its radius: in radial coordinates its speed is at least the absolute radial speed. Thus a path of length less than four cannot escape the larger ball. Inside it the same estimate proves that the incident radial geodesic minimizes length. Applying this to the once-broken path gives $r_y(x)\leq r_y^r(x)<4$. Hence no additional, unmentioned intersection restriction reduces the stated reflected domain.

Apply (R17) to the reflected sum. Equation (R18) cancels each consecutive pair; the term of index zero has no delta source because the reflected chart excludes $S=0$. The only reflected term left is $(P_yV_N)E_N$. Subtract it from (R23) to obtain (R27). The same reasoning covers $N=0$. Causality follows from the supports in (R3). Near $t=0$ the reflected sum vanishes, since its length is bounded below, so the initial point singularity is exactly the incident one. In the Euclidean case with zero lower-order coefficients, (R25) holds for every $\varepsilon$, every higher coefficient vanishes, and (R26) becomes the exact image kernel $$E_0(t,x-e)-E_0(t,x-e^*),\qquad e^*=(0,\ldots,0,-1). \tag{R29}$$ Its reflected source lies outside the interior half-space. This also directly verifies the flat equation, causal support, initial source and boundary trace.

Proof of the error estimate. All coefficients in (R27), including their parameter derivatives, are bounded on each of the stated compact sets. Near the incident center, the map from $(x,y)$ to its incident velocity and $y$ is a smooth diffeomorphism. A smoothly parameterized positive square root of $H_0$ reduces the constant metric there to the Euclidean metric; such a square root is smooth, for instance by the implicit equation $L^2=H_0$ on positive matrices, whose derivative $Z\mapsto LZ+ZL$ is invertible after diagonalizing $L$. The regularity input for $E_N$, composed with this diffeomorphism and multiplied by smooth coefficients, gives (R28) at the incident cone vertex as well as away from it. The reflected cone has no vertex in the chart and uses the same argument or the one-dimensional representation just used for its trace. No factor $\varepsilon^{-1}$ is introduced: these are the scaled variables in which the entire family is smooth.

Here is the additional argument for a nonintegral $\mu=k+\alpha$, with $k\geq0$ integral and $0<\alpha<1$. Set $a=N-(n-1)/2>k+\alpha$. On a fixed annulus about the cone vertex, the $k$-th derivatives of the causal cone function are $\alpha$-Hölder: near its nonzero cone this follows by differentiating $q_+^a$ and using $|u_+^\beta-v_+^\beta|\leq C|u-v|^\alpha$ on bounded intervals when $\beta>\alpha$; away from the cone the function is smooth. This elementary inequality follows from concavity when $0<\beta\leq1$ and the mean value estimate when $\beta\geq1$. By homogeneity, on an annulus of radius $r\leq1$ the Hölder bound for a derivative of order $k$ scales by $r^{2a-k-\alpha}$. This exponent is positive. For two points whose separation is less than half their larger radius, use that rescaled annular estimate. For a pair separated by at least half that radius, the bound $|D^kE_N|\leq Cr^{2a-k}$, together with the zero extension at the vertex, gives the same $C|h|^\alpha$ estimate. These two cases prove local $C^{k,\alpha}$ regularity across the vertex as well. Smooth coordinate changes, tensoring with a constant in the parameters, and smooth coefficient multiplication preserve this class, completing (R28) for every specified $\mu$.

The strict bound in (R28) follows from the behavior in the explicit causal radial distribution away from the vertex, which is a constant multiple of $(t^2-r^2)_+^{N-(n-1)/2}$. Differentiation at the cone loses one power for each derivative; at an integer power the critical derivative can jump. A smooth coefficient does not generally remove that jump. For a noninteger exponent the corresponding Hölder endpoint is attained; for a positive integer exponent the preceding derivative is Lipschitz, while the next derivative can still jump. The complete proof through the vertex is in Section 11 of the linked wave lesson. Increasing $N$ yields any prescribed finite regularity. It does not turn a given finite sum into an exact inverse. $\square$

Returning to physical variables gives the kernel $\varepsilon^{1-n}K_N(t_{\rm phys}/\varepsilon,(q-(p,0))/\varepsilon,y)$ in the corresponding neighborhood. Its error estimates have the powers obtained by this formula and the chain rule; the uniform statement (R28) is in the scaled coordinates. Coefficient extension supplies the incident construction with all centers up to the boundary, while the incident term alone can solve the local mixed problem only before the wave reaches the wall. The reflected construction deals with times less than a fixed multiple of the boundary distance. A construction for times small relative to its square root would need additional control of approaching tangency and is not proved here.

Editorial strengthening of the finite error. Retain (R26)–(R28) and all their hypotheses. Set $a=N-(n-1)/2>0$. For nonintegral $a=k+\alpha$, with integer $k\geq0$ and $0<\alpha<1$, the exact error belongs jointly to $C^{k,\alpha}$ on the stated compact scaled sets and closed parameter collar. For positive integer $a=\ell$, it belongs to $C^{\ell-1,1}$, meaning Lipschitz derivatives at the preceding order, not $C^\ell$. The full cone proof, with its gamma constant, all derivative terms, vertex argument, sharpness and joint parameter transfer, is in Section 11 of Causal kernels, initial data, and short-time geometry. In (R27) the incident term uses the smooth normal chart and the original positive square root; the reflected length stays away from zero. Smooth multiplication preserves the class in the original matrix order. This gives the claimed endpoint for each original error term.

Reflected-error proof 5. The original finite reflected error and its extra distance factor.

Keep the original scaled operator $P_y$ of R21, incident amplitude $A_N$=$\varepsilon ^{2N}$$U_N$((p,0)+εx,y), $y=(p,\varepsilon )$, and all ordered reflected transport coefficients $V_N$. The exact error remains

$$\mathcal R_N=(P_yA_N)E_N(t,r_y(x)) -(P_yV_N)E_N(t,r_y^r(x)). \tag{ER1}$$

Acting on a composed smooth function, the complete operator scaling in R21 gives the following identity,

$$P_y[U_N((p,0)+\varepsilon x,y)] =\varepsilon^2(P_qU_N)((p,0)+\varepsilon x,y),$$

, since y is held fixed during the q differentiation. Consequently

$$\begin{split} P_yA_N&=\varepsilon^{2N+2}(P_qU_N)((p,0)+\varepsilon x,y),\\ P_yV_N(x,p,\varepsilon)&=\varepsilon B_N(x,p,\varepsilon),\\ B_N(x,p,\varepsilon)&=\int_0^1 \left.\partial_s\big[P_{p,s}V_N(x,p,s)\big] \right|_{s=\theta\varepsilon}\,d\theta. \end{split} \tag{ER2}$$

For the second equality, R25 gives $V_0(x,p,0)=I_r$ and $V_N(x,p,0)=0$ for $N>0$, both spatially constant, and R21 at ε=0 is the constant principal operator with no drift or potential. Thus $P_{p,0}V_N(x,p,0)=0$ for every $N\geq 0$. Taylor’s integral formula proves ER2 exactly. All spatial and p derivatives of $B_N$ are bounded on the original compact sets; every fixed ε derivative is also bounded because its integrand is smooth. No claim that $V_N$ has incident vanishing order $\varepsilon ^{2N}$ is made.

Keeping both original terms visible, ER1 therefore also has the exact expression

$$\mathcal R_N=\varepsilon^{2N+2}(P_qU_N)((p,0)+\varepsilon x,y) E_N(t,r_y(x)) -\varepsilon B_N(x,p,\varepsilon)E_N(t,r_y^r(x)). \tag{ER3}$$

For $a>0$ and k,α defined above, the two kernel families have uniform joint endpoint bounds by Section 4. Hence for fixed retained (p,ε), the scaled $C^{k,\alpha }$ norm in (t,x) is at most Cε; the incident term separately is bounded by C$\varepsilon ^{2N+2}$. These bounds also hold after any fixed p derivative permitted by the endpoint class, and spatial derivatives do not consume the explicit ε factor. An ε derivative can consume that factor and can differentiate the cone coordinate; arbitrary ε derivatives are not being asserted to remain functions of this finite regularity.

Reflected-error proof 6. Full physical variables, without suppressing a power.

The physical kernel and error are exactly

$$\begin{split} K_N^{\rm phys}(t_{\rm phys},q,y) &=\varepsilon^{1-n}K_N\left(t_{\rm phys}/\varepsilon, (q-(p,0))/\varepsilon,y\right),\\ \mathcal R_N^{\rm phys}(t_{\rm phys},q,y) &=\varepsilon^{-n-1}\mathcal R_N\left(t_{\rm phys}/\varepsilon, (q-(p,0))/\varepsilon,y\right). \end{split} \tag{ER4}$$

The full wave operator contributes ε^{-2}; the full space/time Dirac mass contributes ε^{n+1} under composition. Therefore the point source remains exactly ρ(y)^{-1}δ(t_phys)δ_y(q)I_r. Nothing in the norm strengthening changes that density or the finite distribution identity.

Let β be a list of physical time/spatial derivatives of length $j\leq k$ and hold $y=(p,\varepsilon )$ fixed. On a physical compact set obtained by scaling a specified original compact set,

$$D_{t_{\rm phys},q}^{\beta}\mathcal R_N^{\rm phys} =\varepsilon^{-n-1-j}(D_{t,x}^{\beta}\mathcal R_N) \left(t_{\rm phys}/\varepsilon,(q-(p,0))/\varepsilon,p,\varepsilon\right), \quad \|D^\beta\mathcal R_N^{\rm phys}\|_\infty \leq C_\beta\varepsilon^{-n-j}. \tag{ER5}$$

For derivatives of order k their Hölder seminorm has the additional distance factor ε^{-α}, giving

$$[D^k\mathcal R_N^{\rm phys}]_\alpha \leq C\varepsilon^{-n-k-\alpha},\qquad \|D^j\mathcal R_{N,\rm incident}^{\rm phys}\|_\infty \leq C_j\varepsilon^{2N+1-n-j},\qquad [D^k\mathcal R_{N,\rm incident}^{\rm phys}]_\alpha \leq C\varepsilon^{2N+1-n-k-\alpha}. \tag{ER6}$$

For nonintegral a, k+α=N−(n−1)/2; for integral a the class is $C^{a-1,1}$. These are norms with the physical metric distance for a fixed source, not joint uniform bounds as ε→0 in physical source coordinates.

For completeness the exact ε derivative at fixed p,t_phys,q retains every term:

$$\partial_\varepsilon\mathcal R_N^{\rm phys} =\varepsilon^{-n-2} \big[-(n+1)\mathcal R_N-t\partial_t\mathcal R_N -\sum_jx_j\partial_{x_j}\mathcal R_N +\varepsilon\partial_\varepsilon\mathcal R_N\big] (t_{\rm phys}/\varepsilon,(q-(p,0))/\varepsilon,p,\varepsilon). \tag{ER7}$$

It is a classical identity when the required derivatives exist and always a distribution identity for ε>0. More generally, after j physical derivatives, write $F=D^\beta R_N$ and $\lambda =n+1+j$. For every h for which a classical claim is made, or distributionally without that restriction,

$$\partial_\varepsilon^h[\varepsilon^{-\lambda}F(t_{\rm phys}/\varepsilon, (q-(p,0))/\varepsilon,p,\varepsilon)] =\varepsilon^{-\lambda-h} \left[\prod_{v=0}^{h-1} (\varepsilon\partial_\varepsilon-t\partial_t -\sum_ix_i\partial_{x_i}-\lambda-v)\right]F (t_{\rm phys}/\varepsilon,(q-(p,0))/\varepsilon,p,\varepsilon). \tag{ER8}$$

Induction proves this identity by differentiating the prefactor and the three composed arguments; each new prefactor exponent is λ+v. All operators in the product differ only by scalar constants and therefore commute. At fixed physical q a p_i derivative also includes $-\varepsilon ^{-1}\partial _{x_i}$, in addition to the explicit $\partial _{p_i}$; thus these bounds are not silently transferred to physical source derivatives. This completes the original scaling comparison, the cone endpoint and the extra error factor, without removing any coefficient, density or exceptional integer endpoint.

8. The first geometric perturbation

The reflected coefficients respond to the first normal derivative of the metric. Fix one boundary point, choose tangential coordinates geodesic for the induced boundary metric, and orthonormalize the tangent frame at that point. Put $$A=-\tfrac12(\partial_{q_n}g_{jk}(0))_{j,k<n},\qquad K(u)=u^tAu,\qquad T=\operatorname{tr}A. \tag{R30}$$ The matrix $A$ is symmetric. Our shape operator is $-dN$, where $N$ points into $X$. Thus $T$ is the sum of the principal curvatures. A definition of mean curvature as their average replaces $T$ by $(n-1)$ times that average. For $n=1$, the tangential matrix is empty and $K=T=0$.

With the base point suppressed from the notation, Taylor expansion of the scaled metric in boundary-normal coordinates gives $$g_\varepsilon=\sum_jdq_j^2-2\varepsilon q_n\,dq'^t A\,dq' +O(\varepsilon^2). \tag{R31}$$ All tangential first derivatives vanish by the choice of the boundary coordinates; the normal cross coefficients vanish identically by (R4). Here and throughout the variation calculation, $O(\varepsilon^2)$ means $\varepsilon^2$ times a function smooth jointly in the spatial variables, the base point and $\varepsilon$. Spatial and base-point derivatives retain the two powers; an $\varepsilon$ derivative can consume one. Taylor’s integral remainder gives precisely this meaning, uniformly on each compact scaled set.

The change $$X'=q'-\varepsilon q_nAq',\qquad X_n=q_n+\tfrac{\varepsilon}{2}K(q') \tag{R32}$$ pulls the Euclidean metric back to (R31), through first order. Indeed the tangential differential square contributes $-2\varepsilon q_n\,dq'^tA\,dq'$ and $-2\varepsilon(Aq'\cdot dq')dq_n$; the normal differential square contributes the opposite cross term. The map fixes $e$, and sends the boundary to the graph $X_n=\varepsilon K(X')/2+O(\varepsilon^2)$. Smooth ODE dependence, the uniform transversality estimate in Section 3, and the inverse function theorem show that perturbing these metric and boundary jets by $O(\varepsilon^2)$ changes each ray, hitting time and inverse chart by $O(\varepsilon^2)$, with all fixed derivative bounds. Consequently the first variation can be computed in this Euclidean graph model.

Write the graph-model initial vector as $S=(u,z)$, with $z<-1$, and the reflected endpoint as $\Phi_\varepsilon(S)$. Both are expressed in the Euclidean $X$ coordinates, with center $\Pi=e$. The parameter in the following calculation is the unit travel-time parameter, so the total path length is $|S|$.

The unreflected path is $\Pi+\lambda S$. The hitting equation is

$$1+\lambda z=\frac{\varepsilon}{2}\lambda^2K(u).$$

At ε=0 its solution is $\lambda_0=-1/z$, and the derivative of the left minus right side with respect to λ is z. Hence the implicit function theorem gives a smooth hitting time and

$$\lambda z=-1+\frac{\varepsilon K(u)}{2z^2}+O(\varepsilon^2). \tag{V3}$$

At the hit point the inward unit normal equals $(\varepsilon Au/z,1)+O(\varepsilon^2)$. Its squared norm differs from 1 only in order ε². Therefore reflection of the tangent gives

$$\widetilde S=S-2(S\cdot N)N =\left(u-2\varepsilon Au, -z-2\varepsilon K(u)/z\right)+O(\varepsilon^2).$$

Adding the remaining segment, $\lambda S+(1-\lambda)\widetilde S+\Pi$, yields

$$\Phi_\varepsilon(S)= \left(u-2\varepsilon(1+1/z)Au, -1-z-\varepsilon K(u)(z^{-2}+2z^{-1})\right) +O(\varepsilon^2). \tag{V4}$$

Composing the endpoint with the fixed isometry $R(X',X_n)=(X',-1-X_n)$ leaves the Euclidean Laplacian unchanged and gives

$$R\Phi_\varepsilon(S)=S+\varepsilon\psi(S)+O(\varepsilon^2),\qquad \psi(S)=\left(-2(1+1/z)Au,(z^{-2}+2z^{-1})K(u)\right). \tag{V5}$$

No diagonalization of A was used, so the calculation applies to a general real symmetric second fundamental form.

9. The coordinate drift fixes the sign

Let $F_\varepsilon(S)=S+\varepsilon\psi(S)$, and put $D=\operatorname{div}\psi$ and $C_{jk}=\partial_j\psi_k+\partial_k\psi_j$. For a test function f in S coordinates, differentiate $f(F_\varepsilon^{-1}(X))$ twice. The inverse derivative is $I-\varepsilon D\psi+O(\varepsilon^2)$, so the scalar pullback of the positive Laplacian is

$$L_\varepsilon f =\Delta f-\varepsilon\sum_{j,k}C_{jk}\partial_j\partial_k f -\varepsilon\sum_k(\Delta\psi_k)\partial_k f +O(\varepsilon^2). \tag{V6}$$

The equivalent divergence form is

$$L_\varepsilon =\Delta-\varepsilon\sum_{j,k}\partial_j(C_{jk}\partial_k) +\varepsilon\sum_k(\partial_k D)\partial_k +O(\varepsilon^2). \tag{V7}$$

Indeed $\sum_j\partial_j C_{jk}=\Delta\psi_k+\partial_kD$. The last term in (V7) is the gradient of the divergence, not the vector whose components are $\Delta\psi_k$. Confusing these two tensors changes the operator.

The pullback of the actual spatial operator $-\Delta$ is

$$Q_\varepsilon=-L_\varepsilon =-\sum_{j,k}\partial_j(a^{jk}\partial_k) -\varepsilon\sum_k(\partial_kD)\partial_k +O(\varepsilon^2),\qquad a^{jk}=\delta_{jk}-\varepsilon C_{jk}. \tag{V8}$$

Let $\mathcal R=S\cdot\partial_S$. A direct differentiation of (V5) gives

$$D=-2(1+1/z)T-2(z^{-3}+z^{-2})K, \qquad \mathcal R D=\frac2z\left(T+\frac K{z^2}\right). \tag{V9}$$

The radial metric condition is $\sum_j C_{jk}S_j=0$. One elementary verification is to calculate $S\cdot\psi=-K/z$ and use $(D\psi+D\psi^T)S=\mathcal R\psi+\nabla(S\cdot\psi)-\psi$. Both tangential and normal components vanish. Consequently the radial drift h for (V8) is

$$h=-\varepsilon\mathcal R D+O(\varepsilon^2) =-\frac{2\varepsilon}{z}\left(T+\frac K{z^2}\right) +O(\varepsilon^2). \tag{V10}$$

This sign can also be obtained without rewriting the operator: the pullback volume Jacobian is $J=1+\varepsilon D+O(\varepsilon^2)$; therefore the radial Laplace–Beltrami drift is $\mathcal R\log J$, and the negative Laplacian has its negative.

The radial identity and the divergence expression must be used together. A coordinate change does not preserve a bare divergence formula relative to an arbitrarily fixed coordinate measure. The extra first-order term in (V8) is therefore present even though the Euclidean operator has no drift in its own coordinates.

10. Every first variation from one recurrence

For the pure Euclidean operator $-\Delta I_r$ in the graph model, define $w_\nu$ by $$v_0=I_r+\varepsilon w_0+O(\varepsilon^2),\qquad v_\nu=\varepsilon w_\nu+O(\varepsilon^2)\quad(\nu>0). \tag{R33}$$ The following formula holds for every integer $\nu\geq0$, with $(-1)!!=1$: $$w_\nu=(-1)^\nu 2^\nu(2\nu-1)!! \left[T+(2\nu+1)\frac{K(u)}{z^2}\right] \frac{1+z}{z^{2\nu+1}}I_r. \tag{V1}$$ The formula concerns each finite index. No convergence of the infinite Hadamard series is needed or asserted.

For $Q=-\partial_j(a^{jk}\partial_k)+b^k\partial_k+c$, with the radial metric condition, the product rule and the flat kernel identities give the coefficient of $E_{\nu-1}$ in $(\partial_t^2+Q)(v_\nu E_\nu)$ as

$$\nu v_\nu+\mathcal Rv_\nu-\frac12 h v_\nu, \qquad h=b\cdot S.$$

The coefficient contributed at the same kernel order by the preceding amplitude is $Qv_{\nu-1}$. Thus the transport system is

$$(\mathcal R+\nu)v_\nu-\frac12 h v_\nu=-Qv_{\nu-1}, \qquad v_{-1}=0. \tag{V11}$$

This also fixes the zeroth transport equation: $2\mathcal Rv_0=hv_0$. For the Euclidean graph model the incident zeroth amplitude is identically I and all positive incident amplitudes vanish. At ε=0 the hitting surface is z=−1; differentiating its shifted location introduces no first-order term because the limiting amplitudes are spatially constant. Hence $w_0|_{z=-1}=0$ and $w_\nu|_{z=-1}=0$ for ν>0.

Equation (V10) now gives

$$\mathcal Rw_0=-\frac1z\left(T+\frac K{z^2}\right),\qquad w_0|_{z=-1}=0.$$

Since $K/z^2$ is homogeneous of degree zero, radial integration gives $w_0=(T+K/z^2)(1+1/z)$, as asserted. There is a second direct check: D=0 on z=−1, so to first order the normalized transport solution is $J^{-1/2}$. At u=0,z=−2, its value is $(1-\varepsilon T)^{-1/2}=1+\varepsilon T/2+O(\varepsilon^2)$. This positive focusing variation is independent of the recurrence calculation.

For ν≥1, $h v_\nu=O(\varepsilon^2)$, and $Q=-\Delta+O(\varepsilon)$. Since the limiting positive amplitudes vanish, (V11) reduces to

$$(\mathcal R+\nu)w_\nu=\Delta w_{\nu-1},\qquad w_\nu|_{z=-1}=0. \tag{V12}$$

For ν=1 the right side is

$$\Delta w_0=\frac{2T(z+2)}{z^3} +\frac{6K(z+2)}{z^5},$$

and the solution is $w_1=-2(T+3K/z^2)(1+z)/z^3$. For the full induction set

$$F_\nu=\left[T+(2\nu+1)K/z^2\right] (1+z)z^{-2\nu-1}.$$

The identities $\Delta'K=2T$ and $u\cdot\partial_uK=2K$ give, for every integer ν≥1,

$$(\mathcal R+\nu)F_\nu =-\frac{\nu z+\nu+1}{z^{2\nu+1}} \left[T+(2\nu+1)K/z^2\right], \tag{V13}$$

$$\Delta F_{\nu-1} =2(2\nu-1)\frac{\nu z+\nu+1}{z^{2\nu+1}} \left[T+(2\nu+1)K/z^2\right] =-2(2\nu-1)(\mathcal R+\nu)F_\nu. \tag{V14}$$

For explicit verification, expand $F_\nu=T(z^{-2\nu-1}+z^{-2\nu}) +(2\nu+1)K(z^{-2\nu-3}+z^{-2\nu-2})$, differentiate each power twice in z, and add $2T$ times the coefficient of K. These operations produce (V13)–(V14) without any dimension-dependent simplification. Writing $w_\nu=c_\nu F_\nu$ therefore requires $c_0=1$ and $c_\nu=-2(2\nu-1)c_{\nu-1}$, whose solution is the coefficient in (V1). Every $F_\nu$ vanishes at z=−1. Uniqueness follows from the scalar radial ODE along each ray: the difference of two solutions is a constant times the inverse ν-th power of the radial parameter, and its value at the nonzero initial parameter is zero.

11. The diagonal derivative retains arbitrary lower-order terms

For the original operator (R1), including every smooth complex matrix $b^j,c$, the reflected amplitudes satisfy $$\left.\partial_\varepsilon V_\nu(e,p,\varepsilon) \right|_{\varepsilon=0} =(-1)^\nu 2^{-\nu-1}(2\nu-1)!!\,T(p)I_r, \qquad \nu\geq0. \tag{V2}$$ Here $x=e$ is held fixed in the scaled boundary-normal chart. In particular the right side is invariant under a change of the bundle frame at the base point. The value for $\nu=0$ is $T/2$, and the signs alternate for positive indices with our negative spatial operator convention.

After stretching coordinates, a smooth zeroth-order term has a factor ε². A smooth first-order matrix coefficient has the form $\varepsilon B^j+O(\varepsilon^2)$ on each fixed compact set, with constant matrices B^j at first order. A general smooth metric, in boundary-normal coordinates tangentially geodesic at the base point, has first jet

$$g_\varepsilon=\sum_j dq_j^2-2\varepsilon q_n\,dq'^T A\,dq' +O(\varepsilon^2). \tag{V15}$$

The explicit map

$$X'=q'-\varepsilon q_nAq',\qquad X_n=q_n+\frac{\varepsilon}{2}K(q') \tag{V16}$$

has Euclidean pullback metric (V15), carries q_n=0 to the graph, and fixes Π. Its Jacobian is $1-\varepsilon q_nT+O(\varepsilon^2)$. Thus the conversion between the divergence-form operator and the Laplace–Beltrami operator adds only a constant first-order term to this precision. More explicitly, $-\Delta_g=-\operatorname{div}(g^{-1}\nabla)+\varepsilon T\partial_n+O(\varepsilon^2)$, so an original coefficient B becomes $B^n-TI$ in the normal direction, with its other components unchanged at first order. All such terms can therefore be included in an effective constant collection $\mathsf B^j$.

Consider the matrix function

$$G_\varepsilon(X)=I+\frac{\varepsilon}{2} \sum_j\mathsf B^j(X_j-\Pi_j). \tag{V17}$$

It is invertible on every fixed compact set for sufficiently small ε, satisfies $G_\varepsilon(\Pi)=I$, and obeys

$$\left(-\Delta I+\varepsilon\sum_j\mathsf B^j\partial_j +O(\varepsilon^2)\right)(G_\varepsilon f) =G_\varepsilon(-\Delta f)+O(\varepsilon^2) \tag{V18}$$

as an equality of differential-operator coefficient expansions. To check it, use $2\partial_jG_\varepsilon=\varepsilon\mathsf B^j$, $\Delta G_\varepsilon=0$, and cancel the two first-order derivative terms. Products or commutators of distinct matrices first appear in order ε²; no commutativity, symmetry, reality or self-adjointness is needed.

Conjugation respects homogeneous Dirichlet data. With source normalization at Π it multiplies the first-order reflected zeroth amplitude by

$$\frac12\sum_j\mathsf B^j(\Phi_0(S)_j-\Pi_j) =\frac12\left(\sum_{j<n}\mathsf B^j u_j-\mathsf B^n(z+2)\right). \tag{V19}$$

This is affine in S, is harmonic, and equals zero at S=(0,−2). Its boundary value at z=−1 is precisely the first-order incident gauge value. Thus it contributes no forcing in (V12). For positive indices it also multiplies a limiting amplitude that is zero, so it cannot change their first variations. Equation (V19) proves both the off-diagonal qualification and the diagonal cancellation for every allowed lower-order operator.

The source-center metric in (V15) differs from the Euclidean one at first order. Its initial tangent is converted by the derivative of (V16) at Π; the normal tangent (0,−2) is unchanged, while any tangential reparameterization affects only first-order functions multiplied by ε. It therefore does not change (V2).

Finally, let $S_\varepsilon(x)=\Phi_\varepsilon^{-1}(x)$. At fixed x=Π, $S_0(\Pi)=(0,-2)$. The chain rule reads

$$\partial_\varepsilon V_\nu(\Pi,0) =\partial_\varepsilon v_\nu(S_0(\Pi),0) +d_Sv_\nu(S_0(\Pi),0)\, \partial_\varepsilon S_\varepsilon(\Pi)|_0.$$

The second term is zero: $v_0(\cdot,0)=I$ and $v_\nu(\cdot,0)=0$ for ν>0. Evaluating (V1) at u=0,z=−2 gives (V2).

For completeness, the density normalization also survives the graph reduction. The coordinate density of (R31) at the center is $1-\varepsilon T+O(\varepsilon^2)$, and the determinant of (R32) at the same point is $1-\varepsilon T+O(\varepsilon^2)$. Under a change of spatial variables, the point source picks up this Jacobian. The two factors in $\rho(e)^{-1}\delta_e$ cancel through first order, giving the Euclidean source normalization. No additional scalar first variation is hidden in (V2).