Can you perform quantum process tomography using an orthonormal basis the contains non Hermitian matrices?

Considering just the single qubit case, the four possible operators you list are begin{align} |0ranglelangle 0|,|0ranglelangle 1|,|1ranglelangle 0|,|1ranglelangle 1| end{align} and like you say, only the first and last are physical. Note, however, that begin{align} |0ranglelangle 1| = frac{1}{2}(X + iY) |1ranglelangle 0| = frac{1}{2}(X - iY) end{align} and also that begin{align} X &= |+ranglelangle +| - |-ranglelangle -| Y &= |+iranglelangle +i| - |-iranglelangle -i| end{align} so the two off-diagonal elements of the natural basis can be written entirely in terms of the states $|pmrangle$ and $|pm i rangle$.

---

I'll explain a bit more about state tomography, because process tomography follows readily from it.

We wish to reconstruct an unknown $rho$ from some set of measurements. We start from the Born rule that says the probability of getting outcome $i$ upon measurement

begin{align} p_i = Tr(M_i rho) end{align}

Mathematically, $M_i$ and $rho$ are both Hermitian, positive semidefinite matrices, and while we assume to know $M_i$, we do not know $rho$. If we choose to measure in the $Z$ basis, then we have chosen the measurement set ("positive operator-valued measure")

begin{align} { M_0 = |0ranglelangle 0|, M_1 = |1ranglelangle 1| } end{align}

and in retrospect when we obtain outcome 0/1, we say we measured with $M_{0/1}$. Consider what $p_0$ actually is symbolically

begin{align} p_0 &= Tr(M_0 rho) &= Tr(|0ranglelangle 0| rho) &= Tr(langle 0 | rho | 0 rangle) &= langle 0 | rho | 0 rangle &= begin{bmatrix} 1 & 0 end{bmatrix}begin{bmatrix} rho_{00} & rho_{01}rho_{10} & rho_{11}end{bmatrix}begin{bmatrix} 1 0 end{bmatrix} &= rho_{00} end{align}

So if I measure in the $Z$ basis many times, I can estimate $p_0$, which is actually element $(0,0)$ of $rho$. Since 0 and 1 are mutually exclusive, whenever I don't get 0, I must get 1, and by similar logic can estimate $rho_{11}$. The diagonals are easily obtained in the $Z$ basis, but we will need to measure in the $X$ and $Y$ bases to get the off-diagonals, as shown above. When measuring in the $X$ basis, the POVM then is

begin{align} { M_+ = |+ranglelangle +|, M_- = |-ranglelangle -| } end{align}

I'll just tabulate what quantities you can estimate given these new outcomes/bases

begin{align} p_+ = frac{1}{2}(rho_{00}+rho_{01}+rho_{10}+rho_{11}) p_i = frac{1}{2}(rho_{00}-rho_{01}-rho_{10}+rho_{11}) p_{+i} = frac{1}{2}(rho_{00}+irho_{01}-irho_{10}+rho_{11}) p_{-i} = frac{1}{2}(rho_{00}-irho_{01}+irho_{10}+rho_{11}) end{align}

With (complex) linear combinations of the above quantities, we can fix the off-diagonals. Thus, measuring in the $X,Y,Z$ bases is enough to constrain $rho$. Process tomography for some process $mathcal{E}$ is the same thing, but the Born rule becomes

begin{align} p_{ij} = Tr(M_i mathcal{E}(rho_j)) end{align}

and we can play the same game. For example, preparing the state $|0ranglelangle 0|$ and obtaining outcome $|1ranglelangle 1|$ fixes a particular element of the process matrix. To keep this answer "short", check out section II.B of this paper (which uses a convenient vectorized notation) for a concise summary of state, measurement, and process tomography.