Here is another way to think about it. It will help to make an analogy with the continuity equation for fluid flow. Suppose you have a fluid which is flowing with some non-uniform velocity $\vec{u}(\vec{r})$. Consider region $\Omega$ of the fluid bounded by a surface $S$, and assume you are told that the total mass in the region $\Omega$ is constant, and moreover the mass density $\rho(\vec{r})$ at each point in $\Omega$ is constant. What could you conclude from this? Well you know there is a mass current in the fluid given by $\vec{J}=\rho \vec{u}$, and since each piece of mass moves continuously, the change in density at a point is equal to the amount of mass flowing into that point. This is expressed mathematically by the "continuity equation" $\dot{\rho} = -\nabla \cdot \vec{J}$. Integrating this over the whole volume, we find $\dot{M} = \int_\Omega \dot{\rho} = \int_\Omega - \nabla \cdot \vec{J} = - \oint_S \vec{J} \cdot \hat{n}\, dA$. On the far right hand side the term $\vec{J} \cdot \hat{n}$ represents the mass flux through the surface.
The stress tensor is in fact not much harder to understand than this. You have two questions: Why is the force on a surface linear in $\hat{n}$? and Why is $\sigma$ symmetric? I will answer these one at a time. Both answers will be done by analogy with the fluid flow.
Why is the force on a surface linear in $\hat{n}$?
Assume we have a large piece of material, and we are told that the momentum density, which I will denote $\vec{p}$, and the angular momentum density, which I will denote $\vec{\ell}$, is constant in some region $\Omega$ with boundary $S$. Let us pick a coordinate system and pick a component of the momentum to look at, say the $i$th component, $p_i$. Then $p_i$ is analogous to $\rho$. Since each little piece of material only exerts forces on its neighbors (there are no long range forces), the $p_i$ must move continuously through the material. Thus the flow of $p_i$ is described by some current $\sigma_{ij}$ which is analogous to $-J_j$ (notice there is a sign convention). Identifying $\dot{p}_i$ with $f_i$, the force per unit volume, the continuity equation gives us $f_i = \partial_j \sigma_{ij}$ (remembering the sign convention). Integrating over the region $\Omega$, with $P_i$ being $i$th component of the total momentum, we find $\dot{P}_i = \int_\Omega \dot{p}_i = \int_\Omega \partial_j \sigma_{ij} = \oint_S \sigma_{ij} n_j \, dA$. Thus the term $\sigma_{ij} n_j$ has the interpretation of momentum flux through the surface, or in other words, the force per unit area on the surface. Thus you the answer to the first part of your question is that the force per unit area is linear in $\hat{n}$ for the same reason that mass flux through an area is linear in $\hat{n}$, and this reason is that there is some current describing how the mass (or momentum) is moving through the medium, and the flux is just the current dotted into $\hat{n}$.
Why is $\sigma$ symmetric?
Now lets address the second part of your question, why must $\sigma_{ij}$ be symmetric if the object is in equilibrium. Let us now consider the $i$th component of the angular momentum density $\ell_i$. We know that an external agent exerting a force per unit area $\vec{f}$ at a point $\vec{r}$ on the surface is exerting a torque whose $i$th component is given by $\tau_i = \epsilon_{ijk} r_j f_k$. However we know from the previous paragraph that $f_k = \sigma_{kh} n_h$. Thus we conclude that $\tau_i$, which is the flux of $\ell_i$, is given by $\epsilon_{ijk}r_j \sigma_{kh} n_h$. We know that this ought to be a current dotted into $\hat{n}$, so the $\ell_i$ current must be $\epsilon_{ijk}r_j \sigma_{kh}$. The change in $\ell_i$ per unit volume, which is $i$th component of the torque per unit volume $\tau_i$, is the divergence of this current: $\dot{\ell}_i = \partial_h \epsilon_{ijk}r_j \sigma_{kh} = \epsilon_{ijk} (\partial_h r_j) \sigma_{kh} + \epsilon_{ijk}r_j (\partial_h \sigma_{kh}) = \epsilon_{ijk} \delta_{hj} \sigma_{kh} + \epsilon_{ijk}r_j f_k = \epsilon_{ijk} \sigma_{kj} + \epsilon_{ijk}r_j f_k.$
The second term is $\vec{r} \times \vec{f}$ as expected, this takes into account the angular momentum produced by uniform translation of the small piece of material. The other term is the antisymmetric part of $\sigma$ and it represents a rotation of the small piece of material about its center of mass. To show that $\sigma$ must be symmetric at an arbitrary point $\vec{r}$ we first move the origin to $\vec{r}$ and then find the expression for $\dot{\ell}_i$, which must be zreo since the object is in equillibrium. We find $0 = \dot{\ell}_i = \epsilon_{ijk} \sigma_{kj}$ where the $\vec{r} \times \vec{f}$ term was dropped because $\vec{r}$ is zero. Thus we find that the antisymmetric part of $\sigma$ must be 0.
You're definitely on the right track but the relationship that you're looking for will depend upon which way you want to model your matter. Dust and radiation are the two models that work the best and are almost equivalent in the answers they produce. Of course the general definition of $ T_{\mu\nu} $ is “the flux of four-momentum $p_{\mu}$ across a surface of constant $ x_{\nu}$”. Since a dust cloud is a collection of moving particles with a somewhat fixed velocity in the inertial reference frame, we can figure out the ( particle ) flux with the velocity and the number of particles by defining the four velocity, $U_{\nu}$, and the number density of the particles $ n $. By this, we can determine the ( four ) flux as $ N_{\nu} = nU_{\nu} $.
Don't forget since our $0^{th}$ index is always full of fun, the $ N_{0} $ component is the number density of the particles, and the non-zero indices correspond to flux in the direction of whichever index you're dealing with.
We're almost there I promise! From that, define your energy density in terms of what we have which is just the particle number density, $n$, and the particle mass to give our $ T_{00} $ term energy density term: $\rho = nm $. Typically, it would be $ nmc^{2} $, but we take units such that $ c = 1 $. The zero index term of the four momentum using these units comes to $ \frac {E} {c^{2}} $ which you should notice to be $ m $.
We now have all the constituents to create the entire stress-energy tensor, we have a component for it, and the relation between our four momentum and flux that brought us that first component. The rest can then be generalized as follows: $$ T_{\mu\nu} = p_{\mu} N_{\nu}=nmU_{\mu}U_{\nu}=\rho U_{\mu}U_{\nu} $$
Where we have either the tensor product between our momentum and flux ( for particles, we can make one massive particle as well and it be just as simple ) vectors, or the tenor product of our velocities with a scale factor $\rho = nm$.
Taking the mixed stress-energy tensor and $ T^{i}_{j} = \nabla_{j}p^{i} $, the connection $ \nabla_{j}$ acting on a vector $ p^{i}$, reduces to $ \frac {\partial p^{i}}{\partial x_{j}}$. As far as I know, that doesn't have any real physical significance.
All in all though, the equation above is the one that you're most likely going to want to go with. Introducing the $\nabla_{j}$ pokes at taking covariant derivatives of tensors, which requires an extra term in the answer from the regular partial differentiation called ( you guessed it ) the affine connection, but when you make $ \nabla_{j}T^{ik} = 0 $ or equivalently $ T^{ik}_{\:\:\:;j} = 0 $ we have conservation of energy and momentum!
All of this was referenced directly from these lecture notes(particularly lecture 1), if I come off as confusing or unclear. The author does a far better job than I.
Best Answer
Here we shall only discuss general relativistic diffeomorphism-invariant matter theories in a curved spacetime in the classical limit $\hbar\to 0$ for simplicity. In particular, we will not discuss the SEM pseudotensor for the gravitational field, but only the stress-energy-momentum (SEM) tensor for matter (m) fields $\Phi^A$. We emphasize that our conclusions will be independent of whether the EFEs are satisfied or not.
I) On one hand, the basic Hilbert/metric SEM tensor-density$^1$ is manifestly symmetric$^2$
$$\tag{1} \sqrt{|g|}T^{\mu\nu}~\equiv~{\cal T}^{\mu\nu}~:=~-2\frac{\delta S_m}{\delta g_{\mu\nu}},$$
cf. a comment by Lubos Motl. However, note that the basic definition (1) is not applicable to e.g. fermionic matter in a curved spacetime, cf. Section II.
Diffeomorphism invariance leads (via Noether's 2nd theorem) to an off-shell identity. Using the matter eqs. of motion (eom)
$$\tag{2} \frac{\delta S_m}{\delta \Phi^A}~\stackrel{m}{\approx}~0, $$
the corresponding Noether's 2nd identity reads$^3$
$$\tag{3} \nabla_{\mu} T^{\mu\nu}~\stackrel{m}{\approx}~0, $$
cf. e.g. Ref. 1. [Here the $\stackrel{m}{\approx}$ symbol means equality modulo matter eoms. The connection $\nabla$ is the Levi-Civita connection.] Eq. (3) serves as an important consistency check. A matter source $T^{\mu\nu}$ to the EFEs should satisfy eq. (3), cf. the (differential) Bianchi identity.
II) In the Cartan formalism, the fundamental gravitational field is not the metric tensor $g_{\mu\nu}$ but instead a vielbein $e^a{}_{\mu}$. The generalized Hilbert SEM tensor-density is defined as
$$\tag{4}{\cal T}^{\mu}{}_{\nu}~:=~{\cal T}^{\mu}{}_a~e^a{}_{\nu}, \qquad {\cal T}^{\mu}{}_a ~:=~- \frac{\delta S_m}{\delta e^a{}_{\mu}}, $$
which is no longer manifestly symmetric, cf. e.g. Ref. 2.
Next we have two symmetries: local Lorentz symmetry and diffeomorphism invariance.
Firstly, local Lorentz symmetry leads (via Noether's 2nd theorem) to an off-shell identity. Using the matter eoms (2), the corresponding Noether's 2nd identity reads
$$\tag{5} {\cal T}^{\mu\nu}~\stackrel{m}{\approx}~{\cal T}^{\nu\mu},$$
i.e the generalized Hilbert SEM tensor-density (4) is still symmetric when the matter eoms are satisfied.
Secondly, diffeomorphism invariance leads (via Noether's 2nd theorem) to an off-shell identity$^4$
$$\tag{6} d_{\mu}{\cal T}^{\mu}{}_{\nu} ~=~{\cal T}^{\mu}{}_a ~d_{\nu}e^a{}_{\mu} -\frac{\delta S_m}{\delta \Phi^A}~d_{\nu}\Phi^A. $$
Not surprisingly, eqs. (5), (6), and $(\nabla_{\nu} e)^a{}_{\mu}=0$ imply eq. (3).
III) On the other hand, the canonical SEM tensor-density
$$\tag{7} \Theta^{\mu}{}_{\nu} ~:=~\frac{\partial {\cal L}_m}{\partial\Phi^A_{,\mu}}\Phi^A_{,\nu} -\delta^{\mu}_{\nu}{\cal L}_m$$
is not always symmetric, cf. e.g. this Phys.SE post. The fact that the Lagrangian density ${\cal L}_m$ has no explicit spacetime dependence leads (via Noether's 1st theorem) to an off-shell identity
$$\tag{8} d_{\mu}\Theta^{\mu}{}_{\nu} ~=~-\frac{\delta S_m}{\delta e^a{}_{\mu}}~d_{\nu}e^a{}_{\mu} -\frac{\delta S_m}{\delta \Phi^A}~d_{\nu}\Phi^A. $$
Recall that the gravitational field, the vielbein $e^a{}_{\mu}$, is not necessarily on-shell. Remember we are doing FT in curved spacetime rather than GR. As a consequence, the first term on the right-hand side of Noether's 1st identity (8) breaks the usual interpretation of Noether's 1st theorem as leading to an on-shell conservation law. [It is comforting to see that it gets restored for a constant vielbein with $d_{\nu}e^a{}_{\mu}=0$.]
IV) Eqs. (4), (6), and (8) imply that
$$\tag{9} d_{\mu}({\cal T}^{\mu}{}_{\nu}-\Theta^{\mu}{}_{\nu})~=~0.$$
One may show that there in general exists a Belinfante improvement tensor-density
$$\tag{10} {\cal B}^{\lambda\mu}{}_{\nu}~=~-{\cal B}^{\mu\lambda}{}_{\nu},$$
such that
$$\tag{11} {\cal T}^{\mu}{}_{\nu}-\Theta^{\mu}{}_{\nu} ~=~d_{\lambda}{\cal B}^{\lambda\mu}{}_{\nu},$$
cf. e.g. my Phys.SE answer here. Note that eqs. (10) and (11) are consistent with eq. (9).
V) Eq. (11) serves as an important consistency check of the Hilbert SEM tensor-density (4) vs. the canonical SEM tensor-density (7). Eq. (11) implies that the two corresponding Noether charges, the energy-momentum $4$-covectors
$$\tag{12} P_{\nu} ~:=~ \int\! d^3x~{\cal T}^0{}_{\nu} \quad\text{and}\quad \Pi_{\nu} ~:=~ \int\! d^3x~\Theta^0{}_{\nu} $$
are equal up to spatial boundary terms
$$\tag{13} P_{\nu}-\Pi_{\nu}~=~\int\! d^3x~d_i{\cal B}^{i0}{}_{\nu}~\sim~0,$$
cf. the divergence theorem.
References:
R.M. Wald, GR; Appendix E.1.
D.Z. Freedman & A. Van Proeyen, SUGRA, 2012; p. 181.
--
$^1$ A tensor-density ${\cal T}^{\mu\nu}= e T^{\mu\nu}$ is in this context just a tensor $T^{\mu\nu}$ multiplied with the density $e=\sqrt{|g|}$.
$^2$ Conventions: In this answer, we will use $(+,-,-,-)$ Minkowski sign convention. Greek indices $\mu,\nu,\lambda, \ldots,$ are so-called curved indices, while Roman indices $a,b,c, \ldots,$ are so-called flat indices.
$^3$ Note that eq. (3) is not a conservation law by itself. To get a conservation law, we need a Killing vector field, cf. e.g. my Phys.SE answer here.
$^4$ Here we have assumed that the matter fields $\Phi^A$ only carry flat or spinorial indices, cf. the setting of my Phys.SE answer here. If $\Phi^A$ also have curved indices, there will be further terms in eq. (6) proportional to the matter eoms.