Since the original paper by Coifman, Lions, Meyer, and Semmes is not easily accessible--I had to physically scan the original journal--I thought I would include in this thread the simplest proof of the result in question. Let $E\in L^{p}(\mathbb{R}^{n}\rightarrow\mathbb{R}^{n})$ and $B\in L^{p'}(\mathbb{R}^{n}\rightarrow\mathbb{R}^{n})$, where $\frac{1}{p}+\frac{1}{p'}=1$ and $1, be div- and curl-free vector fields respectively.
Lemma. Let $E$ and $B$ be as above. Let $h\in C_{c}^{\infty}(\mathbb{R}^{n})$ (wlog supported in unit ball). For all $\alpha,\beta$ such that $\dfrac{1}{\alpha}+\dfrac{1}{\beta}=1+\dfrac{1}{n}, \quad 1\leq \alpha\leq p, \enspace 1<\beta\leq p'$ there exists a constant $C=C(\alpha,\beta)>0$ such that $|\left(h_{t}\ast (E\cdot B)\right)(x)|\leq C\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|E|^{\alpha}\right)^{1/\alpha}\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|B|^{\beta}\right)^{1/\beta},$ $\forall x\in\mathbb{R}^{n}, t>0$.
Proof of Lemma. By density, we may assume briefly that $B$ is Schwartz; this hypothesis is just to justify the following manipulations. Since $B$ is curl-free, a little Fourier analysis ($n\geq 2$) shows that $\widehat{B}=\xi(\dfrac{\xi}{|\xi|^{2}}\cdot\widehat{B})=\widehat{\nabla\phi},\quad \widehat{\phi}:=\dfrac{\xi}{(2\pi i)|\xi|^{2}}\cdot\widehat{B}$ in the sense of tempered distributions.
I claim that $\phi\in L^{p*}(\mathbb{R}^{n})$, if $1, and $\phi\in L_{loc}^{q}(\mathbb{R}^{n})$ for all $q<\infty$, if $p'\geq n$. Indeed, $\frac{\xi}{|\xi|}$ is, up to a multiplicative constant, the symbol of the vector Riesz transform hence bounded on $L^{r}(\mathbb{R}^{n}\rightarrow\mathbb{R}^{n})$, for all $1. $|\xi|^{-1}$ is, up to a multiplicative constant, the Fourier transform of $|x|^{-n+1}$. The claim then follows from the Hardy-Littlewood-Sobolev lemma. If $n=1$, then just use FTC.
Next, I claim that $E\cdot B=\text{div}(E\phi)$ in $\mathcal{D}'$. Formally, this is obvious (just use product rule and the hypothesis that $\text{div}(E)=0$) but to prove this rigorously we use mollification. Let $\psi$ be a $C_{c}^{\infty}$ mollifier. Then \begin{align*} \text{div}((E\ast\psi_{\varepsilon})(\psi_{\varepsilon}\ast\phi))&=\text{div}(E\ast\psi_{\varepsilon})(\psi_{\varepsilon}\ast\phi)+(E\ast\psi_{\varepsilon})\cdot ((\nabla\phi)\ast\psi_{\varepsilon})\\ &=\text{div}(E)\ast\psi_{\varepsilon}+(E\ast\psi_{\varepsilon})\cdot(B\ast\psi_{\varepsilon})\\ &=(E\ast\psi_{\varepsilon})\cdot(B\ast\psi_{\varepsilon})\\ \end{align*} Since $E\ast\psi_{\varepsilon}\rightarrow E$ in $L^{p}$ and $B\ast\psi_{\varepsilon}\rightarrow B$ in $L^{p'}$, the last expression above $\rightarrow E\cdot B$, as $\varepsilon\rightarrow 0$, in $L^{1}$, hence in distribution. Similarly, the LHS $\rightarrow\text{div}(E\phi)$ in $\mathcal{D}'$.
So, we may write \begin{align*} (h_{t}\ast(E\cdot B))(x)&=\int_{\mathbb{R}^{n}}(\nabla h)(\frac{x-y}{t})t^{-n-1}\cdot E(y)\phi(y)dy\\ &=\int_{\mathbb{R}^{n}}(\nabla h)(\frac{x-y}{t})t^{-n-1}\cdot E(y)\left[\phi(y)-\frac{1}{|B_{t}(x)|}\int_{B_{t}(x)}\phi\right]dy \end{align*} where the last line follows by the hypothesis that $\text{div}(E)=0$ and $|B_{t}|^{-1}\int_{B_{t}(x)}\phi$ is a constant. Using H\"{o}lder's inequality with $\beta$, we obtain \begin{align*} |(h_{t}\ast (E\cdot B))(x)|\lesssim_{n}\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|E|^{\beta}dy\right)^{1/\beta}\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}\left|[\phi-\frac{1}{|B_{t}(x)|}\int_{B_{t}(x)}\phi]t^{-1}\right|^{\beta'}dy\right)^{1/\beta'} \end{align*} Using the Poincare inequality, we can estimate the second factor by $\lesssim_{n,\beta,\alpha}\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|\nabla \phi|^{\alpha}\right)^{1/\alpha}$ since $\frac{1}{\alpha}-\frac{1}{n}=1-\frac{1}{\beta}=\frac{1}{\beta'}$ (i.e. $\beta'$ is the Sobolev conjugate of $\alpha$). Recalling that $\nabla\phi=B$ completes the proof of the lemma. $\Box$
Proof of Main Resuklt. We use the maximal function characterization of Hardy spaces with $h$ as our Schwartz function to show that $E\cdot B\in\mathcal{H}^{1}(\mathbb{R}^{n})$. Since $\frac{1}{p}+\frac{1}{p'}=1$, we can find $\alpha, $\beta as in the statement of the lemma. Whence, \begin{align*} \sup_{t>0}\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|E|^{\alpha}\right)^{1/\alpha}\left(\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|B|^{\beta}\right)^{1/\beta}&\leq\left(\sup_{t>0}\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|E|^{\alpha}\right)^{1/\alpha}\left(\sup_{t>0}\dfrac{1}{|B_{t}(x)|}\int_{B_{t}(x)}|B|^{\beta}\right)^{1/\beta}\\ &\leq (M(|E|^{\alpha})(x))^{1/\alpha}(M(|B|^{\beta})(x))^{1/\beta}, \end{align*} where $M$ denotes the Hardy-Littlewood maximal function. Observe that $|E|^{\alpha}\in L^{p/\alpha}(\mathbb{R}^{n})$ and $|B|^{\beta}\in L^{p'/\beta}(\mathbb{R}^{n})$, where $p/\alpha, p'/\beta>1$. Since $M$ is bounded on $L^{r}(\mathbb{R}^{n})$, for $1, we conclude from Holder's inequality that \begin{align*} \int_{\mathbb{R}^{n}}(M(|E|^{\alpha}))^{1/\alpha}(M(|B|^{\beta}))^{1/\beta}dx&\leq\|(M(|E|^{\alpha}))^{1/\alpha}\|_{L^{p}}\|(M(|B|^{\beta}))^{1/\beta}\|_{L^{p'}}\\ &=\|M(|E|^{\alpha})\|_{L^{p/\alpha}}^{1/\alpha}\|M(|B|^{\beta})\|_{L^{p'/\beta}}\\ &\lesssim_{n,p,\alpha}\|E\|_{L^{p}}\|B\|_{L^{p'}}, \end{align*} which completes the proof.