I have found the following claim made very clearly at least once in the published literature (see below):

Let $P$ be a linear partial differential operator defined on an open set $\Omega \subset \mathbb{R}^{n+1}$, strictly hyperbolic with respect to the level surfaces of the first coordinate, which I will denote by $t$. Let $u$ be a distribution in $H^s_\mathrm{loc}(\Omega)$ for some $s \in \mathbb{R}$, such that $Pu=0$. Then the restriction (i.e. pullback) of $u$ to each of the level sets $\Sigma_{\tau} := \Omega \cap \{ t = \tau \}$ is itself in $H^s_\mathrm{loc}(\Sigma_\tau)$.

In fact, a more general statement may be found in Section 2 of:

Bao, Gang; Symes, William W., A trace theorem for solutions of linear partial differential equations, Math. Methods Appl. Sci. 14, No.8, 553-562 (1991). ZBL0754.35023.

There, the authors say that the statement boxed above follows from standard energy estimates, and refer to e.g. Taylor's book on pseudo-differential operators for the arguments. However, I meet the following difficulties when attempting to fill in all details:

• Ostensibly, the "standard energy estimates" discussed in Taylor's book (and elsewhere, e.g. in Hörmander's) pertain distributional solutions which (in the case $\Omega = \mathbb{R}^{n+1}$) are in spaces $\bigcap_{j=0}^{m-1} C^j([\tau_1,\tau_2]; H^{s-j}(\mathbb{R}^n))$, where $m$ is the order of $P$. I, on the other hand, am interested in distributions (locally) in $H^s$ in both time and space.
• In any case, as far as I am aware, the usual energy estimate controls the evolution in time of the sum of (appropriate) "spatial" Sobolev norms of (appropriate) $t$-derivatives of $u$ restricted to the level sets of $t$. Specifically, in the case $Pu=0$ the energy estimate says that the supremum of this quantity over a compact interval $[\tau_0,\tau_1]$ is bounded by a constant times its value at $\tau_0$. To prove that the restriction of $u$ to the $t=\tau_0$ surface is (locally) in $H^s$, presumably we would like to control the $H^s$ norm of this restriction and place an upper bound on it, but this does not seem to be the scenario in the energy inequality.

By the way, I am not doubting that $u$ can be pulled back to a spacelike hypersurface, as a distribution. I know that this is the case because $WF(u)$ is disjoint from the normal bundle of the hypersurface whenever $Pu$ is smooth.

##Update##

I might have the beginning of an argument, but this leads to a further question about Bochner–Sobolev spaces whose answer seems clear in the case of positive Sobolev order $s$, but not so much otherwise. I will explain schematically and in particular assuming for simplicity and definiteness that $P$ is of second order and that $\Omega = \mathbb{R}^{n+1}$ (so that also $\Sigma_\tau \cong \mathbb{R}^n$).

