The reason I'm asking is because characteristic\indicator functions have no smooth derivatives and plus I don't understand in which function space of test functions the authors define the weak formulation of (59).

I think we can agree that whatever the authors mean by "$\psi(u_n(t)\chi_{(0,t)}$", they must realise that it isn't a 'test function' in the traditional (smooth, compactly supported) sense, given its insufficient smoothness in both the $x$ and $t$ variables. However, it is possible that 'test function' means something more general here. The sequences $(u_n)$ and $(f_n)$ being considered here *are* test functions in the traditional sense, so we can still plug things into the natural bilinear pairing between $\mathscr{D}(Q) = C_0^\infty(Q)$ and its dual $\mathscr{D}'(Q)$ (the space of distributions on $Q$) and draw conclusions from this.

The fact that $P(u_n) = f_n$ in $\mathscr{D}(Q)$ implies (tautologically) that

$$
\langle \Psi, P(u_n)\rangle = \langle \Psi,f_n\rangle
$$

for all distributions $\Psi\in \mathscr{D}'(Q)$. We can therefore learn things about $u_n$ and $f_n$ by making clever choices of 'test distribution' $\Psi$. When $\Psi$ is integration against a function, said function may arguably be termed a 'test function', which I think may be what is being proposed here.

This sort of terminology occurs in other settings where one has a non-degenerate bilinear pairing; 'test function' often just means 'thing we're free to choose which will fit into the other slot of the pairing from what we already have'.