How is the stress-energy tensor defined?

I'll share my view of the matter hoping that it helps here. First there is the definition of the stress-energy tensor. We define the stress-energy tensor $T_{munu}$ from the matter action $I_M$ as $$T_{munu}=dfrac{-2}{sqrt{-g}}dfrac{delta I_M}{delta g^{munu}}tag{1}.$$

If you want a definition that is the one. Now you may of course ask about motivations behind it. And here I share my personal view, in hopes that it helps you somehow.

Consider first electrodynamics. Maxwell's equations reads $$nabla_mu F^{munu}=j^nutag{2}$$

and they tell us that the electromagnetic current $j^nu$ sources the electromagnetic field.

What definition (1) is really doing is defining $T_{munu}$ as whatever appears in the right-hand side of the Einstein Field Equations. In other words: in the same way that the electromagnetic current can be defined as that which sources electromagnetic fields the stress-energy tensor $T_{munu}$ can be defined as that which sources the gravitational field, because in the end of the day that is the proposal of GR.

The stress-energy tensor comes from thinking of matter as a fluid - that is, even if you choose to think of matter as a bunch of particles, the stress-energy tensor at a given point in spacetime comes from averaging over all the particles in the vicinity of that point.

First consider the Newtonian version. A single particle only has its momentum, $mvec v$, with components $mv_i$, which is a vector, not a tensor. But now consider a bunch of particles in a small neighborhood. Sure, we could take the sum/average of their momenta, and get another vector, but it wouldn't capture all information about their average motion. In particular, they could have zero average momentum, but that would just mean their momenta cancel out, yet they could be moving really fast in all different directions, ie. have high pressure.

But it turns out we can catch that extra information in the quantity $mv_i v_j$ instead of $mv_i$. For each particle, this is a symmetric tensor product, so when we sum it over all the particles, it's still a symmetric tensor. And when two particles have opposite $x$-momenta, say, its $xx$ component doesn't cancel because $v_x$ gets squared by the product - that's how the $ii$ components give us the pressures.

Then, we divide by the volume under consideration to obtain a density. And with the factor of velocity, this becomes flux:

$frac{dm v_i v_j}{dV} = dm v_i frac{dx_j}{dt dx dy dz} = frac{dm v_i}{dA dt}$

Since the $dx_j$ is the one that canceled, $dA$ is the product of the other two, ie. this is the flux of $i$-momentum through a surface of constant $x_j$. (Note that $dA$ has a sign, according to whether $v_j$ is in the negative or positive $j$-direction.)

Since transfer of momentum is equivalent to force, this captures all forces between the fluid element and the outside world. $T_{ii}$ are the directional pressures, ie. how the element is trying to expand, and $T_{ij}$ are the shear stresses - how it is trying to rotate via friction. In either case, it makes sense that two particles with opposite velocities have the same effect, hence the quadratic form of this quantity.

In relativity, just replace velocity by 4-velocity, and 3-dimensional volume by 4-dimensional hypervolume. Essentially:

$T_{munu} = Sigma{frac{ dm V_{mu} V_{nu}}{domega}} = Sigma{frac{ dm V_{mu} dx_{nu}}{domega dtau}} = Sigma{frac{ dm V_{mu}}{dV dtau}}$

Where the sum is over all particles in your element, $domega$ is hypervolume, $dtau$ is proper time. We can still call this a "momentum flux density", but now it is flux of 4-momentum through a 3D hypersurface $dV = domega / dx_{nu}$ (hypersurface of constant $x_{nu}$) with respect to proper time. Then you get the extra row and column for the time direction, which do correspond to the average 4-momentum - and of course the time component of 4-momentum is energy, hence $T_{tt}$ is energy density.

So in some sense, we're capturing all the macroscopic behavior of the matter around that point, in terms of how anyone could actually interact with it. Ie. you can only sense the fluid by its pressure, shear stress, and energy/momentum. So perhaps the profound statement of special relativity is that the bulk properties of matter take on this simple quadratic form at each point. And the profound mathematical discovery expressed by the field equations is that, in any Riemannian manifold, there is a curvature tensor (the Einstein tensor, $R_{munu} - frac12 R g_{munu}$) with that same form, and whose variation from point to point behaves exactly like matter does, ie. it has zero divergence, ie. it is "conserved". That's why, in retrospect, Riemannian spacetime seems so inevitable: not only does it unify all these concepts, but it has to contain matter as we know it, since the Einstein tensor must be conserved. It's just a question of how that matter flows, and what do you know, it flows in pretty much the simplest way imaginable, along geodesics!

To appreciate the general definition, we need a general perspective.

The action principle is one of the key principles in classical physics. Loosely speaking, it says that if A influences B, then B must also influence A: all influences must go both ways. The general definition of the stress-energy tensor described below relies on appreciating a precise version of action principle, so I'll start by reviewing the intuition behind the action principle. This definition of the stress-energy tensor is abstract (that's a synonym for general), but it's also intuitively efficient, because an appreciation for the action principle is also needed for several other general insights in classical physics anyway.

