The permutations are perhaps best understood as permuting $L$, that is, we want $\mathbb E[x^TL'DL'^Ty]$ with $L'=P^TLP$. We can write $L=I-A$, where $I$ is the identity and $A$ has $1$s on the subdiagonal and $0$ elsewhere. The identity is invariant under the permutation, and where $A$ has $1$s in off-diagonal elements along a chain of indices from $1$ to $n$, $A'=P^TAP$ has them along any of the $n!$ possible chains. Thus
$$
\begin{align}
\mathbb E[x^TP^TIPDP^TI^TPy]&=x^TDy\;,
\\
\\
\mathbb E[x^TP^TIPDP^TA^TPy]
&=\mathbb E[x^TDA'^Ty]
\\
&=\frac1{n!}\sum_{\sigma\in S_n}\sum_{i=1}^nx_iD_{ii}y_{\sigma^{-1}(\sigma(i)+1)}
\\
&=\sum_{i=1}^nx_iD_{ii}\frac1{n-1}\sum_{j\ne i}y_j
\\
&=\frac n{n-1}x^TD\langle y\rangle-\frac1{n-1}x^TDy\;,
\\
\\
\mathbb E[x^TP^TAPDP^TI^TPy]
&=
\frac n{n-1}\langle x\rangle^TDy-\frac1{n-1}x^TDy\;,
\\
\\
\mathbb E[x^TP^TAPDP^TA^TPy]
&=
\mathbb E[x^TA'DA'^Ty]
\\
&=
\frac1{n!}\sum_{\sigma\in S_n}\sum_{i=1}^nx_{\sigma^{-1}(\sigma(i)+1)}D_{ii}y_{\sigma^{-1}(\sigma(i)+1)}
\\
&=\sum_{i=1}^n\frac1{n-1}\sum_{j\ne i}x_jD_{ii}y_j
\\
&=\sum_{j=1}^n\frac1{n-1}\sum_{i\ne j}x_jD_{ii}y_j
\\
&=\frac n{n-1}x^T\langle D\rangle y-\frac1{n-1}x^TDy\;,
\end{align}
$$
so
$$
\mathbb E[x^TP^TLPDP^TL^TPy]=\frac n{n-1}\left(x^TDy-\langle x\rangle^TDy-x^TD\langle y\rangle+x^T\langle D\rangle y\right)\;,
$$
where $\langle x\rangle$ and $\langle y\rangle$ are constant vectors whose entries are the average of the entries of $x$ and $y$, respectively, and $\langle D\rangle$ is a multiple of the identity whose diagonal entries are the average of the diagonal entries of $D$.