[proofplan] The proof is the local real-principal-type Hadamard parametrix construction in the clean paired-Lagrangian setting. First we identify the characteristic set of $P$ with the null cone of the Lorentzian metric and use the Hamilton flow of the scalar principal symbol to obtain the Lagrangian relation $\Lambda_\varepsilon$. Then the Duistermaat-Hörmander real-principal-type parametrix construction is applied microlocally in $\Gamma_\varepsilon$: the eikonal equation fixes the Lagrangian relation, the diagonal normalization fixes the principal symbol, and the transport equations determine all lower homogeneous amplitude terms. Borel summation gives a full amplitude, and the remaining kernels are microlocally smooth in $\Gamma_\varepsilon$; multiplying by proper cutoffs preserves the microlocal construction while making the operator properly supported. [/proofplan]

[step:Identify the characteristic geometry of the normally hyperbolic operator]The scalar principal symbol of $P$ on $T^*O$ is \begin{align*} p(x,\xi)=g_x^{-1}(\xi,\xi). \end{align*} Since $g$ is Lorentzian and $\xi\neq 0$, the differential of the quadratic form $\xi\mapsto g_x^{-1}(\xi,\xi)$ with respect to the fibre variable is not zero on the null cone. Hence $dp\neq 0$ on $\Sigma_O$, so $\Sigma_O$ is a smooth conic hypersurface of $T^*O\setminus 0$. The Hamilton vector field \begin{align*} H_p\in \mathfrak{X}(T^*O\setminus 0) \end{align*} is tangent to $\Sigma_O$ because $H_p p=0$ by [citetheorem:TEMP-66]. Thus the Hamilton flow $\Phi_t$ preserves $\Sigma_O$ wherever it is defined. The time orientation fixes which connected component of the null cone is future-directed; the sign condition $\varepsilon t>0$ selects the corresponding forward branch for $\varepsilon=+$ and backward branch for $\varepsilon=-$. The causal convexity of $O$ ensures that projected null bicharacteristic segments used in the definition of $\Lambda_\varepsilon$ do not leave $O$ and re-enter it, while geodesic convexity supplies a single local geometric branch near the characteristic diagonal. In a sufficiently small conic neighbourhood of the characteristic diagonal, this flowout is an embedded conic Lagrangian and the intersection with $N^*\Delta_O$ is clean with characteristic diagonal \begin{align*} \Lambda_\varepsilon\cap N^*\Delta_O=\{((x,\xi),(x,-\xi)):(x,\xi)\in\Sigma_O\}. \end{align*}[/step]

[guided]The first point is that the operator is of real principal type on the null cone. The principal symbol of $P$ is scalar, namely \begin{align*} \sigma_2(P)(x,\xi)=p(x,\xi)\operatorname{id}_{E_x}, \end{align*} where \begin{align*} p(x,\xi)=g_x^{-1}(\xi,\xi). \end{align*} Because $g_x^{-1}$ is a nondegenerate Lorentzian quadratic form on $T_x^*O$, the fibre derivative of $\xi\mapsto g_x^{-1}(\xi,\xi)$ is the nonzero linear functional $2g_x^{-1}(\xi,\cdot)$ whenever $\xi\neq 0$. Therefore $dp$ cannot vanish on the nonzero null cone \begin{align*} \Sigma_O=\{(x,\xi)\in T^*O\setminus 0:p(x,\xi)=0\}. \end{align*} This proves that $\Sigma_O$ is a smooth conic hypersurface of $T^*O\setminus 0$. The Hamilton vector field $H_p$ is the vector field generated by $p$ using the canonical symplectic form on $T^*O$. Since $H_p p=0$ by [citetheorem:TEMP-66], every integral curve of $H_p$ that starts in $\Sigma_O$ remains in $\Sigma_O$. Thus the flow relation generated by $H_p$ is a null bicharacteristic relation. The theorem assumes the geometric facts needed to avoid caustic and self-intersection complications in the conic region under consideration: $\Lambda_\varepsilon$ is embedded, conic, and Lagrangian in $\Gamma_\varepsilon$, and its intersection with $N^*\Delta_O$ is clean. These assumptions are exactly what permits the use of paired-Lagrangian calculus rather than a more singular multi-branch calculus.[/guided]