The precise version of the action principle can be introduced using a familiar special case, namely a system of $N$ pointlike objects exerting Newtonian forces on each other. Such a system is described by the equations of motion $$ newcommand{bfx}{mathbf{x}} newcommand{bfF}{mathbf{F}} m_nddotbfx_n = bfF_n(bfx_1,bfx_2,...,bfx_N) tag{1} $$ where $bfx_n$ is the location of the $n$th object in space and $m_n$ is its mass. Each overhead dot represents a time-derivative, so equation (1) is just $F=ma$. The force $bfF_n$ on the $n$th object depends on the locations of all the other objects at that time. As it stands, the system of equations (1) does not necessarily satisfy the action principle. To enforce the action principle, we require that the forces $bfF_n$ can all be obtained from a single function $V(bfx_1,...,bfx_N)$ like this: $$ bfF_n = -frac{partial}{partialbfx_n}V. tag{2} $$ This ensures that the forces are all related to each other in a particular way, and by checking a few easy examples, you can convince yourself that it really does enforce the idea that all influences (forces in this case) should go both ways.

If the system does satisfy the action principle (2), then we can derive all of the equations of motion (1) from a single quantity $S$ called the action. The idea is that physically allowed behaviors $bfx_n(t)$ are the behaviors for which $S$ is stationary: $S$ doesn't change to first order if we vary the behavior slightly. This leads to the Euler-Lagrange equations, which in this case are just the equations of motion (1).

General relativity and classical electrodynamics also satisfy the action principle, and in those cases the dynamic entities include fields as well as material objects. In every case, the whole system of equations of motion comes from the requirement that $S$ is stationary under small variations of any physically allowed behavior. In this way, $S$ tells us which behaviors are physically allowed. The fact that all of the equations of motion come from a single action $S$ enforces the idea that all influences go both ways: If the equation of motion for A is affected by B, then the equation of motion for B is affected by A. The stationary-action principle makes this precise.

That general perspective on classical physics is the foundation for a general understanding of conservation laws and for a general definition of conserved currents, including the stress-energy tensor. Here's the idea: Suppose that the action is the sum of two terms, $S=S_f+S_m$. The term $f$ involves only the field, and the term $S_m$ also involves matter. Suppose also that both of these terms are individually invariant under some set of local symmetries — mathematical symmetries that can act differently at each point in space and time. Under these conditions, after making them more precise, we can show very generally that the variation of the matter term $S_m$ with respect to the field is a conserved current. Examples:

• In the case of electrodynamics, the field is the gauge field $A_mu$. The terms $S_f$ and $S_m$ are separately invariant under gauge transformations, and the current $J^mu$ is defined by $delta S_m = int d^4x J^mu delta A_mu$. The current defined by this relation corresponds to the usual electric charge/current density. The current is conserved in the sense that it satisfies $partial_mu J^mu=0$ whenever the matter (and fields) satisfy their equations of motion. In other words, all physically allowed behaviors satisfy the condition $partial_mu J^mu=0$.

• In the case of general relativity, the field is the metric field $g_{munu}$, the symmetry is typically called "diffeomorphism symmetry" (not a great name, but as good as any other name), and the corresponding current $T^{munu}$ — the stress-energy tensor — is defined by $delta S_m propto int d^4x sqrt{|g|}, T^{munu},delta g_{munu}$. The factor of $sqrt{|g|}$ is included so that the measure $d^4x sqrt{|g|}$ is coordinate-invariant. The proportionality constant is specified in the answer by user1620696. The stress-energy tensor is conserved in the sense that $nabla_mu T^{munu}=0$ for all physically allowed behaviors, where $nabla_mu$ is the metric-compatible covariant derivative. This is a consequence of the equations of motion for the matter (which are derived from the same action), but we can also derive it more directly just from the symmetry of $S_m$.

With some effort, we can relate this definition to special cases like the one described in the answer by Adam Herbst. Yes, it does take some effort, but that's usually the price we pay for adopting a general definition: the more general the principles are, the longer the road to practical applications tends to be. To see the whole Grand Canyon, you need to stand a good distance away from it.

Two caveats:

• If the collection of dynamic entities includes spinor fields, then we need to express the metric field $g_{munu}$ in terms of a frame field $e_mu^a$, and the stress-energy tensor is then defined by varying $S_m$ with respect to the frame field.

• Since the definition of $T^{munu}$ only involves $S_m$ and not $S_f$, we can use it even in a model where the field (like the metric field) is just a prescribed background field. But in order for the definition of $T^{munu}$ to be unambiguous, we still need to specify how the theory depends on arbitrary deformations of the background field, because we need that information in order to define the variation with respect to the field. For example, even if we only want to know what the stress-energy tensor is in flat spacetime, we still need to specify how the model would be generalized to slightly curved spacetimes. There are different ways to do that, and they lead to different stress-energy tensors. Thinking through the significance of this is left as an exercise for the reader. Part of my own stumbling journey through that exercise is shamelessly documented in another recent question.