I am sure that since you are asking this question you probably already understand this, but for future & other's reference let me give a quick recap of what we are trying to achieve.
Quantum channels
Any process (in an open quantum system) is some map $\Lambda$ from a space of density matrices to a space of density matrices. I write a, because these spaces are not necessarily of the same dimension (for instance, the tracing out a subsystem does not preserve dimension). Any unitary transformation is such a map as well.
We generally write $\Lambda(\rho_{\mathrm{in}}) = \rho_{\mathrm{out}}$ when our map transforms $\rho_{\mathrm{in}}$ to $\rho_{\mathrm{out}}$. Furthermore, since we always expect any $\rho_{\mathrm{out}}$ to be an actual physical state (it must be positive semidefinite and must have trace $1$), we impose two constraints on $\Lambda$.
- Any map $\Lambda$ should be completely positive. This ensures that $\rho_{\mathrm{out}}$ is always positive semidefinite, even if it a subsystem of a larger whole. This constraint is often abbreviated as "CP".
- Any map $\Lambda$ should be trace preserving: $\mathrm{tr}\big[\Lambda(\rho)\big] = \mathrm{tr}\big[\rho\big]$, $\forall\rho$. This ensures that $\rho_{\mathrm{out}}$ always has unit trace. We abbreviate this constraint as "TP".
Any map $\Lambda$ that is both CP & TP = CPTP, we call a quantum channel. Sometimes we relax the TP constraint to include trace decreasing maps (consider for instance a measurement); some authors refer to these maps as the more general quantum operations.
Different representations of quantum channels
A quantum channel can be represented in different ways; I recap three here.
- The Kraus representation. Nielsen & Chuang refer to this as the operator-sum representation. In mathematical form:
\begin{equation}
\Lambda(\rho) = \sum_{k} A_{k} \rho A_{k}^{\dagger},
\end{equation}
where $\{A_{k}\}$ are known as the Kraus operators and $k$ is always the same or smaller than the system size $d = 4^{n}$. That is, for any valid map $\Lambda$ there can always be found a Kraus representation with at most $d$ operators. The CP constraint is automatically met here, the trace constraint reads: $\sum_{k} A_{k}^{\dagger} A_{k} \leq I$ (with equality for TP).
- The Choi matrix, which is a direct result of the Choi-Jamiolkowski isomorphism. Some intuition on what this is can be found in this previous answer. Consider the maximally entangled state $|\Omega \rangle = \sum_{\mathrm{i}}|\mathrm{i}\rangle \otimes |\mathrm{i}\rangle$, where $\{|\mathrm{i}\rangle\}$ forms a basis for the space on which $\rho$ acts. (Note that we thus have a maximally entangled state of twice as many qubits).
The Choi matrix is the state that we get when on one of these subsystems $\Lambda$ is applied (leaving the other subsystem intact):
\begin{equation}
\rho_{\mathrm{Choi}} = \big(\Lambda \otimes I\big) |\Omega\rangle\langle\Omega|.
\end{equation}
As the Choi matrix is a state, it must be positive semidefinite (corresonding the the CP constraint) and must be unity trace (corresponding to the TP constraint).
- The process- or $\chi$-matrix. We write our map as a double sum:
\begin{equation}
\Lambda(\rho) = \sum_{m,n} \chi_{mn}P_{m}\rho P_{n}^{\dagger},
\end{equation}
where $\{P_{m}\}$ & $\{P_{n}\}$ form a basis for the space of density matrices$^{1}$; we use the Pauli basis $\{I,X,Y,Z\}^{\otimes n}$ (thereby omitting the need for the $\dagger$ at $P_{n}$). The matrix $\chi$ now encapsulates all information of $\Lambda$; the CP constraint reads that $\chi$ must be positive semidefinite, and the trace constraint reads that $\sum_{m,n}\chi_{mn}P_{n}P_{m} \leq I$ (with equality for TP).
The goal for quantum process tomography is now to find a representation of an unknown channel $\Lambda$. We focus on the process matrix.
Standard QPT
Our goal is to find $\chi$ for an arbitrary quantum channel. We give ourselves only the power of inputting different input states $\rho_{\mathrm{in}}$, and measuring the output state $\rho_{\mathrm{out}}$ in different bases with measurement observables $\{M\}$.
We always measure in the Pauli basis, and using slight abusive notation we use the Pauli basis as input states as well$^{2}$. A measurement on $\rho_{out}$ in a basis denoted by $P_{j}$ with an input state $P_{i}$ then has expectation value $\lambda_{ij}$:
\begin{equation}
\begin{split}
\lambda_{ij} &= \mathrm{tr}\big[P_{j}\Lambda(P_{i})\big] \\
&= \mathrm{tr}\big[P_{j}\sum_{mn}\chi_{mn}P_{m}P_{i}P_{n}\big] \\
&= \sum_{mn}\chi_{mn} \mathrm{tr}\big[P_{j}P_{m}P_{i}P_{n}\big] \\
&= \sum_{mn} A_{(ij,mn)} \chi_{mn}.\\
\end{split}
\end{equation}
where $A_{(ij,mn)} = \mathrm{tr}\big[P_{j}P_{m}P_{i}P_{n}\big]$. If we now view all measurement outcomes $\{\lambda_{ij}\}$ as a vector $\overrightarrow{\lambda}$ and if we vectorize $\chi$ to $|\chi\rangle\rangle = \overrightarrow{\chi}$ we get a giant linear system of equations linking the measurement outcomes to the elements of $\chi$:
\begin{equation}
\overrightarrow{\lambda} = A \overrightarrow{\chi}.
\end{equation}
It is now our goal to solve for $\chi$.
Intermediary: some notes on the sets $\{P_{i}\}$ and $\{P_{j}\}$
The set of states from which we built $\{P_{i}\}$, known as the preparation set, needs to have (at first glance) every Pauli eigenstate for every qubit there is, resulting in $6^{n}$ different states. However, building all the Pauli matrices can be done by using any set of states that form a basis for the space of density matrices. A straightforward choice is to use $\{|0\rangle, |1\rangle, |+\rangle, |+i\rangle\}^{\otimes n}$ - both eigenstates of the $Z$ operator and the +1 eigenstates of the X and Y operator. This results in $4^{n}$ different input states.
The set of measurements from which we built all $\{P_{j}\}$ can be as simple as $\{X,Y,Z\}^{\otimes n}$; the "$I$"-measurements can be inferred from the outcomes of those measurements. I link a previous answer of mine on QST where I explain how to built all $\{P_{j}\}$'s from only these measurements; there I explain it in detail for $2$-qubit QST but the generalization to a higher number of qubits is very straightforward.
All in all, we thus need $4^{n} \times 3^{n} = 12^{n}$ different pairs of measurement operators and preparation states to perform QPT.
Solving for $\chi$
Solving for $\chi$ in our system of linear equations can be as straightforward as inverting $A$ ($A$ is indeed invertible). Moreover, by using the Pauli basis, $A$ is also unitary and even Hermitian so $\overrightarrow{\chi}$ is readily calculated as:
\begin{equation}
\overrightarrow{\chi} = A \overrightarrow{\lambda}.
\end{equation}
This, however, does not respect the CPTP constraints in any way.
Luckily, it can be shown that as long as all measurements that were performed had an actual outcome, the TP constraint on $\chi$ is automatically met when using this method.
However, the CP constraint is not automatically met; this means that the calculated $\chi$ might very well have negative eigenvalues. This stems from statistical noise on our estimates of $\lambda_{ij}$ which can be reduced by performing more repeated measurements. Note however that statistical noise will pretty much always persist (the statistical noise scales exponentially with the number of qubits considered in the QPT, so you need to repeat all measurements exponentially often to gain a constant fidelity. Methods to solve this problem are therefore needed.
A very straightforward but also less-then-ideal method of finding a positive semidefinite version of a non-CP $\chi$ is by taking a convex combination with the Identity process such that every eigenvalue of $\chi$ becomes non-negative. Let $\lambda_{\mathrm{min}} < 0$ be the smallest (largest negative) eigenvalue of $\chi$. Then the process matrix $\chi^{*}$ has only nonnegative eigenvalues:
\begin{equation}
\chi^{*} = \frac{1}{\mathrm{tr}[\chi]+2^{2n}|\lambda_{\mathrm{min}}|}\big(\chi + |\lambda_{\mathrm{min}}|I\big)
\end{equation}
The fraction before the sum is a renormalization constant. Of course there are other methods of bringing the eigenvalues to non-negative values, but those will most likely break the TP constraint (this method does in fact not).
(I would expect about 1-100 million repetitions per prepare-measure pair for $3$-qubit QPT will bring negative eigenvalues to something small enough for this method to give a okayish-fidelity. As I mentioned, QPT is hard.)
There exist more elaborate methods that I won't go into to much detail here. If you are familiar with semidefinite linear programming, solving the equation $\overrightarrow{\lambda} = A \overrightarrow{\chi}$ with $\chi$ subjected to our CPTP constraints is exactly a problem that can be optimized by this method. It should be noted that a proper optimization method will reduce the needed number of repeated measurements (often called shots) drastically when compared to the simple linear inversion mentioned above - I can therefore highly recommend optimalization in one way or another if you have a proper interest in QPT.
Another method can be found in this paper, where the authors make use of repeated projection onto the space of CP and TP maps, respectively. Proper proofs that this always converges to a proper minimum I have not yet seen, but I don't rule out the possibility either.
Another approach is to use the Choi-Jamiolkowski isomorphism, mentioned above regarding the Choi matrix. Here one would not optimize for a quantum channel, but for a quantum state (i.e. Quantum State Tomography). QST is much more popular and therefore many more optimization methods exist - I won't go into them here. It should be noted that this approach should be treated very carefully, as general quantum states do not always correspond to a proper quantum channel - which means that the optimization method used for the QST process outputs an estimate that is not a valid quantum channel.
Further reference or reading
My own MSc. thesis (Blatant self-promotion, please excuse me:)) can be found here, where I elaborate more on QPT in chapter 4. Chapter 3 might be a good read as an introduction to the terminology I use in chapter 4.
The text might be a bit convoluted at points but I feel that it introduces most of the basics. I have another text that I like better but I am not sure if I can distribute it; I will check. Furthermore, please feel free to ask me any subsequent questions.
Footnotes
- Note that using a different basis does in fact transform the $\chi$ matrix. We almost always use the Pauli basis though. Also, note that the Paulis are, in fact, not density matrices - but they do form a basis for the space of Hermitian matrices, of which the density matrices are a subset.
- Of course the Pauli operators are not valid density matrices since they are traceless, but using linearity we can combine the eigenstates of a Pauli operator to make that operator. If $\{|\psi_{+}\rangle\}$ & $\{|\psi_{-}\rangle\}$ are the $+1$- & $-1$ eigenstates of a Pauli operator $P$, we can combine them as such:
\begin{equation}
\Lambda(P) = \sum_{+}\Lambda(|\psi_{+}\rangle\langle \psi_{+}|) - \sum_{-}\Lambda(|\psi_{-}\rangle\langle \psi_{-}|).
\end{equation}
