In my research I have stumbled upon the following problem: Find all pairwise orthonormal triples of vectors $\mathbf{x},\mathbf{y},\mathbf{z} \in \mathbb{C}^n$ ($n \geq 4$) satisfying the three `eigenequations'
$(a\mathbf{A}+b\mathbf{B}-c\mathbf{C})\mathbf{x}=(c-b)bc \,\mathbf{x} \\ (a\mathbf{A}+b\mathbf{B}-c\mathbf{C})\mathbf{y}=(c-a)ca \,\mathbf{y} \\ (a\mathbf{A}+b\mathbf{B}-c\mathbf{C})\mathbf{z}=(a+b)ab \,\mathbf{z}$
where - and this is crucial - the $n\times n$ matrices $\mathbf{A},\mathbf{B},\mathbf{C}$ themselves depend on $\mathbf{x},\mathbf{y},\mathbf{z}$, respectively, being defined as
$A_{ij} = \delta_{ij}\sum_k D_{kj}|x_k|^2 - x_i D_{ij} \bar{x}_j \\ B_{ij} = \delta_{ij}\sum_k D_{kj}|y_k|^2 - y_i D_{ij} \bar{y}_j \\ C_{ij} = \delta_{ij}\sum_k D_{kj}|z_k|^2 - z_i D_{ij} \bar{z_j}$
where $\delta_{ij}$ is the Kronecker symbol and $D_{ij}$ ($i,j =1,\ldots,n$) are certain given coefficients satisfying $D_{ii} = 0$, $D_{ij} = D_{ji} > 0$ for $i \neq j$. Notice that the matrices are positive semi-definite and satisfy $\mathbf{A} \mathbf{x} = \mathbf{B} \mathbf{y} = \mathbf{C} \mathbf{z} = 0$ (I think their rank is actually $n-1$) as well as $\mathbf{y}^+\mathbf{C}\mathbf{y} = \mathbf{z}^+\mathbf{B}\mathbf{z}, \quad \mathbf{z}^+\mathbf{A}\mathbf{z} = \mathbf{x}^+\mathbf{C}\mathbf{x}, \quad \mathbf{x}^+\mathbf{B}\mathbf{x} = \mathbf{y}^+\mathbf{A}\mathbf{y}$.
Using orthogonality, one can easily show from the equations that all the `mixed terms' vanish, i.e. $\mathbf{y}^+\mathbf{A}\mathbf{z} = 0$ and so on.
As for the scalars $a,b,c$, they are assumed positive and can be shown to be equal to,
$a = \sqrt{\mathbf{y}^+\mathbf{C}\mathbf{y}} = \sqrt{\mathbf{z}^+\mathbf{B}\mathbf{z}}, \quad b = \sqrt{\mathbf{z}^+\mathbf{A}\mathbf{z}} = \sqrt{\mathbf{x}^+\mathbf{C}\mathbf{x}}, \quad c = \sqrt{\mathbf{x}^+\mathbf{B}\mathbf{x}} = \sqrt{\mathbf{y}^+\mathbf{A}\mathbf{y}}$.
I've been able to notice that any triple of distinct canonical basis vectors (with arbitrary phase factors) solves the above problem. In other words, taking
$x_j = e^{i\phi_1}\delta_{i_1 j}, \quad y_j = e^{i\phi_2}\delta_{i_2 j}, \quad z_j = e^{i\phi_3}\delta_{i_3 j}$
for any distinct $i_1,i_2,i_3 \in \{1,\ldots,n\}$ and any angles $\phi_1,\phi_2,\phi_3$ does the trick. Sometimes there exist other solutions, e.g. if $D_{ij} = 1 - \delta_{ij}$, then any orthonormal triple $\mathbf{x},\mathbf{y},\mathbf{z}$ will do. And so my question is: Assuming all $D_{ij}$'s are distinct (other than $D_{ii} = 0$ and $D_{ij} = D_{ji}$, that is), is it possible to show that the solutions mentioned above are the only ones?
Let me add that with all of the above in mind, one can rewrite the equations in a quite nice form
$(\frac{\mathbf{B}}{\sqrt{\mathbf{x}^+\mathbf{B}\mathbf{x}}}-\frac{\mathbf{C}}{\sqrt{\mathbf{x}^+\mathbf{C}\mathbf{x}}})\mathbf{x}=(\sqrt{\mathbf{x}^+\mathbf{B}\mathbf{x}}-\sqrt{\mathbf{x}^+\mathbf{C}\mathbf{x}}) \mathbf{x} \\ (\frac{\mathbf{A}}{\sqrt{\mathbf{y}^+\mathbf{A}\mathbf{y}}}-\frac{\mathbf{C}}{\sqrt{\mathbf{y}^+\mathbf{C}\mathbf{y}}})\mathbf{y}=(\sqrt{\mathbf{y}^+\mathbf{A}\mathbf{y}}-\sqrt{\mathbf{y}^+\mathbf{C}\mathbf{y}}) \mathbf{y} \\ (\frac{\mathbf{A}}{\sqrt{\mathbf{z}^+\mathbf{A}\mathbf{z}}}+\frac{\mathbf{B}}{\sqrt{\mathbf{z}^+\mathbf{B}\mathbf{z}}})\mathbf{z}=(\sqrt{\mathbf{z}^+\mathbf{A}\mathbf{z}}+\sqrt{\mathbf{z}^+\mathbf{B}\mathbf{z}}) \mathbf{z}$
But I've not been able to move much further. I tried playing for a while with determinants, but came back empty-handed. Any hint would be much appreciated, even for the real case.
EDIT: I've noticed that the original system of three `eigenequations' actually concerns one and the same matrix $a\mathbf{A}+b\mathbf{B}-c\mathbf{C}$ and have now included this fact in the question. Also some other minor errors are now corrected.