The problem can be regarded as a special case of incidence reporting problem with $N$ unit circles and $N$ points in the plane.
If the number of unit circles is $m$, then this problem can be solved in $\tilde{O}(m^2 + N)$ time ($\tilde{O}(\cdot)$ ignores polylogarithmic factors). First, explicitly computing the arrangement (the division of the plane, its edges, and vertices) of circles, then constructing a point location data structure, and query the data structure for each point.
The set of $N$ circles can be partitioned into $\sqrt N$ subsets of each $\sqrt N$ size, and each sub-problem can be solved separately. The total runtime is $\tilde{O}(N^{3/2})$. A more involved approach yields a $\tilde{O}(N^{4/3})$-time randomized algorithm [1].