Let $f$ be a newform (normalized eigenform) of weight $k$ and level $\Gamma_0(N)$. Fix $p$ not dividing $N$ and set $f_p(z)=f(pz)$.
Viewing $f$ and $f_p$ at level $\Gamma_0(pN)$, why are they linearly independent?
Here is what I'm thinking: Suppose there is some nontrivial relation $$ af+bf_p=0. \tag{1}$$ The Hecke eigenvalues for $f$ are the same at level $pN$, except possibly at $p$. Thus, for any $\ell\neq p$, applying $T_\ell-a_\ell(f)$ to equation (1) shows that $T_\ell f_p=a_\ell(f)f_p$. Thus, $f$ and $f_p$ have the same Hecke eigenvalues outside $p$.
Now, at this point, I have heard mumblings that the "multiplicity one theorem" shows that $f=f_p$, but I can't seem to find a suitable reference. Simply Googling "multiplicity one" in this context leads to many high-brow results... Diamond and Shurman do not state the multiplicity one theorm in their book (my main reference thus far), but they do cite Miyake's book as a reference. I looked up the theorem there (I believe it is Theorem 4.6.12), and it seems to require that $f$ and $f_p$ both be of level $N$, which (I am fairly certain) they are not. Can someone clarify this for me and possibly provide a good reference for the classical multiplicity one theorem?
From the Fourier series it is obvious that $f \ne cf(pz)$.
$$f(z)=\sum_n a(n) e^{2i\pi n z},\qquad f(pz)=\sum_n a(n) e^{2i\pi np z}=\sum_n a(\frac{n}p)1_{p|n} e^{2i\pi n z}$$
The sequence $a(n)$ is not a scalar multiple of $a(\frac{n}p)1_{p|n}$.
The Fourier series is uniquely defined by the function : $$a(n)=\int_{iy}^{iy+1}f(z)e^{-2i\pi nz}dz$$