Namely, let $u \in H^s_\mathrm{loc}(\mathbb{R}^{n+1})$ be in the kernel of $P$, and let $\phi \in C_\mathrm{c}^\infty(\mathbb{R}^{n+1})$. Let $\psi \in C_\mathrm{c}^\infty(\mathbb{R}^{n+1})$ be such that $\psi = 1$ on a neighbourhood of $\mathrm{supp}{\phi}$. Then $P(\phi u) = [P,\phi] u = [P,\phi] (\psi u)$, and the commutator is a first-order differential operator. Hence, $\phi u$ solves the PDE with non-zero right-hand side but vanishing Cauchy data on $\Sigma_{T_0}$ for any sufficiently small $T_0$. Now, if $u$ were smooth then the standard energy estimates would say that for any $T_1 > T_0$ there exists a constant $C>0$ such that $$ (*) \quad \mathcal{E}_s(\tau; \phi u) \leq C \int_{T_0}^{T_1}\| [P,\phi] (\psi u) (t,\cdot) \|^2_{H^{s-1}(\mathbb{R}^n)} \, \mathrm{d}t $$ where the energy of order $s \in \mathbb{R}$ and at time $\tau$ is defined as $$\mathcal{E}_s(\tau; \phi u) = \| (\phi u)|_{\Sigma_\tau} \|^2_{H^s(\mathbb{R}^n)} + \| \partial_t (\phi u)|_{\Sigma_\tau} \|^2_{H^{s-1}(\mathbb{R}^n)}.$$ If we could provide a lower bound for the right-hand side of $(*)$, in terms of the norm of $\psi u$ in $H^s(\mathbb{R}^{n+1})$, then presumably we would be close to being done (by a suitable density argument from $C^\infty$). Now I would like to say something along these lines: $$ (**) \quad \int_{T_0}^{T_1}\| [P,\phi] (\psi u) (t,\cdot) \|^2_{H^{s-1}(\mathbb{R}^n)} \, \mathrm{d}t \leq C_1 \| [P,\phi] (\psi u) \|^2_{H^{s-1}(\mathbb{R}^{n+1})} \leq C_2 \| \psi u \|^2_{H^{s}(\mathbb{R}^{n+1})}$$ where the second inequality would be due to $[P,\phi]$ being differential of order $1$. But how can I justify the first inequality? It seems almost obvious in the case of integer $s \geq 1$ – we are simply ignoring those multi-indices with a non-zero entry corresponding to $t$, which would otherwise contribute to the $H^{s-1}(\mathbb{R}^{n+1})$ norm. For general $s \in \mathbb{R}$, I'm not sure. Notice that the left-hand side of $(**)$ is the square of the norm in the Bochner–Sobolev space $L^2([T_0,T_1]; H^{s-1}(\mathbb{R}^{n}))$.

#Update 2#

Perhaps proving the first inequality in $(**)$ is not so difficult when the Sobolev order $s$ is real and $s \geq 1$. For then it seems we can use the Fourier transform, as follows. Set $\ell := s-1$ and $g := [P,\phi](\psi u) \in C^\infty_\mathrm{c}(\mathbb{R}^{n+1})$, and denote partial Fourier transforms with respect to $t$ and the $x$-variables in $\mathbb{R}^n$ by $\mathcal{F}_t$ and $\mathcal{F}_x$ respectively, while the full Fourier transform of $g$ will be $\hat{g}$. Then, using that $\mathcal{F}_t \mathcal{F}_x g = \hat{g}$, and the Plancherel theorem: \begin{align*} \int_{T_0}^{T_1} \| g(t,\cdot) \|^2_{H^{\ell}(\mathbb{R}^n)} \, \mathrm{d}t &\leq \int_{\mathbb{R}} \| g(t,\cdot) \|^2_{H^{\ell}(\mathbb{R}^n)} \, \mathrm{d}t \\ &= \int_{\mathbb{R}} \left( \int_{\mathbb{R}^n} (1 + |k|^2)^\ell |[\mathcal{F}_x g](t,k)|^2 \, \mathrm{d}^n k\right) \mathrm{d}t \\ &= \int_{\mathbb{R}} \left( \int_{\mathbb{R}^n} (1 + |k|^2)^\ell |[\mathcal{F}_t^{-1} \hat{g}](t,k)|^2 \, \mathrm{d}^n k\right) \mathrm{d}t \\ &= \int_{\mathbb{R}^n} (1 + |k|^2)^\ell \left( \int_{\mathbb{R}} |[\mathcal{F}_t^{-1} \hat{g}](t,k)|^2 \, \mathrm{d}t \right) \mathrm{d}^n k \\ &= \int_{\mathbb{R}^n} (1 + |k|^2)^\ell \left( \int_{\mathbb{R}} |\hat{g}(\tau,k)|^2 \, \mathrm{d}\tau \right) \mathrm{d}^n k \\ &\leq \int_{\mathbb{R}^{n+1}} (1 + |\tau|^2 + |k|^2)^\ell|\hat{g}(\tau,k)|^2 \, \mathrm{d}\tau \, \mathrm{d}^n k \\ &= \| g \|^2_{H^{\ell}(\mathbb{R}^{n+1})}. \end{align*} Of course, I used the fact that $\ell \geq 0$ (i.e. $s \geq 1$) in the final inequality above.