Let $Y$ be a smooth second-countable Hausdorff manifold of dimension $n$, let $T>0$, let $I:=(-T,T)$, and write $D_t:=-i\partial_t$. Let $P:C^\infty(I\times Y)\to C^\infty(I\times Y)$ be a properly supported scalar differential operator of order $m\in\mathbb N$ with smooth real principal symbol $p_m(t,y;\tau,\eta)$ in the $D_t$-symbol convention. Let $\Gamma\subset T^*(I\times Y)\setminus 0$ be an open conic set. Let $a_m(t,y,\eta)$ denote the coefficient of $\tau^m$ in $p_m$. Assume that $a_m$ is elliptic on $\Gamma$, and define the normalized principal symbol $\widetilde p_m:=a_m^{-1}p_m$ on $\Gamma$. Assume that on $\Gamma$ this normalized symbol factors as \begin{align*} \widetilde p_m(t,y;\tau,\eta)=\prod_{\ell=1}^m\bigl(\tau-\lambda_\ell(t,y,\eta)\bigr). \end{align*} For each $1\le \ell\le m$, let $\lambda_\ell:I\times (T^*Y\setminus 0)\to\mathbb R$ be smooth and positively homogeneous of degree $1$ in $\eta$, and assume that the roots $\lambda_1,\dots,\lambda_m$ are pairwise distinct on the conic region under consideration. Assume also that $\{0\}\times Y$ is noncharacteristic for $P$ on $\Gamma_0:=\Gamma\cap T^*_{\{0\}\times Y}(I\times Y)$, meaning that $a_m(0,y,\eta)$ is elliptic on $\Gamma_0$.

Fix $(y_0,\eta_0)\in T^*Y\setminus 0$ such that $(0,y_0;\lambda_\ell(0,y_0,\eta_0),\eta_0)\in\Gamma$ for every $1\le \ell\le m$. Let $\Lambda\in\Psi^1_{\mathrm{prop}}(Y)$ be a properly supported positive elliptic tangential pseudodifferential operator whose principal symbol is $|\eta|$ on the chosen conic neighbourhood. Define the unrescaled Cauchy trace map $C:C^\infty(I\times Y)\to C^\infty(Y)^m$ by \begin{align*} C(u):=\bigl(u|_{t=0},\partial_tu|_{t=0},\dots,\partial_t^{m-1}u|_{t=0}\bigr). \end{align*} Define the rescaled companion map $U_\Lambda:C^\infty(I\times Y)\to C^\infty(Y)^m$ by \begin{align*} U_\Lambda(u):=\bigl(\Lambda^{m-1}u|_{t=0},\Lambda^{m-2}\partial_tu|_{t=0},\dots,\partial_t^{m-1}u|_{t=0}\bigr). \end{align*} Microlocally, both maps denote the induced trace maps on distributions for which the displayed Cauchy traces are defined modulo $C^\infty$ over the chosen conic set. In the theorem below the data $g_j$ are unrescaled Cauchy traces for $C(u)$; if rescaled data are used, the corresponding elliptic powers of $\Lambda$ are absorbed into the data operators.

For each branch define $q_\ell:T^*(I\times Y)\setminus 0\to\mathbb R$ by \begin{align*} q_\ell(t,y;\tau,\eta):=\tau-\lambda_\ell(t,y,\eta). \end{align*} Let $H_{q_\ell}$ be its Hamilton vector field with respect to the canonical symplectic form on $T^*(I\times Y)$. For sufficiently small time, let $\Phi_\ell^t$ denote the projected Hamilton flow obtained by lifting $(y,\eta)$ to $(0,y;\lambda_\ell(0,y,\eta),\eta)$, flowing by $H_{q_\ell}$ until the time coordinate is $t$, and projecting to $T^*Y\setminus 0$. Let $C_\ell\subset (T^*(I\times Y)\setminus 0)\times (T^*Y\setminus 0)$ denote the corresponding branch canonical relation from Cauchy data to spacetime covectors, and let $C_\ell^{\mathrm{ret}}\subset (T^*(I\times Y)\setminus 0)\times (T^*(I\times Y)\setminus 0)$ denote the corresponding retarded spacetime branch relation, restricted to source time $\sigma$ lying between $0$ and target time $t$.

Then, after possibly replacing $T$ by a smaller positive number, there exist an open conic neighbourhood $V\subset T^*Y\setminus 0$ of $(y_0,\eta_0)$, an open conic neighbourhood $W\subset T^*(I\times Y)\setminus 0$ of the union of the bicharacteristic flow-outs from $V$ along the Hamilton vector fields $H_{q_\ell}$, properly supported Fourier integral operators $E_j:\mathcal D'(Y)\to\mathcal D'(I\times Y)$ for $0\le j\le m-1$, and a properly supported Cauchy-retarded operator $R:\mathcal D'(I\times Y)\to\mathcal D'(I\times Y)$ with the following property. If $WF(g_j)\cap V$ is the prescribed microlocal support of the unrescaled Cauchy data $g_j\in\mathcal D'(Y)$ and $WF(f)\cap W$ is the prescribed microlocal support of the forcing term $f\in\mathcal D'(I\times Y)$, then $u:=\sum_{j=0}^{m-1}E_jg_j+Rf$ satisfies $Pu-f\in C^\infty$ microlocally on $W$, and $\partial_t^j u|_{t=0}-g_j\in C^\infty$ microlocally on $V$ for every $0\le j\le m-1$.

Each $E_j$ is a finite sum of properly supported Fourier integral operators associated with the branch canonical relations $C_\ell$. Let $(t,y)$ and $(\sigma,z)$ denote the target and source spacetime variables for the Schwartz kernel of $R$. The Schwartz kernel of $R$ is microlocally supported in the Cauchy-retarded region where $\sigma$ lies between $0$ and $t$, meaning $0\le \sigma\le t$ when $t\ge 0$ and $t\le \sigma\le 0$ when $t\le 0$, and its primed wave front relation is contained in the union of the retarded branch relations $C_\ell^{\mathrm{ret}}$. With unrescaled Cauchy traces, for every $s\in\mathbb R$ the parametrix maps $\prod_{j=0}^{m-1}H^{s-j-1/2}_{\mathrm{loc}}(Y)\times H^{s-m}_{\mathrm{loc}}(I\times Y)$ continuously, microlocally over $V\times W$, into $H^s_{\mathrm{loc}}(I\times Y)$ on $W$. If the $j$th Cauchy datum is instead $Q_j\partial_t^ju|_{t=0}$ with $Q_j\in\Psi^{r_j}_{\mathrm{prop}}(Y)$ elliptic on $V$, then the Sobolev exponent for that datum is shifted to $s-j-1/2-r_j$. The distribution produced by the parametrix is unique modulo $C^\infty$ microlocally on $W$. If two parametrices are constructed with the same microlocal cutoffs, then the primed wave front relation of the difference of their Schwartz kernels is disjoint from the relevant microlocal input-output relation over $V\times W$; equivalently, the two parametrices differ by a smoothing operator microlocally on $W$ for those cutoffs.