Is the rectangular function a convolution of $L^1$ functions?

Yes, this works. Let's change the interval to $I=[-1,1]$ for convenience. Then we want to find $f,g\in L^1$ such that $$ \widehat{f}(t)\widehat{g}(t) = \widehat{\chi_I}=\frac{\sin t}{t} . $$ Let's take $\widehat{f}(t) = |t|^{-1/2}$ for $|t|\ge 1$ and then of course $\widehat{g}$ such that $\widehat{f}\widehat{g}=\widehat{\chi_I}$ (so in particular, $\widehat{g}(t)=\textrm{sgn}(t)\,|t|^{-1/2}\sin t$ for $|t|\ge 1$), and then finally, I'm also going to insist that $\widehat{f}, \widehat{g}\in C^{\infty}$.

To show that $f,g\in L^1$, let's first of all discuss the local behavior. This only depends on the large $|t|$ behavior of the Fourier transform, and $|t|^{-1/2}$ has FT $c|x|^{-1/2}$, which is locally integrable (but not locally $L^2$, which is just as well since it's easy to see that at least one of $f,g$ will fail to be in $L^2$ if $f*g=\chi_I$).

The large $x$ asymptotics of $f(x)$ on the other hand depend on the smoothness of $\widehat{f}$, and here we're clearly in great shape since $\widehat{f}''\simeq t^{-5/2}\in L^1$, so $x^2f(x)\in C_0$.

The other function is of the same type, the extra factor $\sin t$ for large $t$ produces essentially a shift in $g$, and $\textrm{sgn}(t)$ amounts to an extra Hilbert transform, which won't affect the local singularities here, so $g\in L^1$ also.

Following the idea of Christian Remling, I have two explicit functions that have the property that $f*g = \mathbf{1}_{[-1,1]}$: \begin{eqnarray} f(x) &=& \frac{K_\frac{1}{4}\left(|x|\right)}{2^{1/4}\pi\Gamma\left(\frac{1}{4}\right)|x|^{1/4}}\\\\ g(x) &=& \frac{\mathrm{sgn}(1-x)}{\sqrt{|1-x|}}\,_1F_2\left(\begin{matrix}-\frac{1}{4} \\ \frac{1}{4}, \frac{3}{4}\end{matrix}; \frac{(1-x)^2}{4}\right) +\frac{\mathrm{sgn}(1+x)}{\sqrt{|1+x|}}\,_1F_2\left(\begin{matrix}-\frac{1}{4} \\ \frac{1}{4}, \frac{3}{4}\end{matrix}; \frac{(1+x)^2}{4}\right)\\ &\,&+ \frac{\Gamma\left(\frac{5}{4}\right)(1-x)}{\sqrt{2}\Gamma\left(\frac{7}{4}\right)}\,_1F_2\left(\begin{matrix}\frac{1}{2} \\ \frac{3}{2}, \frac{7}{4}\end{matrix}; \frac{(1-x)^2}{4}\right)+ \frac{\Gamma\left(\frac{5}{4}\right)(1+x)}{\sqrt{2}\Gamma\left(\frac{7}{4}\right)}\,_1F_2\left(\begin{matrix}\frac{1}{2} \\ \frac{3}{2}, \frac{7}{4}\end{matrix}; \frac{(1+x)^2}{4}\right) \end{eqnarray} Where $K_\frac{1}{4}$ is a modified Bessel function and $\,_1F_2$ is a hypergeometric function. While these aren't the simplest functions, and $g$ is kind of a nightmare, they have the nice property that \begin{eqnarray} \mathcal{F}[f](k) &=& \frac{1}{\sqrt{2\pi}(1+k^2)^{1/4}} \\ \mathcal{F}[g](k) &=& 2\sqrt{2\pi}(1+k^2)^{1/4}\frac{\sin(k)}{k}, \end{eqnarray} so that $f*g = \mathcal{F}^{-1}[\mathcal{F}[f]\cdot\mathcal{F}[g]] = \mathcal{F}^{-1}\left[2\sin(k)/k\right] = \mathbf{1}_{[-1,1]}$.

The convolution integral isn't doable analytically, but it can be evaluated numerically. Doing this in Mathematica gives

where the blue is $f$, the yellow is $g$, and the green is $f*g$. Despite their divergences, both $f$ and $g$ are in $L^1$, as the divergences scale as $|x|^{-1/2}$.

By the Cohen–Hewitt factorization theorem, every $f\in L^1(\mathbb R)$ can be written as $g\star h$ for some $g,h\in L^1(\mathbb R).$