How to write the three-qubit GHZ state in the Pauli basis?

How can one write the GHZ state defined in Ket notation as $$|\psi\rangle= \frac{1}{\sqrt{2}} \left(|000\rangle + |111\rangle\right)$$, in terms of Pauli matrices $$\sigma_{1},\sigma_{2},\sigma_{3}$$?

A direct while brute-force way is by writing state $$|\psi\rangle$$ as density matrix $$\rho$$. Then by noticing that $$\rho=\frac{1}{2^3}\sum_{ijk}t_{ijk}\sigma_{i}\otimes\sigma_{j}\otimes\sigma_{k}$$ for $$i,j,k=1,2,3,4$$ stands for Pauli matrices and identity matrix, find all the coefficients $$t_{ijk}$$ by $$Tr(\rho\sigma_{i}\otimes\sigma_{j}\otimes\sigma_{k})$$. If you are familiar with this method, you might find it's nothing more than Quantum State Tomography.

But maybe for special states like GHZ state, we may not use this method and directly find the rule of those coefficients while I can't see it directly.

Summary: The expression you're looking for is:

$$\frac{1}{4} \left[ (III + IZZ + ZIZ + ZZI) + (XXX - XYY - YXY - YYX)\right]$$

where Pauli string notation like $$XYX$$ denotes $$\sigma_1 \otimes \sigma_2 \otimes \sigma_1$$, for example.

To start, we need to write the $$n$$-qubit GHZ state as an operator, namely

$$$$|\psi^{(n)}\rangle\langle\psi^{(n)} | = \frac{1}{2}\left( |0\rangle \langle 0|^{\otimes n} + |0\rangle \langle 1|^{\otimes n} + |1\rangle \langle 0|^{\otimes n} + |1\rangle \langle 1|^{\otimes n}\right) \tag{1}$$$$ where we use the notation $$$$A^{\otimes n} = \underbrace{A \otimes A \otimes \dots \otimes A}_{\text{n times}}\tag{2}$$$$ Taking the pauli matrices to be $$\{I, X, Y, Z\}$$, we will break this into its diagonal component and off-diagonal component. Looking at the general $$n$$-qubit case we can rewrite each term of $$|\psi\rangle \langle \psi|$$ as one of these: \begin{align} |0\rangle \langle 0|^{\otimes n} &= \frac{1}{2^n}(I + Z)^{\otimes n} \tag{3a-d}\\ |1\rangle \langle 1|^{\otimes n} &= \frac{1}{2^n}(I - Z)^{\otimes n}\\ |0\rangle \langle 1|^{\otimes n} &= \frac{1}{2^n}(X + iY)^{\otimes n}\\ |1\rangle \langle 0|^{\otimes n} &= \frac{1}{2^n}(X - iY)^{\otimes n} \end{align}

The terms on the right each look a bit like $$(a + b)^n$$ and $$(a - b)^n$$, which we would be able to expand using the Binomial theorem for $$a, b \in \mathbb{R}$$. However since we are dealing with non-commuting operators, we will have to settle for something that just looks like the binomial theorem and see how much we can simplify.

Permutations for expressing $$(A + B)^{\otimes n}$$

Consider a permutation operator $$S_\pi$$ for some permutation $$\pi: \{1, \dots, n\}\rightarrow \{1, \dots, n\}$$ which has the effect of rearranging the subsystems of its input according to a given permutation $$\pi$$. For example if $$(\pi(a), \pi(b), \pi(c)) = (b, c, a)$$ then the action of $$S_\pi$$ on a separable operator would be $$$$S_\pi (A \otimes B \otimes C) = B \otimes C \otimes A \tag{4}$$$$

Then for two operators $$A, B$$ with $$[A, B] \neq 0$$ in general we can write \begin{align} (A + B)^{\otimes n} = A^{\otimes n} + \sum_{\pi} S_\pi\left(A^{\otimes (n-1)} \otimes B\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)}\otimes B^{\otimes 2}\right) + \cdots \tag{5} \\ (A - B)^{\otimes n} = A^{\otimes n} - \sum_{\pi} S_\pi\left(A^{\otimes (n-1)} \otimes B\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)}\otimes B^{\otimes 2}\right) - \cdots \tag{6} \end{align}