[step:Choose local phase functions and write the paired ansatz] By the assumed local parametrizability of $\Lambda_\varepsilon\setminus N^*\Delta_O$, equivalently by the local parametrization theorem for conic Lagrangians [citetheorem:TEMP-43] applied in each conic chart, every point of $\Lambda_\varepsilon\setminus N^*\Delta_O$ has a conic neighbourhood on which $\Lambda_\varepsilon$ is represented by a nondegenerate positively homogeneous phase function \begin{align*} \phi:W\times \Theta\to \mathbb{R}, \end{align*} where $W\subset O\times O$ is open and $\Theta\subset \mathbb{R}^N_0$ is open conic for some $N\in\mathbb{N}$. The associated critical set is \begin{align*} C_\phi:=\{(x,y,\theta)\in W\times \Theta:\partial_{\theta_j}\phi(x,y,\theta)=0\text{ for every }1\leq j\leq N\}. \end{align*} On this chart, the parametrized Lagrangian is \begin{align*} \{(x,y;\partial_x\phi(x,y,\theta),\partial_y\phi(x,y,\theta)):(x,y,\theta)\in C_\phi\}. \end{align*} Near $N^*\Delta_O$ we use the standard pseudodifferential diagonal phase in local coordinates on $O$. A conic partition of unity in $\Gamma_\varepsilon$ means a locally finite family of smooth degree-zero homogeneous cutoffs on $T^*(O\times O)\setminus 0$ whose supports lie in the chosen diagonal and flowout charts and whose sum is $1$ on $\Gamma_\varepsilon$ after shrinking to the microlocal region under construction. The paired ansatz is the distribution obtained by summing the diagonal conormal oscillatory representation and the oscillatory representations attached to the phase charts for $\Lambda_\varepsilon\setminus N^*\Delta_O$ after applying this partition. Its coefficients are classical symbol-valued amplitudes, namely asymptotic sums of homogeneous smooth sections of $E\boxtimes E^*$ on the relevant conic phase domains, tensored with the half-density bundle on $O\times O$ identified using the product density $d\mu(x)d\mu(y)$; this is the half-density convention used in the real-principal-type paired-Lagrangian calculus. The eikonal equation is exactly the statement that the $x$-covector of the phase lies on the characteristic set: \begin{align*} p(x,\partial_x\phi(x,y,\theta))=0 \end{align*} for every $(x,y,\theta)\in C_\phi$. This holds because the parametrized Lagrangian is contained in $\Lambda_\varepsilon$, and every point of $\Lambda_\varepsilon$ has first covector lying in $\Sigma_O$. [/step]

[step:Normalize the diagonal symbol so that the leading term gives the identity kernel] Let \begin{align*} I_O:C_c^\infty(O;E)\to C^\infty(O;E) \end{align*} denote the identity operator. Its Schwartz kernel is the diagonal distribution \begin{align*} \delta_{\Delta_O}\in \mathcal{D}'(O\times O;E\boxtimes E^*), \end{align*} defined using the density $d\mu$ and the fibre pairing $E_x^*\times E_x\to\mathbb{C}$ by the identity \begin{align*} I_O f=f \end{align*} for every $f\in C_c^\infty(O;E)$. In the diagonal pseudodifferential part of the ansatz, choose the principal symbol $q_{-2}$ of the parametrix on the elliptic part of $T^*O\setminus 0$ by \begin{align*} q_{-2}(x,\xi)=p(x,\xi)^{-1}\operatorname{id}_{E_x} \end{align*} where $p(x,\xi)\neq 0$. Along the characteristic diagonal, use the standard $\varepsilon$-side boundary-value normalization for the inverse $p(x,\xi)^{-1}\operatorname{id}_{E_x}$. In the paired-Lagrangian principal-symbol exact sequence for the clean pair $(N^*\Delta_O,\Lambda_\varepsilon)$, this means that the homogeneous degree $-2$ diagonal coefficient is the normalized boundary value whose product with $p(x,\xi)\operatorname{id}_{E_x}$ is $\operatorname{id}_{E_x}$ modulo the flowout symbol term. With the half-density factor fixed by $d\mu(x)d\mu(y)$, this is precisely the principal symbol of the identity kernel $\delta_{\Delta_O}$. Thus the leading diagonal coefficient is fixed by the condition that the order-zero part of $PG_O^\varepsilon$ equals the identity microlocally near $N^*\Delta_O$. This normalization also fixes the initial data for the transport equations along the selected branch of $\Lambda_\varepsilon$ issuing from the characteristic diagonal. Therefore no additional homogeneous solution term may be added on that component of $\Lambda_\varepsilon$ without changing either the identity normalization or the prescribed transport initial data. [/step]

