I'll drop the bracket notation for mathematical clarity.
There are an important caveats for the adiabatic theorem. First of all, $H$ is assumed to be non-degenerate, and must stay non-degenerate, ie level crossing are forbidden. Intuitively, you need to match the spectra of $H(t)$ at different times. This is a first sufficient condition.
Also, in general, even when you impose normalisation, the $\phi_n(t)$ are not uniquely specified. There is a gauge invariance which consists of multiplying each of them by a time dependent phase $e^{i\theta_n(t)}$. This is why only closed loops in parameter space give rise to the geometric object known as Berry curvature. This means that when constructing such a basis, you need to constantly resolve this ambiguity. It also means that when you find a candidate, you can construct many more assuming the $\theta_n$ are differentiable.
To make the discussion more simple, I'll assume a finite dimensional Hilbert space, the general case is more technical. One way to construct a solution, only assuming the differentiability of $H$ and the preserved non-degeneracy. The simplest way is to show that they are the solution of an ordinary differential equation (ODE) that you can easily obtain using non degenerate perturbation theory.
Assume you've found such a solution. First of all, you need to define the $E_n(t)$, and check their differentiability. 976
This is done by using Cayley-Hamilton and seeing the spectrum as roots of varying polynomial. You can now use a general property of polynomials depending differentiably on a parameter which says that as long as its roots stay non degenerate, the roots can be matched at all values of the parameter and are differentiable under the parameter. You can easily convince yourself of this fact for low degrees since there are some explicit formulas.
Now that you have the spectrum, say you found a candidate normalised $\phi_n(t)$. Taking the derivative of their defining property $H\phi_n= E_n\phi_n$, you get: $$ (H-E_n)\dot\phi_n= -(\dot H -\dot E_n)\phi_n $$ and you want to invert the $H-E_n$ to get the ODE, but you can't since $E_n$ is an eigenvalue. The equation actually defines $\dot\phi_n$ up to added multiple of $\phi_n$ which is the consequence of the previously mentioned gauge invariance. You can choose to resolve the ambiguity by enforcing $\langle \dot \phi_n,\phi_n\rangle=0$ which has the advantage of being compatible with the normalisation conservation $\Re \langle \dot \phi_n,\phi_n\rangle=0$, and is typically the implicitly chosen convention in perturbation theory. I now defined $n$ perfectly valid ODE's whose solutions satisfy you problem by construction.
Hope this helps.
Edit
I'll just clarify the mathematical reasoning. The French call it a demonstration by "analyse-synthèse" (couldn't find a translation), ie you assume you have the solution and accumulate enough necessary conditions until you are all set to construct the solution from scratch.
For your question, I assumed I found a basis depending differentiably with $t$, and followed through the equations to find a ODE that it would satisfy. The only tricky part is due to the non-unique nature, you need to additionally make a consistent choice so that the ODE is well defined. This is the "analysis." The "synthesis" which I left out is simply checking that the ODE I found actually has a solution. This can easily be seen using classic mathematical theorems such as Cauchy-Lipschitz. In addition, using the unicity of the ODE's solution once specifying the initial conditions (here, assume you start out in a eigenbasis), it is guaranteed that the constructed solution is a solution of the original problem. 976