Let $f= \sum_{n=1}^{\infty}a_nq^n \in S_2^{new}(N)$ be a normalized new form of weight $2$ with respect to $\Gamma_0(N)$ and assume $p|N$ is a prime. Then must $a_p=0$ if $p^2|N$ and belongs to $\{-1,1\}$ if $p|N$?
If $f$ has rational coefficients, then this can be seen from the corresponding modular elliptic curve as $N$ is the conductor.
So here I don't assume $f$ has rational coefficients so then the definition field $K_f$ may be large than $\mathbb Q$. Then again there is a modular abelian variety, but it's subtle because the $L$-function of the abelian variety matches the product of $L$-functions of Galois conjugation of $f$ rather than $f$, and I don't know much about Neron model.
Can it be proved only using knowledge of modular forms? I am also happy with a proof involving algebraic geometry of abelian varieties.
Yes, this is a standard property of modular forms which you can prove "by hand" using a computation with double cosets; suitably stated, the property holds for all weights, including weight 1 (whereas modular abelian varieties are a weight 2 thing).
In the $p^2 \mid N$ case, the idea is to check that the double coset $\Gamma_0(N) \begin{pmatrix} 1 & 0 \\ 0 & p \end{pmatrix} \Gamma_0(N)$ giving the Hecke operator $U_p$ is actually stable under right-multiplication by $\Gamma_0(N/p)$, so if $f \in S_k(N)$, then $f \mid_k U_p$ is actually in $S_k(N/p)$. However, $f \mid_k U_p$ also has the same Hecke eigenvalues away from $p$ as $f$, so if $f$ is new of level $N$, it had better be 0.
For $p \mid\mid N$ this needs to be modified slightly because the index of $\Gamma_0(N)$ in $\Gamma_0(N/p)$ is $p+1$ instead of $p$. This gives an extra term in the formula, and you end up deducing that $f\mid_k U_p + f \mid_k w_p$ is zero, where $w_p$ is the Atkin-Lehner operator. Since $w_p$ is an involution, its eigenvalue had better be $\pm 1$ and this gives the result. (For general weights $k$ there is an extra normalisation factor coming out here, and you get that the $U_p$ eigenvalue is $\pm p^{(k-2)/2}$ instead.)
For newforms of $\Gamma_1(N)$ levels instead of $\Gamma_0(N)$ levels, there is a slightly more complicated statement, where you have to keep track of the powers of $p$ dividing both $N$ and the conductor of the character of $f$.