More details: $E_{ij}$ is the average of $A_{\sigma(i)\sigma(j)}$ over $\sigma$. When $i=j$ you get $E_{ii} = p A_{ii} + \frac{1-p}{n-1}\sum_{k\neq i} A_{kk}$ where $p$ is the probability that $\sigma(i)=i$. To calculate $p$ note that when $\sigma(n)=n$, the restriction of $\sigma$ to $[1,n-1]$ has precisely one less cycle. It follows that $\sum_{\sigma(n)=n}w_n(\sigma) = \frac{\theta}{\theta+n-1} \sum_{\tau\in S_{n-1}} w_{n-1}(\tau)$ where $w_n$ is the weighing above on $S_n$.

Similarly, for $i\neq j$ we have a sum over $A_{kl}$ with different weights. Easy cases include $k=i, l=j$ (this is $\frac{\theta^2}{(\theta+n-1)(\theta+n-2)}$ for the same reason as the diagonal), $k=j, l=i$ (this is $\frac{\theta}{(\theta+n-1)(\theta+n-2)}$ since now there is only one more cycle than the restriction to $[1,n-2]$), $k=i, l\neq i,j$ or $k\neq i,j, l=j$ (this is $\frac{\theta}{\theta+n-1}\left(1-\frac{\theta}{\theta+n-2}\right)$ since we fix one co-ordinate but not the other).

A bit more difficult is the case $k=j, l\neq i,j$ and $i\neq i,j, l=i$. Again for $i=n, j=n-1$ this means that $n, n-1$ are consecutive on a cycle of length at least $3$. The probability that $n$ is on a cycle of length $3$ is the complement of the probability that it is on a cycle of length $1$ or $2$ (that is $\frac{\theta}{\theta+n-1} + (n-1)\frac{\theta}{(\theta+n-1)(\theta+n-2)}$. Given that, the successor to $n$ on the cycle is uniformly distributed, so the probability that $k=j, l\neq i,j$ is $\frac{1}{n-1} \left(1 - \frac{\theta}{\theta+n-1} + (n-1)\frac{\theta}{(\theta+n-1)(\theta+n-2)}\right)$.

Finally, the probability that $k,l$ are distinct from $i,j$ is the complement of the above cases, and if so then the pair $k,l$ is uniformly distributed on the (n-2)(n-3) possibilities. It turns out each of these possibilities has probability $\frac{1}{(\theta+n-1)(\theta+n-2)}$.

2 added a line break.

You can calculate this directly. All you need to know is the probability that $(\sigma(i),\sigma(j))=(k,l)$ for each pair $(k,l)$. Write $e_i$ for the $i$th basis vector, $v=\sum_i e_i$. Let $E = \mathbb{E}(A)$. Then:

$$e_i^T E e_i = \frac{\theta-1}{\theta+n-1} e_i^TAe_i + \frac{1}{\theta+n-1}\mathrm{Tr} A$$ and $$(\theta+n-1)(\theta+n-2)e_i^T \begin{equation*} \begin{split} &(\theta+n-1)(\theta+n-2)e_i^T E e_j = e_j\\ &= v^TAv - \mathrm{Tr} A + (\theta-1) \left[e_i^TAv+v^TAe_j+e_j^TAe_i-e_i^TAe_i+e_j^TAe_j\right] + (\theta-1)^2 e_i^TAe_j$$\end{split} \end{equation*}

1

You can calculate this directly. All you need to know is the probability that $(\sigma(i),\sigma(j))=(k,l)$ for each pair $(k,l)$. Write $e_i$ for the $i$th basis vector, $v=\sum_i e_i$. Let $E = \mathbb{E}(A)$. Then:

$$e_i^T E e_i = \frac{\theta-1}{\theta+n-1} e_i^TAe_i + \frac{1}{\theta+n-1}\mathrm{Tr} A$$ and $$(\theta+n-1)(\theta+n-2)e_i^T E e_j = v^TAv - \mathrm{Tr} A + (\theta-1) \left[e_i^TAv+v^TAe_j+e_j^TAe_i-e_i^TAe_i+e_j^TAe_j\right] + (\theta-1)^2 e_i^TAe_j$$