We can combine (5)-(6) to find \begin{align} \frac{1}{2} &\left((A + B)^{\otimes n} + (A - B)^{\otimes n} \right) \\&= A^{\otimes n} + \sum_{\pi} S_\pi\left(A^{\otimes (n-2)} \otimes B^{\otimes 2}\right) + \sum_{\pi} S_\pi\left(A^{\otimes (n-4)} \otimes B^{\otimes 4}\right) + \cdots \\&= \sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(A^{\otimes (n-2t)} \otimes B^{\otimes 2t}\right) \tag{7} \end{align} where we observe that every term $$B^{\otimes k}$$ with odd $$k$$ has been eliminated. This lets us write the following:

\begin{align} |0\rangle \langle 0|^{\otimes n} + |1\rangle \langle 1|^{\otimes n} &= \frac{1}{2^{n}}\left[(I + Z)^{\otimes n} + (I - Z)^{\otimes n}\right] \\&= \frac{1}{2^{n-1}} \sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(I^{\otimes (n-2t)} \otimes Z^{\otimes 2t}\right) \tag{8} \end{align}

and$$^1$$ \begin{align} |0\rangle \langle 1|^{\otimes n} + |1\rangle \langle 0|^{\otimes n} &= \frac{1}{2^n} \left[(X + iY)^{\otimes n} + (X - iY)^{\otimes n} \right] \\&= \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} \sum_\pi S_\pi\left(X^{\otimes (n-2t)}\otimes (iY)^{\otimes 2t}\right) \\&= \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} (-1)^t \sum_\pi S_\pi\left(X^{\otimes (n-2t)}\otimes Y^{\otimes 2t}\right) \tag{9} \end{align}

with the final expression becoming the sum of these two parts: $$$$|\psi^{(n)}\rangle\langle\psi^{(n)} | = \frac{1}{2^{n-1}}\sum_{t=0}^{\lfloor n/2\rfloor} \left[ \sum_\pi S_\pi\left(I^{\otimes (n-2t)} \otimes Z^{\otimes 2t}\right) + (-1)^t \sum_\pi S_\pi\left(X^{\otimes (n-2t)} \otimes Y^{\otimes 2t}\right)\right] \tag{10}$$$$

Solution for $$n=3$$

Now we can write out 3-qubit case by just enumerating over Pauli strings with the correct parity. In the sum of (8) we will keep the terms $$III$$, $$IZZ$$, $$ZIZ$$, $$ZZI$$. In the sum of (9) we will keep the terms $$XXX$$, $$XYY$$, $$YXY$$, $$YYX$$. Plugging into (10) we get

\begin{align} |\psi^{(3)}\rangle\langle\psi^{(3)} | = \frac{1}{4} \left[ (III + IZZ + ZIZ + ZZI) + (XXX - XYY - YXY - YYX)\right] \tag{11} \end{align}

Direct Fidelity Estimation

Looking at Equations (10) and (11) we can see that as the system size grows, the GHZ state will continue to have support on a tiny subset of the $$4^n$$ Pauli strings: There are only terms containing $$Z$$ and $$I$$ together or $$X$$ and $$Y$$ together, and even among these a large chunk are missing for the parity constraints shown above.

One consequence of this is that you can estimate the fidelity of a GHZ on $$n$$ qubits using only $$O(n)$$ measurement configurations (Gühne, 2007) compared to something exponential in $$n$$ that you might expect. This is a special case of the more general idea of Direct Fidelity Estimation (DFE) (Flammia, 2011) which demonstrated how fidelity estimation could be improved based on how many of the $$4^n$$ Pauli strings were missing from the Pauli representation of your state. Of course, (Huang, 2020) brought this way down to $$O(1)$$ scaling using classical shadows but this is slightly less intuitive than the basic idea.

$$^1$$ This recovers Equation (5) in (Gühne, 2007)