We will derive the necessary and sufficient conditions for a local maximum of the functional
$$ g(L) := \sum_{P \in \mathcal{P}} z^T L L^T(P L L^T + I)^{-1} z = z^T L(L^T P L + I)^{-1} L^T z,$$
where $L$ is a real $n\times n$ matrix with $\mathrm{trace}(L L^T) = 1$ (you essentially proved the above equality in the comments, the equations simplify somewhat further by using this).
The main trick is to view $L$ as a vector in $\mathbb R ^{n^2}$ lying on the $(n^2-1)$-dimensional unit sphere. The special orthogonal group $SO(n^2)$ acts transitively on this sphere, so if we fix some vector $v \in \mathbb R ^{n^2}$ on the unit sphere then we can instead of $g$ maximize the functional on $SO(n^2)$ given by
$$ h(R):=\sum_{P \in \mathcal{P}} z^T [Rv]([Rv]^T P [Rv] + I)^{-1} [Rv]^T z,$$
where $R \in SO(n^2)$ and where the square brackets turn a column vector of length $n^2$ in a $n\times n$ matrix (say row by row starting at the upper left corner).
Since the Lie algebra of $SO(n^2)$ is the set of $n^2\times n^2$ skew-symmetric matrices (see for example 'Introduction to smooth manifolds' by John Lee), the necessary and sufficient conditions for a local maximum $R$ become
$$\left.\frac{d}{dt}\right|_{t=0}h(e^{tA}R)=0,$$
$$\left.\frac{d^2}{dt^2}\right|_{t=0}h(e^{tA}R)\le 0,$$
for any $n^2\times n^2$ skew-symmetric matrix $A$. We will work out the first order conditions and leave working out the second order conditions to the interested reader.
Let $\overline{E}_{ab}$ be a $n^2\times n^2$ elementary matrix (so zeros everywhere except at position $ab$ where it is $1$). Because the first order conditions are linear in $A$ we only need to consider $$A_{pqrs}:=\overline{E}_{((p-1)n+q)((r-1)n+s)}-\overline{E}_{((r-1)n+s)((p-1)n+q)}$$ in them, for $1\le p,q,r,s \le n$.
We can compute that $$[A_{pqrs}Rv]= E_{pr}[Rv]E_{sq}-E_{rp}[Rv]E_{qs}$$ (where $E_{ij}$ is a $n\times n$ elementary matrix).
We get
\begin{align*}
\left.\frac{d}{dt}\right|_{t=0}h(e^{tA_{pqrs}}R)=\sum_{P \in \mathcal{P}} & z^T [A_{pqrs}Rv]([Rv]^T P [Rv] + I)^{-1} [Rv]^T z + z^T[Rv]([Rv]^T P [Rv] + I)^{-1} [A_{pqrs}Rv]^T z
\\& - z^T[Rv]([Rv]^T P [Rv] + I)^{-1} ([A_{pqrs}Rv]^T P [Rv]+[Rv]^T P [A_{pqrs}Rv])([Rv]^T P [Rv] + I)^{-1} [Rv]^T z.
\end{align*}
Finally writing $L=[Rv]$ and using the above expression for $[A_{pqrs}Rv]$ we can simplify the above to
\begin{align*}
0=\sum_{P \in \mathcal{P}} z^T (I- L(L^T P L + I)^{-1} L^T P) (E_{pr}L E_{sq}-E_{rp} L E_{qs})(L^T P L + I)^{-1}L^T z.
\end{align*}