Why is the (non-relativistic) stress tensor linear and symmetric?

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.

Concerning the symmetry: if the stress tensor were not symmetric there would be a net torque on the object and it would rotate. For an explanation see the image below. The components $\sigma_{ij}$ represent the shear forces on the $i$th face in the $j$th direction:

Since the stress tensor describes objects in equilibrium, the object doesn't rotate (that is, the net torque is zero -- this is the "conservation of angular momentum part) so $\sigma_{ij}=\sigma_{ji}$.