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$.