There are at least two interpretations: one concerns the actual points generated by this process and the other concerns the process itself.
If a realization of the Poisson process is given and pairs of points are to be chosen from that realization, then there is nothing to be done except systematically compare all distances to all other distances (a double loop over the points).
Otherwise, if the procedure is intended to consist of (i) creating a realization of the process and then (ii) selecting a pair of points at random, then the assumptions imply the two points are selected uniformly and independently from the circle. The calculation for this situation can be performed once and for all.
Notice that the squared distances $r_1 = d_1^2$ and $r_2 = d_2^2$ are uniformly distributed, whence the desired probability is
$$p(a,b) = \Pr\left(d_1^2 \lt \frac{d_2^2}{a(1 + b d_2^2)}\right) = \int_0^1 d r_2 \int_0^{\max(0, \min(1, r_2 / (a(1 + b r_2))))} d r_1.$$
The $\max$ and $\min$ can be handled by breaking into cases. Some special values of $a$ and $b$ have to be handled. Because the integration is a square window over a region generically bounded by lines and lobes of a hyperbola (with vertical axis at $1/(ab)$ and horizontal axis at $-1/b$), the result is straightforward but messy; it should involve rational expressions in $a$ and $b$ and some inverse hyperbolic functions (that is, natural logarithms). I had Mathematica write it out:
$$\begin{array}{ll}
\frac{b+1}{b} & \left(-1\leq a<0\land \frac{1}{a}-b\leq 1\land b<-1\right)\\ &\lor \left(a<-1\land \frac{1}{a}-b<1\land b<-1\right) \\
-\frac{1}{b (a b-1)} & \frac{1}{a}-b=1\land a<-1 \\
\frac{a^2 b+2 a b+a-2}{2 (a b-1)} & b=0\land a>0\land \frac{1}{a}-b>1 \\
\frac{b-\log (b+1)}{a b^2} & a>0\land \frac{1}{a}-b\leq 1\land b>-1 \\
\frac{a b^2+a b-a b \log (b+1)-b+\log (b+1)}{a b^2 (a b-1)} & a>0\land \frac{1}{a}-b\leq 1\land b\leq -1 \\
\frac{\log (1-a b)}{a b^2} & a>0\land \frac{1}{a}-b>1\land b\leq -1 \\
\frac{a b^2+a b+\log (1-a b)}{a b^2} & \left(-1<b<0\land a>0\land \frac{1}{a}-b>1\right) \\
& \lor \left(b>0\land a>0\land \frac{1}{a}-b>1\right) \\
\frac{b-\log ((-b-1) (a b-1))}{a b^2} & a<0\land \frac{1}{a}-b>1
\end{array}$$
Numeric integration and simulation over the ranges $-2 \le a \le 2$ and $-5 \le b \le 5$ confirm these results.
Edit
The modified question asks to replace $d_i^2$ by $d_i^\alpha$ and assumes $a$ and $b$ both positive. Upon making a substitution $r_i = d_i^\alpha$, the region of integration remains the same and integrand becomes $(2/\alpha)^2(r_1 r_2)^{2/\alpha-1}$ instead of $1$. Writing $\theta = \alpha/2$, we obtain
$$\frac{1}{2} a^{-1/\theta } \, _2F_1\left(\frac{1}{\theta },\frac{2}{\theta };\frac{\theta +2}{\theta };-b\right)$$
when $(a>0\land a<1\land a b+a\geq 1)$ or $a\geq 1$ and otherwise the result is
$$-a^{\frac{1}{\theta }} \left(\frac{1}{1-a b}\right)^{\frac{1}{\theta }}+\frac{1}{2} a^{\frac{1}{\theta }} (1-a b)^{-2/\theta } \, _2F_1\left(\frac{1}{\theta },\frac{2}{\theta };\frac{\theta +2}{\theta };1+\frac{1}{a b-1}\right)+1.$$
Here, $_2F_1$ is the hypergeometric function. The original case of $\alpha=2$ corresponds to $\theta=1$ and then these formulae reduce to the fourth and seventh of the eight previous cases. I have checked this result with a simulation, letting $\theta$ range from $1$ through $3$ and covering substantial ranges of $a$ and $b$.
I have never met your exercice before, but this one is interesting and leads to funny calculation.
First, let me say that your computation are not totally correct. In fact you have computed the law of the random variable M | Y=y. Here are a complete answer.
Let one determine the pdf of the random variable M | Y=y. I emphasise that here, Y is fixed. It is easy to see that
\begin{align}
\mathrm{p}(M \leq x | Y=y) =& p(U_1, \ldots,U_Y \leq x | Y=y) \\
=& \prod_{y=1}^y p(U_1 \leq x) \\
=& \prod_{y=1}^y \int_0^x 1 \mathrm{d}t \\
=& x^y \\
=& \int_0^x yt^{y-1} \mathrm{d}t
\end{align}
So the pdf of $M|Y=y$ is the function $t\rightarrow yt^{y-1}$.
To find the pdf of $M$, a simple chain rule is enough. Hence
\begin{align}
\mathrm{p}(M \leq x) =& \sum_{y=1}^{+\infty} \mathrm{p}(M\leq x, Y=y) \\
=& \sum_{y=1}^{+\infty} \mathrm{p}(M\leq x | Y=y) \mathrm{p}(Y=y) \\
=& \int_0^x \sum_{y=1}^{+\infty} \frac{1}{(e-1)y!} yt^{y-1} \mathrm{d}t \\
=& \int_0^x \sum_{y=1}^{+\infty} \frac{1}{(e-1) (y-1)!} t^{y-1} \mathrm{d}t \\
=& \int_0^x \frac{e^t}{e-1} \mathrm{d}t.
\end{align}
As a consequence, the pdf of $M$ is the function $t \rightarrow e^t / (e-1)$.
From the pdf of M, it is easy to derive the two first moments of the distribution. Two IPP lead to
\begin{align}
\mathbb{E}[M] =& \frac{1}{e-1} \\
\mathrm{Var}(M) =& \frac{e^2-3e+1}{(e-1)^2}
\end{align}
Here it is. I hope I have answered your question.
Best Answer
Some guide:
Notice that $$\mathbb{E}(A) \neq \pi (\mathbb{E}(R))^2$$ in general. Hence we can't say that since $\mathbb{E}(R)= \frac12$, hence $\mathbb{E}(A)=\frac{\pi}{2^2}$, this is not the right approach.
We have, $$\mathbb{E}(A) = \pi\mathbb{E}(R^2)$$ and to evaluate $\mathbb{E}(R^2)$, you might like use for the formula of $Var(X)=\mathbb{E}(X^2)-\mathbb{E}(X)^2$.
To find variance of $A$, you might want to use the same formula as well.