Elliptic regularity of harmonic forms in $L^1$ $\newcommand{\M}{M}$
This is a cross-post. I am looking for a reference for the regularity of harmonic forms which belong to $L^1(M)$.
Explicitly, let $\M$ be a smooth oriented Riemannian manifold.
Let $\sigma$ be a differential $k$-form on $\M$ with coefficients in $L^1(\M)$.  We say $\sigma$ is weakly harmonic if 
$$
\int_{\M} \langle \sigma , \Delta \alpha \rangle=0 \, \text{ for every smooth compactly supported $\alpha \in \Omega^{k}(\M)$}. 
 $$
Theorem: Weakly harmonic forms (in $L^1$) are smooth.
I know a reference for this theorem, when we assume $\sigma \in L^2$ instead of $L^1$ (Warner's book). I am looking for a reference which deals with the $L^1$ case, which naively looks harder.
(Or, if possible, for an easy way to reduce the $L^1$ theorem to the $L^2$ theorem).
If that matters, I don't care about global topology here, i.e. we can assume $M=\mathbb{R}^n$ topologically. (but not "Riemannialy", that is, I don't assume the metric is Euclidean).
I only care about the smoothness, not about any generalized version of Hodge decomposition.

In my case, I actually know that $\sigma \in L^1$ is weakly closed and weakly co-closed (which is stronger than being weakly harmonic, if $M$ is not closed), but I am not sure if it should make the proof easier. For completeness, these are the definitions I am using for these properties: 
$$
 \text{ weakly closed if } \,  \int_{\M} \langle \sigma, \delta \alpha \rangle =0 \, \text{ for every compactly supported $\alpha \in \Omega^{k+1}(\M)$}, 
$$
$$
\text{ weakly co-closed if } \, \int_{\M} \langle \sigma ,d \alpha \rangle=0 \, \text{ for every compactly supported $\alpha \in \Omega^{k-1}(\M)$}, 
$$
 A: Let me turn my comment into a general discussion for which a special case is an answer to your question.
First, on an open domain $D \subset \mathbb{R}^n$, there is a standard elliptic regularity result that says if $u$ is a distribution on $D$ satisfying weakly $$ a^{ij}\partial^2_{ij}u + b^k\partial_ku + cu= f, $$ where $a^{ij}, b^i, c, f$ are assumed to be smooth and, for any $x \in D$ and $\xi \in \mathbb{R}^n\backslash\{0\}$, $a^{ij}(x)\xi_i\xi_j \ne 0$, then $u$ is smooth. The easiest way to prove this is by constructing a parametrix using pseudodifferential operators. Details can be found in the book by Treves, Introduction to Pseudodifferential and Fourier Integral Operators, as well as the one by Chazarain and Piriou, Introduction to the Theory of Linear Partial Differential Equations. I'm surprised that I can't find a more recent exposition of this, since I learned it all from these books over 30 years ago.
Second, the exact same proof still works, if $u, f$ are vector-valued and $a^{ij}, b^i, c$ are matrix-valued and, for any $x \in D$ and $\xi \in \mathbb{R}^n\backslash\{0\}$, $a^{ij}(x)\xi_i\xi_j$ is invertible. This may or may not be stated explicitly in the books above but is easily verified.
Finally, using local coordinates and a partition of unity on a manifold and local trivializations of a vector bundle, the regularity theorem extends to a bundle-valued distribution $u$ satisfying $$Lu = f,$$ where $L$ is an elliptic PDO and $f$ is a smooth section of the bundle. Based on the table of contents, Treves' book appears to discuss this.
In your specific situation, a harmonic differential form $u$ satisfies such a PDE with $f = 0$. Moreover, if $u$ is in $L^1$, it is a distribution. Therefore, $u$ is smooth.
I'm also pretty sure Warner's proof can be adapted to prove what you want. His proof is essentially the same as the pseudodifferential one. However, the latter is actually easier to learn and understand, even though it might not seem so at first.
