You can try $\mathcal B(\varphi_1,\varphi_2)=\int \int\varphi'_1(x)\ \varphi_2'(y) (|x-y|\ln|x-y|)\ dx\ dy$ , i.e. $K$ is the second derivative $\partial_x\partial_y$ of (the distribution defined by) the locally integrable function $H(x,y)=|x-y|\ln|x-y|$ .
This extension is clearly non unique (since any distribution applied to $\varphi_1\varphi_2$ vanishes if the two functions have disjoint supports) but a somewhat standard solution (different from the one above) would use Schwartz's pseudo-function Pf. $r^m$ defined e.g. in Théorie des distributions, formula (II,3;4), a well-defined distribution whose restriction to $\{0\}^c$ coincides with the locally summable function $|x|^m$. This leads (via the Fourier transform of Pf. $r^{-1}$, formula (VII,7;18)) to $$\mathcal B(\varphi_1,\varphi_2)=-2[\int \hat{\varphi_1}(\xi)\ \hat{\varphi_2}(-\xi)\ln|\xi|\ d\xi\ +\ (\mathcal C+\ln(2\pi))\int\varphi_1\int\varphi_2]$$where $\mathcal C$ is Euler's constant. Or equivalently (through the definition of Pf. $r^{-1}$)$$\mathcal B(\varphi_1,\varphi_2)=\lim_{\varepsilon\to0}[\int\int_{|x-y|>\varepsilon}\frac{\varphi_1(x)\varphi_2(y)}{|x-y|}\ dx\ dy\ -\ I(\varepsilon)]$$where $I(\varepsilon)\propto \ln\varepsilon$ , i.e. the finite part of the possibly diverging integral.