[step:Solve the transport equations along the null flowout]Let \begin{align*} K_\phi\in \mathcal{D}'(W;E\boxtimes E^*) \end{align*} denote the local oscillatory kernel distribution on the open set $W\subset O\times O$ defined by the phase chart $\phi$ and by the oscillatory integral \begin{align*} K_\phi(x,y)=\int_{\Theta} e^{i\phi(x,y,\theta)}a(x,y,\theta)\,d\mathcal{L}^N(\theta). \end{align*} Here \begin{align*} a:W\times \Theta\to E\boxtimes E^* \end{align*} is a classical symbol-valued amplitude in the half-density convention, and the displayed formula is understood after choosing a local trivialization of $E\boxtimes E^*$ over $W$. Apply $P$ in the left variable $x$ to this local kernel. The principal part obtained by differentiating the exponential is multiplication by \begin{align*} p(x,\partial_x\phi(x,y,\theta)). \end{align*} This term vanishes on $C_\phi$ by the eikonal equation. The next homogeneous term is a first-order transport operator along the Hamilton vector field $H_p$ restricted to $\Lambda_\varepsilon$, together with the subprincipal symbol of $P$ and the natural connection term induced by the chosen local trivialization of $E$ and the half-density bundle. The leading transport equation, with the initial value fixed by the diagonal identity normalization, determines the principal amplitude on the chosen component of $\Lambda_\varepsilon$. Inductively, after homogeneous terms of degrees higher than a given degree have been chosen, the coefficient of the next degree gives an inhomogeneous linear transport equation along the same null bicharacteristics. Since the bicharacteristic segments under consideration lie in the geodesically convex set $O$ and the initial data are prescribed on the characteristic diagonal, each such transport equation has a unique smooth solution on the selected local branch. Borel summation of the resulting homogeneous symbols gives a classical amplitude $a$ whose left residual has symbol vanishing to infinite order in the conic region. The right action of $P$ on the $y$ variable is not imposed by assuming formal self-adjointness of $P$. Instead, it is the transpose-kernel action specified in the statement using the $E^*$ pairing and the density $d\mu$. The repaired hypotheses include the two-sided paired-parametrix compatibility theorem for this convention: the left transport system fixes the left residual, the transpose right transport system fixes the right residual, and the residual correction theorem removes the remaining finite-order errors without changing the already normalized principal symbols. Hence the resulting paired kernel is a two-sided microlocal parametrix.[/step]

