I posted (nearly) this on MSE, so if it doen't belong here plese remove.
Let $\Phi\in\mathcal{L}(\mathbb{C}^{n \times n})$ and suppose that for any unitary
conjugation operator $\mathcal{U}\in\mathcal{L}(\mathbb{C}^{n \times n})$,
$\mathcal{U}\Phi=\Phi\mathcal{U}$. Let $\{e_1,e_2,\dots,e_n\}$ denote the standard
orthonormal basis for $\mathbb{C}^n$, let $e_{ij}$ denote the $n\times n$ matrix with $1$ in the $i$th row
$j$th column and zeros elswhere, and let $E_{ij}=\Phi(e_{ij})$. Claim,
there exists $r,s\in\mathbb{C}$ such that,
\begin{equation}
E_{ij}=\begin{cases} re_{ij} & \text{ if } i\neq j, \text{ and }\\
%
re_{ij}+sI & \text{ if } i= j,
\end{cases}\tag{1}
\end{equation}
so that for $B\in\mathbb{C}^{n\times n}$,
$\Phi(B)=rB+\mathrm{tr}(B)sI$. Accordingly, fix $1\leq i,j\leq n$ with
$i\neq j$, let $K_{ij}$ denote the span of $\{e_i,e_j\}$, and let
$P_{ij}$ denote the orthogonal projection onto $K_{ij}$. Notate the
compressions to $K_{ij}$ by
\begin{gather*}\epsilon_{ii}=P_{ij}E_{ii}\Bigm|_{K_{ij}}=
\begin{bmatrix} a_{ii} & b_{ii} \\
c_{ii} & d_{ii} \end{bmatrix}\qquad
%
\epsilon_{ij}=P_{ij}E_{ij}\Bigm|_{K_{ij}}=
\begin{bmatrix} a_{ij} & b_{ij} \\
c_{ij} & d_{ij} \end{bmatrix}\\
%
\epsilon_{ji}=P_{ij}E_{ji}\Bigm|_{K_{ij}}=
\begin{bmatrix} a_{ji} & b_{ji} \\
c_{ji} & d_{ji} \end{bmatrix}\qquad
%
\epsilon_{jj}=P_{ij}E_{jj}\Bigm|_{K_{ij}}=
\begin{bmatrix} a_{jj} & b_{jj} \\
c_{jj} & d_{jj} \end{bmatrix}
%
\end{gather*}
Let $U_{1}$, $U_{2}$, and $U_{3}$ be the unitary matrices which fix
the orthogonal complement of $K_{ij}$ with action on $K_{ij}$ given by
$u_{1}=\begin{bmatrix} i & 0 \\ 0 & 1\end{bmatrix}$,
$u_{2}=\begin{bmatrix} 0 & 1 \\ 1 & 0\end{bmatrix}$, and
$u_{3}=\frac{1}{\sqrt 2}\begin{bmatrix} 1 & 1 \\
-1 & 1\end{bmatrix}$ respectively. Note that for
$A=\begin{bmatrix} a & b \\ c &
d\end{bmatrix}\in\mathbb{C}\times\mathbb{C},$
$$u_1Au_1^\dagger=\begin{bmatrix} a & -ib \\ ic & d\end{bmatrix}\quad
u_2Au_2^\dagger=\begin{bmatrix} d & c \\ b & a\end{bmatrix}\quad
u_3Au_3^\dagger= \frac12\begin{bmatrix}
a+b+c+d & -a+b-c+d \\ -a-b+c+d & a-b-c+d \end{bmatrix}$$
For $k=1,2,3$, $U_kP_{ij}=P_{ij}U_k$ so that
$\mathcal{U}_{k}\Phi= \Phi\mathcal{U}_{k}$ implies for
$\ell,m\in\{i,j\}$ ,
\begin{equation}
\label{eq:compression}
u_k\epsilon_{\ell m}u_k^\dagger=
P_{ij}\Phi(U_ke_{\ell m}U_k^\dagger)\Bigm|_{K_{ij}}\tag{2}
\end{equation}
With $k=1$, equation (2) shows that the off-diagonal entries of
$\epsilon_{ii}$ and $\epsilon_{jj}$ equal zero and that all entries of
$\epsilon_{ij}$ and $\epsilon_{ji}$ except for $b_{ij}$ and $c_{ji}$
must equal zero. Since $\mathcal{U}_2(e_{ii})=e_{jj}$, $a_{ii}=d_{jj}$
and $d_{ii}=a_{jj}$, and since $\mathcal{U}_2(e_{ij})=e_{ji}$,
$b_{ij}=c_{ji}$. With these identities, it follows that
$2u_3\epsilon_{ij}u_3^\dagger=
\begin{bmatrix}b_{ij} & b_{ij} \\ -b_{ij} &
-b_{ij} \end{bmatrix}$. Further, since
$2U_3e_{ij}U_3^\dagger=e_{ii}+e_{ij}-e_{ii}-e_{ii}$,
$$\begin{bmatrix}b_{ij} & b_{ij} \\ -b_{ij} &
-b_{ij} \end{bmatrix}=
\begin{bmatrix} a_{ii}-d_{ii} & b_{ij} \\
-b_{ij} & d_{ii}-a_{ii}\end{bmatrix},$$
so that $a_{ii}-d_{ii} = b_{ij}$. Letting $r=b_{ij}$ and $s=d_{ii}$ one has,
$$\epsilon_{ii}=\begin{bmatrix} s+r & 0 \\
0 & s \end{bmatrix}\quad
%
\epsilon_{ij}=\begin{bmatrix} 0 & r \\
0 & 0 \end{bmatrix}\qquad
%
\epsilon_{ji}=\begin{bmatrix} 0 & 0 \\
r & 0 \end{bmatrix}\qquad
%
\epsilon_{jj}=\begin{bmatrix} s & 0\\
0 & s+r \end{bmatrix}\qquad
%
$$
Letting $i,j$ run through all unequal pairs yields equation (1).
Remark. Using similar techniques one can show the following.
Let $\Phi\in\mathcal{L}(\mathbb{C}^{n\times n})$ and suppose that for any
orthogonal conjugation operator
$\mathcal{O}\in\mathcal{L}(\mathbb{C}^{n\times n})$,
$\mathcal{O}\Phi=\Phi\mathcal{O}$. There exists $r,s,t\in\mathbb{C}$
such that, for $B\in\mathbb{C}^{n\times n}$,
$\Phi(B)=rB+sB^\top+\mathrm{tr}(B)tI$.