# Why is stress linearly dependent on strain in elastic materials?

The displacement vector $\mathbf{u}(\mathbf{x}, t) = \mathbf{r}(\mathbf{x}) - \mathbf{x}$ is used to keep track of the motion of the material points in a material. Firstly, we take the gradient of the displacement vector, yielding the Jacobian, and split the result into symmetric and anti-symmetric components: $$\nabla \mathbf{u} = \mathbf{J} = \frac{1}{2} (\mathbf{J} + \mathbf{J}^\top) +\frac{1}{2} (\mathbf{J} - \mathbf{J}^\top)$$ Secondly, we use the fact that $\frac{1}{2}(\mathbf{J} - \mathbf{J}^\top) \in \mathfrak{so}(3)$ to conclude that the anti-symmetric component represents a rotation. Finally, the strain tensor is defined: $$\varepsilon = \frac{1}{2} (\mathbf{J} + \mathbf{J}^\top)$$

My question is: why do we define strain in that way? Is there some reason that stress is linearly dependent on strain in elastic materials (EM forces)? Why use velocity (Jacobian) instead of position (displacement vector)?

As lemon already pointed out $\mathbf{J}$ is not velocity. It a second-order tensor called the displacement gradient and commonly denoted by $\mathbf{H}$. It has nine components: $\frac{\partial u_i}{\partial x_j}$. Here $x_i$ are coordinates of a material point in the undeformed configuration and $y_i$ are the coordinates of the same point in the current configuration and $u_i = y_i - x_i$.
A good way to motivate $\mathbf{\varepsilon}$ is by looking at the exact kinematics of deformation. One can show that the changes in lengths of curves and changes in angles between curves are characterized by the the right-Cauchy Green deformation tensor $\mathbf{C} = \mathbf{F}^T\mathbf{F}$. Here $\mathbf{F}$ is the deformation gradient with the components $\frac{\partial y_i}{\partial x_j}$. From this, we define the Lagrangian strain tensor $\mathbf{E}$ as $$\mathbf{E} = \frac12\left(\mathbf{C} -\mathbf{I}\right).$$ In terms of the displacement gradient, $$\mathbf{E} = \frac12\left(\mathbf{H}^T + \mathbf{H} + \mathbf{H}^T \mathbf{H}\right) =\frac12\left( \mathbf{\nabla u}^T + \mathbf{\nabla u} + \mathbf{\nabla u}^T \mathbf{\nabla u} \right).$$ Assuming that $\left|\mathbf{\nabla u}\right| << 1$, we can approximate $\mathbf{E}$ by $\mathbf{\varepsilon}$ where $$\mathbf{\varepsilon} = \frac12\left(\mathbf{H}^T + \mathbf{H}\right) = \frac12\left( \mathbf{\nabla u}^T + \mathbf{\nabla u}\right).$$ This captures changes in lengths and angles when the displacement gradient is small. This also captures changes in volumes. And that is the reason, we use this. In many practical applications, this is more than sufficient since the displacement gradients tend to be very small compared to unity.
There is no reason why stress should depend linearly on deformation for elastic materials. However, for many practical applications such as when we deal with the elastic behavior of metallic materials (where the displacment gradients are small compared to unity), we can take stress to be linearly proportional to $\mathbf{\varepsilon}$ assuming initial stress is zero.
And, generally, stress only depends linearly on strain for very small strains (see nonlinear theory). So it's just a semi-empirical approximation. But you can see that this approximation should work in the limit $\epsilon\to 0$ by taking the series expansion of any function of strain...