[guided]Here is the core Hadamard calculation. In a phase chart away from the diagonal, define \begin{align*} K_\phi\in \mathcal{D}'(W;E\boxtimes E^*) \end{align*} to be the local oscillatory kernel distribution represented, after choosing a local trivialization of $E\boxtimes E^*$ over $W$, by \begin{align*} K_\phi(x,y)=\int_{\Theta} e^{i\phi(x,y,\theta)}a(x,y,\theta)\,d\mathcal{L}^N(\theta), \end{align*} where $\Theta\subset\mathbb{R}^N_0$ is open conic and \begin{align*} a:W\times \Theta\to E\boxtimes E^* \end{align*} is a classical amplitude. The measure in the oscillatory integral is explicitly the $N$-dimensional Lebesgue measure $\mathcal{L}^N$ on the phase variable. When the second-order operator $P$ acts in the $x$ variable, the highest-order derivatives fall on the exponential factor. This produces the scalar factor \begin{align*} p(x,\partial_x\phi(x,y,\theta)) \end{align*} multiplying the leading amplitude. Therefore the leading obstruction to solving $PK=0$ away from the diagonal is precisely the eikonal equation \begin{align*} p(x,\partial_x\phi(x,y,\theta))=0 \end{align*} on $C_\phi$. This equation holds because the phase parametrizes $\Lambda_\varepsilon$, and $\Lambda_\varepsilon$ is contained in the null Hamilton flowout. Once the eikonal term vanishes, the next coefficient is no longer algebraic; it is a transport equation. Its differentiation direction is the Hamilton direction $H_p$ along $\Lambda_\varepsilon$, and its lower-order coefficients come from the subprincipal part of $P$, the bundle action on $E$, and the half-density convention. The initial condition is not arbitrary: at the characteristic diagonal, the diagonal pseudodifferential normalization requires that applying $P$ produce the identity kernel. This fixes the principal amplitude. After the principal amplitude is fixed, the lower-order amplitudes are obtained recursively. Suppose all homogeneous amplitude components above degree $r$ have already been chosen so that the corresponding residual vanishes through those degrees. The degree $r$ residual is a known smooth inhomogeneous term. Setting the next amplitude component to cancel it gives a linear transport equation along the same null bicharacteristics. Because the relevant bicharacteristic segment stays inside the chosen geodesically convex neighbourhood $O$, and because the initial value is specified at the characteristic diagonal, the transport equation is an ordinary linear differential equation along a smooth Hamilton integral curve with smooth coefficients. The local existence and uniqueness theorem for linear ordinary differential equations therefore gives a unique smooth solution on the local branch. Solving these equations for all homogeneous degrees gives a formal classical symbol, meaning a sequence of smooth homogeneous amplitude terms with decreasing degrees. The classical Borel summation theorem for symbols gives an actual classical amplitude with that prescribed asymptotic expansion. The remaining left residual has all homogeneous symbol coefficients equal to zero in the conic region; by the microlocal smoothing criterion for oscillatory kernels with symbols vanishing to infinite order, this residual is microlocally smooth there. Finally, the theorem requires both $PG_O^\varepsilon-I$ and $G_O^\varepsilon P-I$ to be microlocally smooth in $\Gamma_\varepsilon$. The right-sided statement is obtained by the same transport construction for the right action of $P$ on the $y$ variable. The density $d\mu$ and the pairing between $E^*$ and $E$ are exactly what make this right action well-defined at the level of kernels. Thus the amplitude is fixed, modulo smoothing amplitudes, by the diagonal normalization and the two systems of transport equations.[/guided]

[step:Invoke the clean paired-Lagrangian parametrix theorem] We use the standard real-principal-type paired-Lagrangian parametrix theorem in the vector-bundle half-density convention. In this form, the theorem applies to a differential operator with scalar real-principal-type principal symbol, an embedded conic Hamilton flowout $\Lambda$ from the characteristic diagonal, a clean pair $(N^*\Delta,\Lambda)$ locally represented by nondegenerate homogeneous phase functions, the diagonal principal-symbol exact sequence with the chosen boundary-value normalization, and compatible left and transpose-right transport systems with residual correction. Its conclusion is that the eikonal equation, diagonal identity normalization, and left and right transport equations constructed above produce a paired Lagrangian kernel $K$ associated with $(N^*\Delta,\Lambda)$ such that both left and right compositions with the operator equal the identity modulo kernels whose primed wave front sets are disjoint from the prescribed conic neighbourhood. The hypotheses of this theorem are satisfied here. The symbol $p$ is real and of principal type on $\Sigma_O$ by the first step. The flowout $\Lambda_\varepsilon$ is embedded conic Lagrangian in $\Gamma_\varepsilon$ by hypothesis. The pair $(N^*\Delta_O,\Lambda_\varepsilon)$ intersects cleanly by hypothesis. The local nondegenerate homogeneous phase parametrizations required to write the oscillatory representatives exist by hypothesis, and are compatible with the general parametrization result [citetheorem:TEMP-43]. Therefore the theorem supplies a kernel \begin{align*} K_{0,\varepsilon}\in \mathcal{D}'(O\times O;E\boxtimes E^*) \end{align*} which is microlocally a paired Lagrangian distribution associated with $(N^*\Delta_O,\Lambda_\varepsilon)$ in $\Gamma_\varepsilon$ and satisfies \begin{align*} \operatorname{WF}'(P K_{0,\varepsilon}-\delta_{\Delta_O})\cap\Gamma_\varepsilon=\varnothing \end{align*} and \begin{align*} \operatorname{WF}'(K_{0,\varepsilon} P-\delta_{\Delta_O})\cap\Gamma_\varepsilon=\varnothing. \end{align*} Here $P K_{0,\varepsilon}$ means that $P$ acts in the left $O$ variable of the kernel, while $K_{0,\varepsilon}P$ means the corresponding right kernel action defined using $d\mu$ and the $E^*$ pairing. [/step]

