Is there a constant $c>0$ such that for $X,Y$ two iid variables supported by $[0,1]$, $$ \liminf_\epsilon \epsilon^{1}P(XY<\epsilon)\geqslant c $$ I can prove the result if they have a density, of if they have atoms, but not in the general case.
If $\epsilon \geqslant \tfrac{1}{n}$, then $$ \mathbb{P}(XY<\epsilon) \geqslant \sum_{i=1}^n \mathbb{P}(X, Y \in [\tfrac{i1}{n}, \tfrac{i}{n}]) = \sum_{i=1}^n (\mathbb{P}(X \in [\tfrac{i1}{n}, \tfrac{i}{n}]))^2 . $$ It follows that $$ \mathbb{P}(XY<\epsilon) \geqslant \frac{1}{n} \biggl(\sum_{i=1}^n \mathbb{P}(X \in [\tfrac{i1}{n}, \tfrac{i}{n}])\biggr)^2 = \frac{1}{n} \, . $$ If we choose $n$ so that additionally $\epsilon < \frac{1}{n1}$, then we obtain $$ \mathbb{P}(XY<\epsilon) > \frac{\epsilon}{\epsilon + 1} \, . $$ This leads to the desired result with $c = 1$.

2$\begingroup$ Interestingly, while the estimate given above is sharp at least when $\epsilon = \tfrac{1}{n}$ (consider a uniform distribution over $\frac{i}{n  1}$, $i = 0, 1, \ldots, n  1$), my guess is that the optimal value of $c$ is in fact $2$. However, I fail to see a proof. $\endgroup$ – Mateusz Kwaśnicki Mar 10 at 21:44
Mateusz Kwaśnicki's guess, that the optimal value of $c$ is 2, is correct! In fact:
Theorem: Suppose $X,Y$ are i.i.d. random variables on $\mathbb{R}$. Then the limit $\lim_{\varepsilon \rightarrow 0} \varepsilon^{1} \mathbb{P}[X  Y < \varepsilon]$ exists. It is $+\infty$ unless $X$ has a density function $f$ satisfying $f_2^2 = \int f^2\ \text{d}x < \infty$, in which case the limit is equal to $2 f_2^2$.
Restricting to distributions on $[0,1]$, CauchySchwarz says $$ 1 = \int_0^1 f\ \text{d}x \le \left ( \int_0^1 f^2\ \text{d}x \right )^{1/2} \left ( \int_0^1 1\ \text{d}x \right )^{1/2} = f_2, $$ with equality iff $f = 1$. That is, the constant $c = 2$ is valid for every distribution, and it is sharp only for the uniform distribution.
Here are a few examples where the limit is infinite:
Example 1: If $X$ has an atom then $\mathbb{P}[XY < \varepsilon]$ is bounded away from 0. So the limit in question diverges like $\varepsilon^{1}$.
Example 2: The Cantor ternary function is the CDF of a probability distribution on $[0,1]$ which is nonatomic but not absolutely continuous. A random sample from this distribution is given by $\sum_{k \ge 1} c_k \tfrac{2}{3^k}$ where $\{c_k\}$ are i.i.d. Bernoulli random variables. (Thought of in terms of the "repeatedly take out the middle third" construction of the Cantor set, the $c_k$'s are the choices of left versus right third.) If $\varepsilon = 3^{n}$, then $XY < \varepsilon$ iff the first $n$ $c_k$'s agree. Thus $$ \mathbb{P}[XY < \varepsilon] = 2^{n} = \varepsilon^{\log_3 2}. $$ So the limit in question diverges like $\varepsilon^{1+\log_3 2} = \varepsilon^{0.37}$.
Example 3: The power law $f(x) = \tfrac{1}{2\sqrt{x}}$ is an example of a density function which is not in $L^2$. In this case we can evaluate $\mathbb{P}[XY < \varepsilon]$ exactly (well, if you believe in inverse hyperbolic trig functions). The asymptotic behavior is $$ \mathbb{P}[XY < \varepsilon] = \tfrac{1}{2} (\log \varepsilon) \varepsilon + (\tfrac{1}{2} + \log 2) \varepsilon + O(\varepsilon^2). $$ So the limit in question diverges logarithmically.
Mateusz's slick argument worked on the block diagonal, a sum of $n$ squares of width $1/n$. This covers about half of the area of the strip $XY < \tfrac{1}{n}$, which is why the resulting constant is half of optimal. There's probably a handson way to extend it, but you start using words like "convolution" and I found it easier to reason in Fourier space. This question is asking for a bound on the density of $XY$ at 0. This random variable has nonnegative Fourier transform (characteristic function, if you prefer), and that condition alone is enough to guarantee positive density (shameless selfcitation: Lemma 3.1 in this paper). In general, when it makes sense, the density at 0 of a random variable is equal to the integral of its Fourier transform. So our job is to make that concrete and translate it back to a statement about PDFs.
Convention: The Fourier transform (characteristic function) of a random variable $X$ is the function $t \mapsto \mathbb{E}[e^{2 \pi \mathrm{i} X t}]$. The Fourier transform of a function $f$ is $t \mapsto \hat f(t) = \int e^{2 \pi \mathrm{i} x t} f(x)\ \text{d}x$. If $X$ is absolutely continuous (has a density function), then its Fourier transform is equal to the Fourier transform of its density function. With this convention the Plancherel identity reads $\int f \bar g\ \text{d}x = \int \hat f \overline{\hat g}\ \text{d}t$ (no normalization constant).
Lemma: Let $\psi$ denote the Fourier transform of $X$. Then $\psi \in L^2$ iff $X$ is absolutely continuous and its density function $f$ is in $L^2$. Moreover, if this holds, then $\psi_2 = f_2$.
Proof of Lemma: The Fourier transform is an isometry of the space $L^2$. The ($\Leftarrow$) direction is immediate: if $X$ has a density function in $L^2$, then certainly its Fourier transform is in $L^2$. For ($\Rightarrow$), if $\psi$ is in $L^2$ then it is the Fourier transform of some function $f \in L^2$. But then $f$ defines the same distribution as $X$, so it is indeed the density function. $\square$
(This is probably a basic result in some probability textbook?)
Proof of Theorem: Let $\psi$ denote the Fourier transform of $X$. Then the Fourier transform of $XY$ is $t \mapsto \psi(t) \psi(t) = \psi(t) \overline{\psi(t)} = \psi(t)^2$. Let $T_\varepsilon$ denote the "triangle filter" $T_\varepsilon(x) = \varepsilon^{1} (1  x/\varepsilon)$ for $x \le \varepsilon$ and $T_\varepsilon(x) = 0$ otherwise. It is a standard calculation that $\hat T_\varepsilon(t) = \operatorname{sinc}^2(\pi t \varepsilon)$. (Here $\operatorname{sinc} x = \tfrac{\sin x}{x}$.) This is integrable, so we can compute the expected value of $T_\varepsilon$ using the Fourier transform. We find $$ \varepsilon^{1} \mathbb{P}[XY < \varepsilon] \ge \mathbb{E}[T_\varepsilon(XY)] = \int \psi(t)^2 \operatorname{sinc}^2(\pi t \varepsilon)\ \text{d}t. $$ The integrand is nonnegative and, as $\varepsilon \to 0$, it converges (pointwise) to its (pointwise) supremum $\psi(t)^2$. So, by an appropriate incantation of the monotone convergence theorem (the one in this comment applies exactly), the integral converges to $\int \psi(t)^2\ \text{d}t = \psi_2^2$. If this is infinite, then also $\lim_{\varepsilon \to 0} \varepsilon^{1} \mathbb{P}[XY < \varepsilon] = +\infty$.
Suppose now that $\psi_2^2 < \infty$, i.e., $\psi \in L^2$. Consider the "box filter" $B_\varepsilon$ defined by $B_\varepsilon(x) = (2 \varepsilon)^{1}$ for $x < \varepsilon$ and $B_\varepsilon(x) = 0$ otherwise. This has $\hat B_\varepsilon(t) = \operatorname{sinc}(2 \pi t \varepsilon)$. This is bounded and $\psi^2$ is integrable. So we can again compute expected values in Fourier space: $$ (2\varepsilon)^{1} \mathbb{P}[XY < \varepsilon] = \mathbb{E}[B_\varepsilon(XY)] = \int \psi(t)^2 \operatorname{sinc}(2 \pi t \varepsilon)\ \text{d}t. $$ The integrand converges pointwise to $\psi(t)^2$. We don't have nonnegativity this time, but now we know integrability of the bounding function $\psi^2$. So we can apply the dominated convergence theorem, concluding that $$ \lim_{\varepsilon \to 0} (2 \varepsilon)^{1} \mathbb{P}[XY < \varepsilon] = \int \psi(t)^2\ \text{d}t = \psi_2^2. $$ The lemma completes the proof. $\square$