Time Dependent Perturbation Theory: How do we know the system stays in the same Hilbert space? In introductory time dependent perturbation theory, I have read that given a system with a Hamiltonian of the form:
$H = H^{(0)} + \delta H(t) \tag{0} $
where the $H^{(0)}$ is independent of time and has known eigenstates given by $|n \rangle $.
Now, if the perturbation is 'small' then the wave function at some later time (after the perturbation has been switched on) can be given as:
$| \Psi(t) \rangle  = \Sigma_{n}c_n(t)e^{{-iH^{(0)}t}/\hbar} |n \rangle     \tag{1}       $
where now the coefficients can be expressed in a perturbative series of the form:
$c_n(t) = c_n^{(0)}(t) + c_n^{(1)} (t) + c_n^{(2)} (t) + .....          \tag{2} $
with the superscripts meaning the order of the correction.
My question is as follows. In writing $(1)$ we are assuming that the vector space spanned by the eigenkets of the Hamiltonian are going to be the Hilbert space occupied by the system after the perturbation as well. Why is this always justified?
It seems to me that if the perturbation sort of enables new degrees of freedom, then the Hilbert space must necessarily change. For instance, if my $H^{(0)}$ is for a $1D$ SHO along Y, but the perturbation is some driving force along Y. Wouldn't such a perturbation necessarily change the Hilbert space that the problem is now set in? In such a case is $(1)$ still valid?
 A: The perturbation theory, like any approximation technique, has its range of applicability. A perturbation that changes the Hilbert space of the system would likely violate the assumptions of the perturbation theory.
Kondo effect and the Anderson orthogonality catastrophe
In practice the applicability of PT is usually judged on a case-by-case basis, and it is its failure that signals that one needs to apply more sophisticated methods. A well-known example is the Kondo effect - Jun Kondo calculate the probability of scattering from magnetic impurities up to the third order in PT and discovered that in this order the PT diverges (although the 2nd order didn't contain any anomalies). This spurred several decades of activity, with notable contributions from the people such as Phil Anderson and Abrikosov, and a range of techniques including renormalization group, large-N expansion, bosonization, and finally the exact solutions of both Kondo and Anderson models using the Bethe ansatz. The core of the problem is the so-called Anderson orthogonality catastrophe, when the ground state of the perturbed system is orthogonal to the ground state of the unperturbed Hamiltonian - this is probably similar to what the OP had in mind.
A: Already when writing your equation (0), you assume that $H$, $H^{(0)}$ and $\delta H(t)$ are all operators on the same Hilbert space (otherwise, how could you take their sum).

For instance, if my $H^{(0)}$ is for a $1D$ SHO along Y [I assume you meant X here], but the perturbation is some driving force along Y. Wouldn't such a perturbation necessarily change the Hilbert space that the problem is now set in?

If your $H^{(0)}$ describes a 1D harmonic oscillator, I would assume that $H^{(0)} = \frac{p_x^2}{2m} + \frac{m\omega^2}{2} x^2$ without a kinetic term in $y$-direction $\frac{p_y^2}{2m}$. A perturbation with driving force in $y$-direction would thus have to include such a term, which is not small. The question is therefore a bit problematic.
Let us consider instead a 1D harmonic oscillator and a perturbation that depends on some other additional degree of freedom, e.g. spin. For example, assume $H^{(0)} = \frac{p_x^2}{2m} + \frac{m\omega^2}{2} x^2$ as above and $\delta H = \lambda\, x \otimes \sigma_z$ or so. (I don't think it makes a difference for your question whether $\delta H$ is time-dependent or not.)
In this case, it looks at first as if $H^{(0)}$ was an operator acting on one-dimensional position space $L^2(\mathbb R)$ only, but $\delta H$ is acting on the larger space $L^2(\mathbb R) \otimes \mathbb C^2$. We have to interpret $H^{(0)}$ as acting on $L^2(\mathbb R) \otimes \mathbb C^2$ as well, that is, the operator is actually
$$ H^{(0)} = \bigl( \frac{p_x^2}{2m} + \frac{m\omega^2}{2} x^2 \bigr) \otimes 1_{\mathbb C^2} . $$
It can be a bit confusing since this part is usually left implicit. Just from seeing "$H^{(0)} = \frac{p_x^2}{2m} + \frac{m\omega^2}{2} x^2$" we can not determine whether it is supposed to act on $L^2(\mathbb R)$ or $L^2(\mathbb R) \otimes \mathbb C^2$, that must be inferred from the context. Your point of view should however be that it was already the same operator before the perturbation was added.
(And yes, the "extended" $H^{(0)}$ then has twice as many eigenstates, they are $|n, \uparrow\rangle$ and $|n,\downarrow\rangle$ where $|n\rangle$ are the eigenstates of just the 1D harmonic oscillator.)
