Physical picture
The situation you are describing are gravitational waves (small perturbations to the metric) propagating on either a flat space metric $\eta_{\mu\nu}$ or some curved geometry $g^B_{\mu\nu}$. These waves do carry momentum and energy, and so there is a contribution to the full stress tensor $T_{\mu\nu}$ from these waves. However, since these waves are small, we can treat the contribution of the stress-energy tensor of the waves on the geometry (in the lingo, this is referred to as backreaction of the GWs on the metric) in a perturbative fashion. In particular, if we are working to leading order (purely linear in the metric perturbation), it is self-consistent to ignore the effect of the stress-energy tensor of the GWs, on the geometry.
Mathematical details
Recasting the exact equations
The exact Einstein's equations (in units with $c=1$) are
\begin{equation}
G_{\mu\nu} = 8\pi G_N T_{\mu\nu}
\end{equation}
where $G_N$ is Newton's constant.
Now let's consider just the left hand side. Let's write the metric as flat space $\eta_{\mu\nu}$ plus a small perturbation $h_{\mu\nu}$
\begin{equation}
g_{\mu\nu} = \eta_{\mu\nu} + \sqrt{G_N} h_{\mu\nu}
\end{equation}
The reason for putting the $\sqrt{G_N}$ will become apparent later, but keep in mind that $G_N$ is a small number and we will use it as our expansion parameter. By the way can go through all the subsequent steps with any background metric $g^B_{\mu\nu}$ instead of $\eta_{\mu\nu}$; the content is the same, but you will need to keep track of some extra details like using the covariant derivative associated with the background metric, and defining the background stress-energy tensor.
Then we can write
\begin{equation}
G_{\mu\nu} = \sqrt{G}_N \hat{\mathcal{E}}^{(1)} h_{\mu\nu} + G_N \hat{\mathcal{E}}[h]^{(2)} h_{\mu\nu} + \cdots + G_N^{N/2} \hat{\mathcal{E}}[h]^{(N)} h_{\mu\nu}
\end{equation}
where $\hat{\mathcal{E}}[h]^{(N)}$ is a differential operator containing two derivatives, and $h$ raised to the power $N-1$. Schematically, $\hat{\mathcal{E}}[h]^{(2)} h \sim h^2 \partial^2 h + h \partial h \partial h$. Note that the leading order operator, $\mathcal{E}^{(1)}$, does not depend on $h$ at all.
In Minkowski spacetime, $T_{\mu\nu}=0$. So Einstein's equations can be written as
\begin{equation}
\sqrt{G}_N \hat{\mathcal{E}}^{(1)} h_{\mu\nu} + G_N \hat{\mathcal{E}}[h]^{(2)} h_{\mu\nu} + \cdots + G_N^{N/2} \hat{\mathcal{E}}[h]^{(N)} h_{\mu\nu}= 0
\end{equation}
Note that up to this point, we actually have not made any approximations. We have written the exact Einstein equations with $T_{\mu\nu}=0$, but just in some weird choice of variables, and without explicitly computing the operators $\hat{\mathcal{E}}^{(N)}$. (Of course you can find these in the literature if you want, although my notation may not be standard).
Iterative solution
Now we develop an iterative solution for the metric perturbation $h_{\mu\nu}$. We write
\begin{equation}
h_{\mu\nu} = h^{(1)}_{\mu\nu} + \sqrt{G_N} h^{(2)}_{\mu\nu} + ... + G_N^{N/2} h^{(N)}_{\mu\nu}
\end{equation}
We then follow the following iterative procedure:
- Plug $h^{(1)}_{\mu\nu}$ into Einstein's equations and work to leading order in the small parameter $G_N$.
- Plug the solution $h^{(1)}_{\mu\nu}$ into Einstein's equations and solve the result for $h^{(2)}_{\mu\nu}$.
- Repeat as many times as desired.
Leading order
At leading order $\mathcal{O}(G_N^{1/2})$, we find
\begin{equation}
\hat{\mathcal{E}^{(1)}} h^{(1)}_{\mu\nu} = 0
\end{equation}
In de Donder gauge, this equation becomes
\begin{equation}
\square h^{(1)}_{\mu\nu} = 0
\end{equation}
This is a standard wave equation, and has standard wave solutions. If you are careful about keeping track of the gauge conditions and constraints, you will end up finding there are two independent solutions, corresponding to the $+$ and $\times$ polarizations.
Next to leading order
Now that we have solved this equation, we work to the next order $\mathcal{O}(G_N)$ for $h^{(2)}_{\mu\nu}$.
\begin{equation}
\hat{\mathcal{E}^{(1)}} h^{(2)}_{\mu\nu} + \hat{\mathcal{E}^{(2)}}[h^{(1)}]h^{(1)}_{\mu\nu} = 0
\end{equation}
We can rearrange this equation as
\begin{equation}
\hat{\mathcal{E}^{(1)}} h^{(2)}_{\mu\nu} = T^{(\rm eff)}[h^{(1)}]_{\mu\nu}
\end{equation}
where $T^{\rm eff}[h^{(1)}]_{\mu\nu}$ is an effective stress-energy tensor that depends on the first order perturbation $h^{(1)}_{\mu\nu}$. Therefore, you can think of the stress-energy associated with gravitational waves as sourcing a higher order correction to the metric itself. This leads to some very interesting phenomena, such as non-linear gravitational memory.
There is an important caveat about thinking of $T^{\rm eff}[h^{(1)}]_{\mu\nu}$ as a real stress energy tensor: it is not gauge invariant. In particular, the stress energy tensor for gravitational waves depends on the choice of coordinates. Therefore it is not truly a stress-energy tensor (it is just some part of the Einstein tensor we have moved to the right hand side). Often, the physical insight you get from thinking of this term as if it were stress-energy tensor outweighs this mathematical fact, but if you aren't careful this can lead to confusion.
Having said that, as @AndrewSteane pointed out in the comments, you can obtain a gauge-invariant non-local stress-energy tensor by averaging the effective stress-energy tensor over a region consisting of several wavelengths.
What's generating the waves?
Your final question was: how can there be waves if nothing has generated the waves?
In fact this situation is not special to gravity. Mathematically the point is that solutions to the equations of motion are only unique if you specify asymptotic boundary conditions. If you specify that the metric must be asymptotically flat, then these boundary conditions will rule out plane wave solutions, and the only way to get GWs is to have a source somewhere in the spacetime (in which case $T_{\mu\nu}$ is not zero). However, mathematically it is possible to consider more general boundary conditions where plane waves can come in from infinity and go out to infinity. Physically, you can think of this as an idealization where an observer lives in some local frame and experiences a passing gravitational wave from very far away. This exact same situation can happen in any theory described by a wave equation, such as electromagnetism.
In more detail, you may notice that the curvature tensor $R_{\mu\nu\rho\sigma}$ has 20 independent components in 4 dimensions, but Einstein's equations $G_{\mu\nu}=8\pi G_{N} T_{\mu\nu}$ are only 10 equations relating components of the curvature to the stress-energy tensor. Anyway, the point is that there are components of the curvature which are not fixed by the matter distribution. These components can be written in terms of the Weyl tensor. The fact that the Weyl tensor is not fixed by the matter distribution, is what mathematically allows for gravitational waves to propagate through the spacetime even when $T_{\mu\nu}=0$.
