The Direct Answer.
Here's the short and direct answer. On an oriented $n$-dimensional smooth manifold $M$, we can define the integrals of $n$-forms. So, if $\omega$ is an $n$-form on $M$, then we can define the number $\int_M\omega\in\Bbb{R}$, called the integral of the differential $n$-form over $M$ (ok I'm glossing over some integrablity conditions for the sake of brevity). So, the 'things' we are supposed to integrate on $n$-dimensional oriented smooth manifolds are differential $n$-forms. You write
which is just a tensor field. Why does it make sense to put a $\int$ before a tensor field? I don't think $\int$ is defined to operate on tensor fields.
and yes, $\phi\cdot\epsilon$ is a tensor field, but it's not some random tensor field; it is an alternating $(0,n)$ tensor field (i.e an $n$-form), and this is what allows us to define its integral.
Brief definition of Integrals of $n$-forms.
The $n$-form $\epsilon$ is what a mathematician would call a volume form (and in this case it's the volume form induced by the Lorentzian metric, and it defines the orientation). So, if you have a function $\phi:M\to\Bbb{R}$, then $\phi\cdot\epsilon$ is again an $n$-form, and this is the 'correct' thing to be integrating on $M$ (because it is oriented).
The topic of integration on manifolds is a basic part of differential geometry, and is covered in any good textbook (e.g Spivak, Lee). The definition of the integral $\int_M\omega$ is slightly subtle, and the details are covered in these textbooks. The key thing is that if you have a coordinate chart $(U,\alpha=(x^1,\dots, x^n))$ which is positively oriented, then you can write $\omega=f\,dx^1\dots\wedge x^n$ for some unique function $f:U\to\Bbb{R}$, and we have
\begin{align}
\int_U\omega&:=\int_{\alpha[U]}(f\circ\alpha^{-1})\cdot d\lambda_n\equiv \int_{\alpha[U]}f\circ\alpha^{-1},
\end{align}
where the symbols on the right are the standard multivariable calculus integrals of functions defined on subsets of $\Bbb{R}^n$ (i.e the usual integral with respect to $n$-dimensional Lebesgue measure... or you can just use Riemann integrals throughout). In words, this is saying the most obvious thing you can do: to integrate an $n$-form in a coordinate patch, write the $n$-form as a function times a wedge of $dx^i$'s and then integrate that chart-dependent function over the corresponding subset of $\Bbb{R}^n$. To then define $\int_M\omega$ from these various $\int_U\omega$, one uses what's called a partition of unity, and the purpose of me emphasizing orientation is so that in this global step, we choose the correct signs... i.e we have to put humpty dumpty back together correctly.
More Context, and Reconciling Various Ideas.
Perhaps the following answer of mine will be helpful in some details. In reality, the correct things to integrate on an arbitrary smooth manifold are scalar densities. However, if your manifold is oriented, then you can construct an isomorphism from the vector space of scalar densities to the vector space of $n$-forms. So, one then, by 'transport of structure', define integrals for $n$-forms.
I think what you're most confused by is that we're integrating complicated types of objects (scalar densities/$n$-forms) as opposed to smooth functions. But this shouldn't be surprising at all! When doing integrals, we don't just integrate functions $f$. We integrate functions with respect to some measure $\mu$ (in $\Bbb{R}$, this is the integral with respect to $1$-dimensional Lebesgue measure, which for continuous functions is equal to the basic vanilla Riemann integral we all know and love). In $\Bbb{R}^n$, we integrate functions $f$ with respect to $n$-dimensional Lebesgue measure: $\int_{\Bbb{R}^n}f\,d\lambda_n$. Or more generally on an abstract measure space $(X,\mathfrak{M},\mu)$, we integrate functions $f:X\to\Bbb{R}$ with respect to the measure $\mu$, to get the number $\int_Xf\,d\mu$. The question now becomes: on an abstract manifold $M$, what is the measure? Is there a natural measure? In the general case the answer is no, but for GR, the answer is yes! The metric tensor $g$ gives rise to a unique measure $\mu_g$, also denoted $\lambda_g$ or $dV_g$, which a mathematician might call the Riemann-Lebesgue volume measure on $M$ (see here for details of the definition and construction). So, we can now talk about integrating functions on smooth manifolds $\int_Mf\,dV_g$.
So, the point is that there are many ways to develop the subject of integration on manifolds:
- If you're more analysis-minded, you might go the Lebesgue route and ask 'what measure should I use'? In the case of Riemannian/Lorentzian geometry, you'd end up using the measure $dV_g$. In classical mechanics (particularly the Hamiltonian formulation on phase space, you'd use the measure $dV_{\Omega}$ induced by the symplectic form). Once you have the measure, you can start integrating functions with respect to this measure.
- If you were more differential-geometry minded, you'd say 'hmmm...what are the objects I should be integrating? Well, a manifold comes with a bunch of local coordinate charts, so maybe I can use this local information to somehow put together a global definition?' Pursuing this line of thought leads you to the notion of integrating scalar densities.
- Going one step further, you'd say scalar densities aren't so nice to work with. Also, from single variable calculus, $f(t)\,dt$ can be thought of as a $1$-form, and differential $n$-forms capture the notion of signed-volumes, and differential forms are very easy to compute with. Can I somehow make a definition for integrating $n$-forms? Pursuing this line of thought, you'd see that if you impose the orientedness condition on $M$, then several pesky minus signs disappear and you now have a notion of integration for $n$-forms (particularly you can now integrate the $n$-form $\phi\cdot\epsilon$).
Finally, you can check that regardless of which route you take, all these approaches give you the same answer given correct hypotheses (oriented $n$-dimensional pseudo-Riemannian manifold). So, the fact that you have to end up considering more complicated objects ($n$-forms) is not really an issue, and is in fact a strength in disguise (they're very computationally flexible, and we have the famous Stokes' theorem).