Edit: Justification for the choice of expansion in weird time derivative.
I think I understand your struggle from your comment. I struggled a lot with this thing, too, because no author (except Chapman himself) bothers to explain where this weird derivative comes from.
The punchline is that, you don't need it. Tong does not use it, Kardar does not use it... Hell, even Landau and Lifshitz do not use it. If you read Kirkpatrick, he uses it; if you read literature on granular systems, they use it. Mostly, academic articles use it because it is convenient, but pedagogical texts don't. They don't use it because this weird derivative is very useful (as a formal tool) but washes off the details of the perturbation. On top of that, you never go practically beyond the first order because the second and third order give you the so-called Burnett and super-Burnett equations that are unstable (no stable homogeneous state). They include the third and fourth derivatives in the space of the hydrodynamical field.
So let me explain (in 1D for notational convenience). The Chapman-Enskog method assumes three things:
The spatial and temporal dependence of the microscopic distribution function can be fully incorporated in the macroscopic fields $f(r, v, t) \equiv f(\rho(r, t), T(r, t), u(r, t)|v)$. The dependence on the macroscopic term is "functional", meaning that you can have a term proportional to $\nabla T$, for example.
$f$ can be perturbatively expanded in gradients of the field: $$f(\rho(r, t), T(r, t), u(r, t)|v)=f^{(0)}+\epsilon f^{(1)} + \epsilon^2 f^{(2)}+\ldots $$
$f^{(0)}$ depends only directly on the fields and is assumed to be the local equilibrium approximation of $f$. $f^{(1)}$ contains term like $\nabla \rho$, $f^{(2)}$ terms like $\nabla\cdot(\nabla T)$. $\epsilon$ is a formal parameter. The validity of the method stands on the smallness of the higher order terms, hence on the smallness of the gradients, and hence on a small Knudsen number. That's why we usually say that $\epsilon$ is related to the Knudsen number (and indeed, you can do an expansion in the Knudsen number). But here, $\epsilon$ is just a term that orders the series in number of gradients.
$\int f^{0} \mathrm{d}v=\int f \mathrm{d}v = \rho$, $\, \int f^{0} v \mathrm{d}v=\int f v \mathrm{d}v = \rho u$ and $\int f^{0} \, 0.5 \, mv^2 \mathrm{d}v=\int f 0.5 \, mv^2 \mathrm{d}v = 0.5\rho T$. This means that the hydrodynamical field is given by the integral of the local equilibrium term.
And that's all. A useful thing one can do (I don't really have a source on that, but I find it really nice) is to rescale the gradient as follows: $\nabla \to \epsilon \nabla$.
Now, let's quickly go over the method. The Boltzmann equation reads:
$$\partial_t f + \epsilon v\partial_x f = J(f, f).$$
If you inject the perturbation series, you are mostly done (after pages of computation):
$$\partial_t (f^{(0)} + \dots) + \epsilon\partial_x (f^{(0)}+\dots) = J(f^{(0)}+\dots, f^{(0)}+\dots).$$
The zeroth order is given by terms without any $\epsilon$. Which you might think is given by:
$$\partial_t f^{(0)} = J(f^{(0)}, f^{(0)}).$$
Please stop reading here and ask yourself if you are really not including any higher order term in $\epsilon$ here. We know that $f^{(n)}$ is of order $\mathcal{O}(\epsilon^0)$ $\forall n$ . But, is $\partial_t f^{(n)}$ of order $\mathcal{O}(\epsilon^0)$? Said differently, can the derivative create terms proportional to $\epsilon^m$?
Of course, I'm asking so you are indeed including higher-order terms. Let's see how. We recall that $f$ has time dependence through the hydrodynamical fields:
$$\partial_t f^{(0)} = (\partial_t \rho) \partial_\rho f^{(0)}+(\partial_t u) \partial_u f^{(0)}+(\partial_t T) \partial_T f^{(0)}.$$
But, the time derivative of a field is just a hydrodynamical equation. Take $\rho$ for example:
$$\partial_t\rho = -\epsilon\partial_x (\rho u).$$
This equation will itself have different gradients, and somehow, if you want your perturbation to be consistent, you also have to correctly keep track of the order of this derivative. Which can, therefore, be expressed as an expansion in terms of gradients:
\begin{split}
\partial_t\rho &= 0 + \epsilon(-\partial_x (\rho u)) + \epsilon^2 0 + \dots\\
& \equiv \partial_t^{(0)}\rho + \epsilon(\partial_t^{(1)}\rho) + \epsilon^2(\partial_t^{(2)}\rho) + \dots
\end{split}
Where in the second line, we defined the symbol $\partial_t^{(i)}$, which you will agree is not a derivative but just the term of a "series" expansion. With this insight, we can go back to our expansion of the derivative of the microscopic velocity distribution:
\begin{split}
\partial_t f^{(0)} &= (\partial_t^{(0)} \rho + \epsilon (\partial_t^{(1)} \rho) + \dots) \partial_\rho f^{(0)}\\
&+(\partial_t^{(0)} u + \epsilon (\partial_t^{(1)} u) + \dots) \partial_u f^{(0)}\\
&+(\partial_t^{(0)} T + \epsilon (\partial_t^{(1)} T) + \dots) \partial_T f^{(0)}.
\end{split}
So, now, let's order the terms in a series of $\epsilon$.
\begin{split}
\partial_t f^{(0)} &= (\partial_t^{(0)} \rho) \partial_\rho f^{(0)}
+(\partial_t^{(0)} u) \partial_u f^{(0)}+(\partial_t^{(0)} T) \partial_T f^{(0)} \\
&+\epsilon\left[(\partial_t^{(1)} \rho) \partial_\rho f^{(0)}
+(\partial_t^{(1)} u) \partial_u f^{(0)}+(\partial_t^{(1)} T) \partial_T f^{(0)}\right] + \dots
\end{split}
More simply, we see that $\partial_t^{(i)}$ is not a "real-time" derivative but still follows a kind of chain rule (and other properties of the regular derivative) such that we can rewrite the above expression as:
$$\partial_t f^{(0)} = \partial_t^{(0)} f^{(0)} + \epsilon\partial_t^{(1)} f^{(0)} + \dots $$
Which is what we wanted to see. Note then that the correct zero-order Boltzmann equation is not:
$$ \partial_t f^{(0)} = J(f^{(0)}, f^{(0)}),$$
but:
$$ \partial_t^{(0)} f^{(0)} = J(f^{(0)}, f^{(0)}).$$
Just to recapitulate:
We can perfectly work with only $\partial_t$. When we do so, we have to be careful because even if $f$ does not contain terms proportional to $ \epsilon$, $\partial_t f$ does contain terms proportional to $\epsilon$. Hence, we must expand the derivative as a series in power of $\epsilon$ to consistently, order by order, arrange all terms in $\epsilon$.
The method above is equivalent to writing $\partial_t g = (\partial_t g)_0 + \epsilon
^1(\partial_{t} g)_1 + \dots$. But a nice thing is that $(\partial_{t} g)_m$ behaves as a derivative. Hence, it suggests we write it down in the form: $\partial_t g = \partial_t^{(0)} g + \epsilon\partial_t^{(1)}g + \dots$. Using the chain rule, it will create terms like $(\partial_t^{(1)}T)\partial_T g$, where $\partial_t^{(1)}T$ has to be understood as the contribution of order $\epsilon$ in the time derivative of $T$ (or equivalently linear in the gradients, note the discrepancy between what I say here and what is said in your article. In your article, $\epsilon^0$ corresponds to a linear gradient; in my case, $\epsilon^0$ corresponds to zero gradient. This is because I rescaled the gradient by $\epsilon$ for simplicity.)
End of edit.
Old post (still relevant):
Yes, I think too that this notation is a bit unfortunate (exactly as for your question on $\partial/\partial_{t_{\text{coll}}}$).
First, the $t_1$, $t_2$, ... notation used in the article is nonsense because you only have one time $t$, and you do not create additional "fake time" as in typical Multiple-scale pertubation theory. What is written in your article as:
$$\partial_t g = \partial_{t_0}g + \epsilon \partial_{t_1}g + \ldots$$
is usually written as:
$$ \partial_t g = \partial_t^{(0)} g + \partial_t^{(1)}g + \dots $$
See, for example, Contemporary Kinetic Theory of Matter (by Kirkpatrick) or The Mathematical Theory of Non-uniform Gases (by Chapman itself, available on the Internet Archive). This notation makes clear the fact that the "derivative symbols are not really derivative (although they act a bit like it). Our master Chapman tells it himself (p. 113 or 137 of the PDF with $\partial_r\equiv \partial_t^{(r)})$:
These time-derivatives are divided into parts, as follows
$$\frac{\partial n}{\partial t} = \sum \frac{\partial_r n}{\partial t}, \quad \frac{\partial \mathbf{c_0}}{\partial t} = \sum \frac{\partial_r \mathbf{c_0}}{\partial t}, \quad \frac{\partial T}{\partial t} = \sum \frac{\partial_r T}{\partial t},$$
where the quantities on the right are not themselves time-derivatives [...].
I would say it's better to write the derivative this way:
$$\partial_t g = (\partial_t g)_0 + \epsilon
^1(\partial_{t} g)_1 + \dots,$$
where we assumed that the result of a time derivative can just be expanded into a series in power of $\epsilon$. Each terms in $(\dots)_i$ is of order $\epsilon^0$. But they will contribute to a different physical process. The zeroth order term will give you the Euler equation (Eq. 15 in your pd); hence, they correspond to phenomena happening on a time scale similar to the inverse wavelength since these are described by waves or advection ($w=ck$). The first-order term will give you the dissipative contributions (viscosity, heat diffusion, ...) to Navier-Stokes, which happen on a diffusive timescale ($w=Dk^2$).
To understand a bit more clearly what it all means. I will just paraphrase your article. Take the Navier-Stokes equation for the momentum:
$${\displaystyle \partial_t(\rho \mathbf {u} )=-\nabla \cdot \left(\rho \mathbf {u} \otimes \mathbf {u} +[p-\zeta (\nabla \cdot \mathbf {u} )]\mathbf {I} -\mu \left[\nabla \mathbf {u} +(\nabla \mathbf {u} )^{\mathrm {T} }-{\tfrac {2}{3}}(\nabla \cdot \mathbf {u} )\mathbf {I} \right]\right).}$$
Then simply:
$${\displaystyle \partial^{(0)}_t(\rho \mathbf {u} )=-\nabla \cdot \left(\rho \mathbf {u} \otimes \mathbf {u} + p \mathbf {I}\right).}$$
$${\displaystyle \partial_t^{(1)}(\rho \mathbf {u} )=\nabla \cdot \left([\zeta (\nabla \cdot \mathbf {u} )]\mathbf {I} +\mu \left[\nabla \mathbf {u} +(\nabla \mathbf {u} )^{\mathrm {T} }-{\tfrac {2}{3}}(\nabla \cdot \mathbf {u} )\mathbf {I} \right]\right).}$$
This somehow justifies the idea that the author had of writing these things in terms of derivative with different "timescales", as in regular multiscale perturbation (but again, that's just completely nonsense mathematically and confuses more than it helps). Indeed, $\epsilon$ keeps track of the gradients. The zeroth order corresponds to the term proportional to a single gradient and the first order to the term proportional to a double gradient (Laplacian) in the transport equations. Hence, $\partial_{t_0}$ corresponds to the change in time of some variable due to the single gradients terms, $\partial_{t_1}$ corresponds to the change in time due to double gradients, and so on.
Concerning your last question, $f^{(0)}$ is the local equilibrium approximation of the probability distribution (which creates the Euler equation and the single gradient expansion), while $f^{(1)}$ contains terms proportional to the gradient of the macroscopic fields (which create the additional terms necessary to obtain the full Navier-Stokes equation).
