When $f \in L^1(B_R)$ we may verify that the suggested function $$u(x) = - \frac{1}{2\pi}\int_{B_R} f(y) \log |x-y|\,dy$$ is a distributional solution at the very least. Note that $u \in L^1_{\text{loc}}(\mathbb{R}^2)$, indeed for any $B_r(x_0) \subset \mathbb{R}^2$ by splitting the integral into $B_R\setminus B_{r+1}(x_0)$ and $B_{r+1}(x_0)$ we have \begin{align} \int_{B_r(x_0)} |u(x)| \,dx & \le \frac{1}{2\pi}\int_{B_R \setminus B_{r+1}(x_0)} \left(\int_{B_r(x_0)} \log |x-y| \,dx\right) |f(y)| \,dy \\& \quad + \frac{1}{2\pi} \int_{B_{r+1}(x_0)} \left(\int_{B_r(x_0)} |\log |x-y|| \,dx\right) |f(y)|\,dy = I_1 + I_2.\end{align} For the first integral $I_1$ note that for $y \not\in B_{r+1}(x_0)$ the function $x \mapsto \log |x-y|$ is a positive harmonic in $B_r(x_0)$ therefore by mean value property we have $$\int_{B_r(x_0)} \log |x-y| \,dx = |B_r|\log |x_0 - y| < |B_r|\log (|x_0| + R).$$ As for the second integral $I_2$ for $y \in B_{2r+1}(x_0)$ we have $B_{r}(x_0) \subset B_{|x_0 - y| + r}(y) \subset B_{2r+1}(y)$ so that $$\int_{B_r(x_0)}|\log |x-y||\,dx \le \int_{B_{|x_0 - y| + r}(y)} |\log |x - y||\,dx \le \int_{B_{2r+1}} |\log |x||\,dx < +\infty.$$

Now it's straightforward to verify that $u$ is a distributional solution, by Fubini's theorem we may interchange the order of integration in the usual way to show that $$\int_{\mathbb{R}^2} \left(-\int_{B_R} \frac{1}{2\pi}\log |x-y| f(y) \,dy\right) \Delta \varphi (x) \,dx = \int_{\mathbb{R}^2} f \varphi \,dx$$ for all $\varphi \in C_c^\infty(\mathbb{R}^2)$.

N.B.: The $\log 2R$ is just there to make $\overline{u}$ positive with $2R \ge |x-y|$ for $x,y \in B_R$. You can refer to some standard Potential theory book for the fact that distributional sub/super-harmonic functions actually equal a.e. to a u.s.c. (resp. l.s.c.) sub/super-harmonic function so that the stated maximum principle in the paper may apply (otherwise let me know, I'll try to dig up some reference).