[step:Insert proper cutoffs without changing the microlocal construction] Choose an exhaustion of $O$ by relatively compact open sets $O_j\subset O$ with $\overline{O_j}\subset O_{j+1}$ and $\bigcup_{j=1}^{\infty}O_j=O$. Choose smooth functions \begin{align*} \rho_j:O\times O\to [0,1] \end{align*} with support contained in $O_{j+1}\times O_{j+1}$, equal to $1$ on $O_j\times O_j$, and such that the differences $\rho_{j+1}-\rho_j$ form a locally finite proper-support decomposition on $O\times O$. For each compact conic subset $C\subset\Gamma_\varepsilon$, choose $j$ so large that the base projection of $C$ is contained in $O_j\times O_j$; then $\rho_j=1$ on an open neighbourhood of that base projection. Apply the preceding microlocal construction on each conic part of $\Gamma_\varepsilon$ lying over $O_j\times O_j$, and sum the resulting locally finite series of cutoff kernels. Denote the resulting distribution by \begin{align*} K_O^\varepsilon\in \mathcal{D}'(O\times O;E\boxtimes E^*). \end{align*} Multiplication by each smooth cutoff does not enlarge wave front sets, and local finiteness permits wave front statements to be checked on compact conic subsets. If $C\subset\Gamma_\varepsilon$ is compact conic and $\rho_j=1$ near the base projection of $C$, then every derivative of $\rho_j$ vanishes near that base projection. Therefore the extra terms produced when $P$ differentiates the cutoff-supported kernel have primed wave front set disjoint from $C$. Since every point of $\Gamma_\varepsilon$ lies in such a compact conic subset after shrinking to a small conic neighbourhood, the cutoff replacement changes the local parametrix only by residual kernels whose primed wave front sets are disjoint from all of $\Gamma_\varepsilon$. The exhaustion construction gives support proper over both factors, so $K_O^\varepsilon$ defines a properly supported continuous linear operator \begin{align*} G_O^\varepsilon:C_c^\infty(O;E)\to \mathcal{D}'(O;E). \end{align*} The preceding microlocal equalities are therefore preserved: \begin{align*} \operatorname{WF}'(PG_O^\varepsilon-I)\cap\Gamma_\varepsilon=\varnothing \end{align*} and \begin{align*} \operatorname{WF}'(G_O^\varepsilon P-I)\cap\Gamma_\varepsilon=\varnothing. \end{align*} The same construction also preserves the paired-Lagrangian description in $\Gamma_\varepsilon$, the Fourier integral description away from $N^*\Delta_O$, the pseudodifferential diagonal description away from $\Lambda_\varepsilon$, the eikonal equation in each phase chart, and the transport determination of the full amplitude modulo smoothing amplitudes. Thus $G_O^\varepsilon$ has all asserted properties, and no support condition beyond proper support has been used or obtained. [/step]