Why is modern control theory in the time domain instead of the S-domain?
I would say that the difference is not about time domain or frequency domain, is more about using transfer function or state space equations.
The state equation for a system with a single input and a single output is
$$\sum_{n=0}^{N} a_n \left(\frac{d}{d t}\right)^n y(t) = \sum_{n=0}^{N} b_n \left(\frac{d}{d t}\right)^n u(t)$$
The state space representation for a system
$$\dot{x}(t) = A \, x(t) + B \, u(t), \\ y(t) = C \, x(t) + D \, u(t)$$
where $y(t)$, $x(t)$ and $u(t)$ all can be vectors, extending this easily to multiple inputs and multiple outputs. For a high order system the number of parameters in state space representation is much larger, and there are multiple possible representations, however there are some canonical representations that are unique to a given system, and may have as few parameters as the transfer function.
Linking the two representations
This can be represented in the Laplace domain as
$$s \cdot X(s) = A \, X(s) + B \, U(s), \\ Y(s) = C \, X(t) + D \, U(s)$$
If you want to understand the dynamic of the system you can manipulate the state update equation as follows
$$s X(s) - A\, X(s) = B \, U(s) \\
 (s \mathbb{I} - A) X(s) = B \, U(s)$$
From the above equation you start to notice that when $X(s)$ is an eigenvector of You have $(s \mathbb{I} - A) v = s v - \lambda v$, and this will define one mode of the system $x(t) = x(0) e^{\lambda t}$, this is why you can study the dynamic of the system in terms of the eigenvalues of $A$.
Proceeding to the transfer function we can solve the above equation for $X(s)$
$X(s) = (s \mathbb{I} - A)^{-1} B \, U(s)$
and replace in the outptut equation
$$Y(s) = C \, X(s) + D \, U(s) \\
C \, (s \mathbb{I} - A)^{-1} B \, U(s) + D \, U(s) \\
(C \, (s \mathbb{I} - A)^{-1} B + D) U(s)
$$
And finally the transfer function is
$$H(s) = \frac{Y(s)}{U(s)} = (C \, (s \mathbb{I} - A)^{-1} B + D) U(s)$$
The result is so difficult to write that I will skip (and that is what is called classical control).
This equation will have rational functions (ratio of polynomials) in $s$, If a system has multiple inputs and multiple outputs it will be a matrix of scalar transfer functions mapping the effect of each input to each output.
Apart from notation compactness, linear algebra methods scale well to higher dimensions. One example of other advantage that I can cite is related to my master's degree. You can study system stability using Lyapunov stability condition, and this leads to linear matrix inequalities that can be applied to uncertain systems (with uncertainty in the model matrices).