Consider the maximally entangled state in $d$ dimensions, $|\Psi\rangle:= \frac{1}{\sqrt{d}} \sum_{i=0}^{d-1} |i,i\rangle^{AB}$, where $|i\rangle^{AB} := |i\rangle^{A}\otimes|i\rangle^{B}$ and $\{|i\rangle^{A}\}_{i=0}^{d-1}$ and $\{|i\rangle^{B}\}_{i=0}^{d-1}$ are the computational bases in Hilbert space $\mathcal{H}_A$ and $\mathcal{H}_B$, respectively.
It is known (I think the construction was found by Wooters) that if we choose bases $|c_j\rangle^{A} := \frac{1}{\sqrt{d}} \sum_{k=0}^{d-1} \omega^{jk} |k\rangle$ and $|c_j\rangle^{B} := \frac{1}{\sqrt{d}} \sum_{k=0}^{d-1} \overline{\omega}^{jk} |k\rangle$, where $\omega := e^{\frac{2\pi i}{d}}$ and $\overline{.}$ denotes the complex conjugate, that the maximally mixed state has the nice form $|\Psi\rangle = \frac{1}{\sqrt{d}}\sum_{j=0}^{d-1} |c_j\rangle^{A} \otimes |c_j\rangle^{B}$.
This can be shown by direct calculation;
Representing the 'old' basis vectors by the new ones and inserting into the maximally entangled state yields
$|\Psi\rangle = \frac{1}{d\sqrt{d}} \sum_{i,j,k=0}^{d-1} \omega^{(k-j)i}|c_j\rangle^{A}\otimes|c_k\rangle^{B}$.
Now we observe that $\sum_{i=0}^{d-1} \omega^{i(k-j)} = \begin{cases} 0, & \text{if } k\neq l\\d, & \text{if } k=l\end{cases}$ (since $\omega$ is a $d$-th root of $1$), which leads to $|\Psi\rangle = \frac{1}{\sqrt{d}} \sum_{j=0}^{d-1} |c_j\rangle^{A}\otimes|c_j\rangle^{B}$.
I was told that a similar result holds for arbitrary orthonormal bases in system $A$ and $B$, as long as the coefficients of the basis in $B$ are the complex conjugate of the coefficients in $A$. For some reason, to me, this isn't as obvious as it might be, so I tried to prove this statement. However, I'm having troubles showing that this indeed is true.
But, first things first; consider the orthonormal bases $\{|b_j\rangle^{A}\}_{j=0}^{d-1}$ and $\{|b_j\rangle^{B}\}_{j=0}^{d-1}$. We can express the 'old' basis vectors by the new ones
$|i\rangle^{A} = \sum_{j=0}^{d-1} \beta_{ij} |b_j\rangle^{A}$ with $\sum_{j=0}^{d-1} |\beta_{ij}|^2 = 1$ and $0 = \langle i | j\rangle^{A} = \sum_{k,l=0}^{d-1} \overline{\beta_{ik}} \beta_{il} \langle c_k|c_l\rangle^{A}=\sum_{k=0}^{d-1} \overline{\beta_{ik}} \beta_{ik}$,
$|i\rangle^{B} = \sum_{j=0}^{d-1} \gamma_{ij} |b_j\rangle^{B}$ with $\sum_{j=0}^{d-1} |\gamma_{ij}|^2 = 1$ and $0 = \langle i | j\rangle^{B} = \sum_{k,l=0}^{d-1} \overline{\gamma_{ik}} \gamma_{il} \langle c_k|c_l\rangle^{B}=\sum_{k=0}^{d-1} \overline{\gamma_{ik}} \gamma_{ik}$.
Inserting this yields
$|\Psi\rangle = \frac{1}{\sqrt{d}} \sum_{j,k=0}^{d-1} \sum_{i=0}^{d-1} \beta_{ij}\gamma_{ik} |b_j\rangle^{A}\otimes |b_k\rangle^{B}$.
If the claim is indeed true, there must be a way to show that $\sum_{i=0}^{d-1} \beta_{ij}\gamma_{ik} = \delta_{jk}$ if and only if $\gamma_{ij} = \overline{\beta_{ij}}$. So, we have to prove two directions.
Assume $\gamma_{ij} = \overline{\beta_{ij}}$.Then for $j=k:$ $\sum_{i=0}^{d-1} \beta_{ij}\overline{\beta_{ij}} = \sum_{i=0}^{d-1} |\beta_{ij}|^2 = 1$, since the basis vectors are assumed to be normalized. For $j \neq k:~ \sum_{i=0}^{d-1} \beta_{ij}\overline{\beta_{ij}} = 0$, as the basis vectors are assumed to be orthogonal. So, the reverse direction is easy.
Assume $\sum_{i=0}^{d-1} \beta_{ij}\gamma_{ik} = \delta_{jk}$. Then we obtain $d^2-d$ equations with RHS $0$ and $d$ equations with RHS $1$, altogether a system of $d^2$ equations for $d^2$ unknowns $\gamma_{ik}$ (the $\beta_{ij}$ are fixed, as we fix the first basis and try to find the second one).
I'm having troubles finishing the second direction, as it isn't obvious to me that the only solution to this system of equations is $\overline{\gamma_{ik}} = \beta_{ik}$. Am I missing any information? Is there any trick?
Any suggestions and hints are highly appreciated!
