Let $(M,g)$ be a globally hyperbolic Lorentzian manifold with Lorentzian volume density $dV$ and $\Sigma$ a Cauchy hypersurface, i.e. each inextendible causal curve hits $\Sigma$ exactly once. Furthermore let $P: C^\infty(M) \rightarrow C^\infty(M)$ be a formally self-adjoint wave operator, that is a $2^{nd}$-order differential operator with principal symbol given by the metric $g$ and $$(P\varphi,\psi)_{L^2(dV)} = (\varphi,P\psi)_{L^2(dV)}, \qquad \varphi,\psi \in C^\infty_c(M).$$ For each fixed $q\in M$ consider the Cauchy problem $$\left\{ \begin{array}{cl} Pu(\cdot,q) = 0 \qquad & \text{on} \enspace M \\[2mm] u(\cdot,q) = u_0(\cdot,q) & \text{on} \enspace \Sigma \\[2mm] \nabla_\nu u(\cdot,q) = u_1(\cdot,q) & \text{on} \enspace \Sigma \end{array} \right., \qquad (*)$$ which by well-posedness determines the function $u$ on $M\times M$ uniquely. Here $\nabla_\nu$ denotes the covariant normal derivative along $\Sigma$. My question now is:

When I have symmetry of the initial data on $\Sigma$ in the sense $$u(\sigma,q) = u(q,\sigma) \quad \text{and} \quad \nabla_\nu^{(1)} u(\sigma,q) = \nabla_\nu^{(2)} u(q,\sigma), \qquad \sigma \in \Sigma,~q \in M, \quad (**)$$ does this imply symmetry of $u$, i.e. $u(p,q)=u(q,p)$ for all $p,q\in M$? Btw $\nabla_\nu^{(1)}, \nabla_\nu^{(2)}$ stand for differentiating w.r.t the first and the second entry, respectively.

An equivalent question would be whether solutions of $(*)$ are solutions of the corresponding Cauchy problem for $Pu(q,\cdot)=0$ with initial data given by $(**)$.

All cases I checked, where one can give a certain solution operator propagating the solution from the initial data, symmetry directly holds true, which makes me quite optimistic. Thanks!