Let \(d\geqslant 2\) . In this paper, we investigate the following stochastic differential equation (SDE) in \({{\mathbb {R}}}^d\) driven by Brownian motion \( \textrm{d} X_t=b(t,X_t)\textrm{d} t+\sqrt{2}\textrm{d} W_t, \) where b belongs to the space \({{\mathbb {L}}}_T^q \textbf{H}_p^\alpha \) with \(\alpha \in [-1, 0]\) and \(p,q\in [2, \infty ]\) , which is a distribution-valued and divergence-free vector field. In the subcritical case \(\frac{d}{p}+\frac{2}{q}<1+\alpha \) , we establish the existence and uniqueness of a weak solution to the integral equation: \( X_t=X_0+\lim _{n\rightarrow \infty }\int ^t_0b_n(s,X_s)\textrm{d} s+\sqrt{2} W_t. \) Here, \(b_n:=b*\phi _n\) represents the mollifying approximation, and the limit is taken in the \(L^2\) -sense. In the critical and supercritical case \(1+\alpha \leqslant \frac{d}{p}+\frac{2}{q}<2+\alpha \) , assuming the initial distribution has an \(L^2\) -density, we show the existence of weak solutions and associated Markov processes. Moreover, under the additional assumption that \(b=b_1+b_2+\mathord {\textrm{div}}a\) , where \(b_1\in {{\mathbb {L}}}^\infty _T{{\textbf{B}}}^{-1}_{\infty ,2}\) , \(b_2\in {{\mathbb {L}}}^2_TL^2\) , and a is a bounded antisymmetric matrix-valued function, we establish the convergence of mollifying approximation solutions without the need to subtract a subsequence. To illustrate our results, we provide examples of Gaussian random fields and singular interacting particle systems, including the two-dimensional vortex models.
