Let $M$ be a smooth $d$-dimensional oriented Riemannian manifold, and let $p \in M$. Let $0 \le k \le d$ be fixed.
Does there exist an open neighbourhood $U$ of $p$, which admit a non-zero harmonic $k$-form? i.e $\omega \in \Omega^k(U)$ satisfying $d\omega=\delta \omega=0$?
Since a form $\omega$ is harmonic if and only if $\star \omega$ is harmoinc ($\star$ is the Hodge dual operator), the answers for a given $k,d-k$ are the same.
For $k=d$, one can take $\omega$ to be the Riemannian volume form.
For $k=1$, we can take $\omega=df$ where $f$ is a harmonic function. Then $\delta \omega=\delta df=0$. Locally, there are always harmonic functions- we can solve the Dirichlet problem on a small ball with boundary, that is finding a harmonic function which is zero on the boundary.
This solves the cases $k=0,1,d-1,d$.
So, we are left with the cases $2 \le k \le d-2$.
Yes, locally there is always an infinite-dimensional space of harmonic $k$-forms (as long as the dimension $d$ is at least 2). One way to see this is a variational approach, since harmonic forms are critical points of an energy functional.
Assume that $U \subset M$ is a small ball around $p$ and that $\alpha$ is some smooth $(k-1)$-form defined in a neighborhood of $\partial U$. Let $\beta$ be a $(k-1)$-form minimizing $\int_U \langle d\beta, d\beta \rangle$ with boundary values $\beta = \alpha$ on $\partial U$, and let $\omega = d\beta$. Then $d\omega = 0$, so for harmonicity we only need to show that $\delta \omega = 0$, which turns out to be exactly the Euler-Lagrange equation for this functional: If $\gamma$ is a $(k-1)$-form vanishing on $\partial U$, we get that $$ 0 = \frac{d}{dt}\int_U \langle \omega + td\gamma, \omega + t d\gamma) \rangle = 2 \int_U \langle d\gamma,\omega \rangle = 2 \int_U \langle \gamma, \delta \omega \rangle, $$ which then implies $\delta \omega = 0$. So $\omega$ is a harmonic $k$-form with boundary values $\omega = d\alpha$ on $\partial U$. Since $\alpha$ was an arbitrary smooth $(k-1)$-form in a neighborhood of $\partial U$, the space of boundary values, and thus the space of locally harmonic forms is infinite-dimensional.
Obviously, a lot of details have to be filled in, in particular the smoothness of the minimizer. These are treated in textbooks on elliptic differential equations.