How is solving Proca equation equivalent to scalar field equation? My prof. told me that using differential forms proca equation reduces to solving for scalar field equation. How is that? I can’t see how does one relate to Scalar equation using differential forms.
Proca equation: $$\mathcal{L} = \frac{-1}{16}F^{\mu v}F_{\mu v} + \frac{1}{8\pi}m^2A_\mu A^\mu.$$
Equation of motion for Proca: $$\partial_\mu F^{\mu v} + m^2 A^v = 0.$$
 A: They reduce to 4 Klein-Gordon equation Because, as noted in a comment by the OP, the 4-divergence of $A^\nu$ is identically 0, this follow directly from the equation of motion. So we have
$$
\partial_\mu A^\mu = 0 \tag{1}
$$
Now, we start from the equation of motion you propose
$$
\begin{aligned}
\partial_\mu F^{\mu\nu} + m^2 A^\nu = 0  
\end{aligned}
\tag{2}
$$
We insert the definition of the tensor $F^{\mu\nu} = \partial^\mu A^\nu - \partial^\nu A^\mu $ and we obtain
$$
\partial_\mu (\partial^\mu A^\nu - \partial^\nu A^\mu) + m^2 A^\nu  = 0 \tag{3} 
$$
From there we get
$$
\partial^\mu \partial_\mu A^\nu - \partial_\mu \partial^\nu  A^\mu + m^2A^\nu = \Box A^\nu -  \partial^\nu \partial_\mu  A^\mu + m^2A^\nu = 0 \tag{4}
$$
Where in the first equality in equation $(4)$, in the second term, I just swapped the derivatives $\partial_\mu$ with $\partial^\nu$
Now if we use $(1)$ in $(4)$ the second term vanishes and we obtain
$$
(\Box + m^2)A^\nu  = 0 \tag{5}
$$
Which is just a set of 4 Klein-Gordon equations, one for each component of the $4-$vector $A^\nu$. 
Notice that those are not $4$ equation for $4$ scalar fields, but for the components of a $4-$ vector.
A: The Lagrangian that you use for the Proca equation looks a bit unusual, I will factor out $\frac{1}{8}$ and change the factor  ($\frac{m}{\sqrt{\pi}} \rightarrow \frac{m}{\sqrt{2}}$): Then written out in differential forms it looks like:
$$L = 8{\cal{L}} =-\frac{1}{2} dA \wedge \star dA + \frac{1}{2} m^2   A\wedge \star A$$
We will derive the L-E equations via variation (actually I did most of it already in another post, so I will shorten the derivation a little bit, the details can be looked up in the post https://physics.stackexchange.com/a/432941/30506 ):
$$\delta S =\int\delta L = -\frac{1}{2} \int (d\delta A \wedge \star dA + dA \wedge \star d\delta A) +  \frac{1}{2} m^2  \int (\delta A \wedge \star A + A\wedge \star \delta A) =   -\frac{1}{2} \int 2 d\delta A \wedge \star dA + m^2 \int \delta A \wedge \star A $$
In the next manipulation we will use the product rule: 
$$d(\delta A\wedge \star d A) = d\delta A \wedge\star dA - \delta A \wedge d\star dA $$ therefore we have:  $$-d\delta A\wedge \star dA = - d(\delta A\wedge \star dA ) -\delta A\wedge d\star dA $$.
We will this substitute in the first term of the varied action:
$$\delta S =\int  - d(\delta A\wedge \star dA )  -\int  \delta A\wedge d\star dA  + m^2  \int \delta A \wedge \star A $$
Finally the first term is an integral over an absolute derivative which be transformed into a surface integral on whose surface the variation $\delta A=0$.   So we get finally:
$$0=\delta S  = - \int \delta A \wedge  (d\star dA  - m^2  \star A )$$
As the last expression has to be zero for all variations  $\delta A$, the result of the variation is:
$$d\star dA  - m^2   \star A =0 $$  or a bit more nicely written (actually at the moment I don't know if $\star \star =1$ or $\star \star =-1$, but I will check that up).
$$ \star d\star dA = m^2 A$$ 
In some books $ \delta: =  \star d\star$ (this $\delta$, however, has nothing to do with the variation) and with $F=dA$ we get:
$$\delta F =m^2 A $$
The definition of the hodge operator can be looked up in the other post  https://physics.stackexchange.com/a/432941/30506  I mentioned at the beginning. 
