  In this paper, we introduce and analyze arbitrarily high-order quadrature
rules for evaluating the two-dimensional singular integrals of the forms
\begin{align}
  I_{i,j} = \int_{\mathbb{R}^2}\phi(x)\frac{x_ix_j}{|x|^{2+\alpha}} \d x, \quad
0< \alpha < 2
  \end{align} where $i,j\in\{1,2\}$ and $\phi\in C_c^N$ for $N\geq 2$. This
type of singular integrals and its quadrature rule appear in the numerical
discretization of fractional Laplacian in non-local Fokker-Planck Equations in
2D. The quadrature rules are trapezoidal rules equipped with correction weights
for points around singularity. We prove the order of convergence is
$2p+4-\alpha$, where $p\in\mathbb{N}_{0}$ is associated with total number of
correction weights. Although we work in 2D setting, we formulate definitions
and theorems in $n\in\mathbb{N}$ dimensions when appropriate for the sake of
generality. We present numerical experiments to validate the order of
convergence of the proposed modified quadrature rules.
