Starting off in the same manner as Yanior Weg, assume $b_1$ and $b_2$ are fixed. Then $P(r_2 \to b_1|r_1 \to b_1) = \frac{P(r_2 \to b_1 \bigcap r_1 \to b_1)}{P(r_1 \to b_1)} = \frac{P(r_2 \to b_1 \bigcap r_1 \to b_1)}{\frac{1}{2}} = 2 \cdot P(r_2 \to b_1 \bigcap r_1 \to b_1) = 2(P(r_1 \to b_1))^2$
In this Desmos plot, the depiction of this can be seen (I'll be referring to this later, so it may be useful to have it open). $P(r_1 \to b_1)$ is the area of the shaded region inside the square. By limiting $x_2, y_2$ such that $x_1 < x_2 < 1$, $y_1 < y_2 < 1$, later calculations can be simplified. The final answer is then $$\int_0^1\int_0^1\int_{x_1}^1\int_{y_1}^14\cdot2A^2 dy_2dx_2dy_1dx_1 = 8\int_0^1\int_0^1\int_{x_1}^1\int_{y_1}^1A^2 dy_2dx_2dy_1dx_1$$, where $A$ is the area of the shaded region inside the square. $2A^2$ is multiplied by $4$ to account for $x_2<x_1$ and $y_2<y_1$.
The line separating the shaded region has equation $$y = f\left(x\right)=-\frac{x_{1}-x_{2}}{y_{1}-y_{2}}x+\frac{\left(x_{1}-x_{2}\right)\left(x_{1}+x_{2}\right)}{2\left(y_{1}-y_{2}\right)}+\frac{y_{1}+y_{2}}{2}$$
This intersects $y = 1$ at $$x = I_1 = \frac{\left(y_{1}-y_{2}\right)\left(y_{1}+y_{2}-2\right)}{2\left(x_{1}-x_{2}\right)}+\frac{x_{1}+x_{2}}{2}$$
and $y = 0$ at $$x = I_2 = \frac{\left(y_{1}-y_{2}\right)\left(y_{2}+y_{1}\right)}{2\left(x_{1}-x_{2}\right)}+\frac{x_{1}+x_{2}}{2}$$
There are now four cases: $(1)\ I_1 < 0, I_2 < 1, \ (2)\ I_1 < 0, I_2 > 1, \ (3)\ I_1 > 0, I_2 < 1$, and $(4)\ I_1 > 0, I_2 > 1$.
Letting $$F(x) = \int_0^x f(t)dt = -\frac{x_{1}-x_{2}}{2\left(y_{1}-y_{2}\right)}x^{2}+\left(\frac{\left(x_{1}-x_{2}\right)\left(x_{1}+x_{2}\right)}{2\left(y_{1}-y_{2}\right)}+\frac{y_{1}+y_{2}}{2}\right)x$$
here $A_i$ represents the area of case $i$:
$A_1 = F(I_2)$
$A_2 = F(1)$
$A_3 = I_1 - F(I_1) + F(I_2)$
$A_4 = I_1 - F(I_1) + F(1)$
From here, $$J = \underbrace{\int_{(1)}A_1^2dy_2dx_2dy_1dx_1}_{J_1}+\underbrace{\int_{(2)}A_2^2dy_2dx_2dy_1dx_1}_{J_2}+\underbrace{\int_{(3)}A_3^2dy_2dx_2dy_1dx_1}_{J_3}+\underbrace{\int_{(4)}A_4^2dy_2dx_2dy_1dx_1}_{J_4}$$, where $(1)$ is the region in $x_1, y_1, x_2, y_2$ such that case 1 happens (restricted to $0 < x1 < x2 < 1$ and $0 < y1 < y2 < 1$), etc.
To simplify, making the substitution $x_s = x_2 + x_1, x_d = x_2 - x_1$ and $y_s = y_2 + y_1, y_d = y_2 - y_1$ helps a lot. The integral would then need to be multiplied by the Jacobian of $\frac{1}{4}$.
For case $1$, the integral can be written out as $$J_1 = \frac{1}{4}\int_0^1 \int_{x_d}^{2-x_d} \left(\int_{1-\sqrt{1-x_{d}x_{s}}}^{x_{d}}\int_{0}^{-\frac{x_{d}x_{s}}{y_{d}}+2}A_1^2dy_{s}dy_{d}+\int_{x_{d}}^{\sqrt{2x_{d}-x_{d}x_{s}}}\int_{0}^{\frac{2x_{d}-x_{d}x_{s}}{y_{d}}}A_1^{2}dy_{s}dy_{d}-\int_{1-\sqrt{1-x_{d}x_{s}}}^{\sqrt{2x_{d}-x_{d}x_{s}}}\int_{0}^{y_{d}}A_1^{2}dy_{s}dy_{d}\right)dx_s dx_d = \frac{1}{32}-\frac{1}{24}\ln(2)$$
For case $2$: $$J_2 = \frac{1}{4}\int_0^1 \int_{x_d}^{2-x_d} \left(\int_{x_{d}}^{\sqrt{x_{s}x_{d}}}\int_{0}^{2-\frac{x_{s}x_{d}}{y_{d}}}A_2^{2}dy_{s}dy_{d}+\int_{\sqrt{x_{s}x_{d}}}^{1}\int_{0}^{2-y_{d}}A_2^{2}dy_{s}dy_{d}-\int_{x_{d}}^{\sqrt{2x_{d}-x_{d}x_{s}}}\int_{0}^{\frac{2x_{d}-x_{d}x_{s}}{y_{d}}}A_2^{2}dy_{s}dy_{d}-\int_{\sqrt{2x_{d}-x_{d}x_{s}}}^{1}\int_{0}^{y_{d}}A_2^{2}dy_{s}dy_{d}\right) dx_s dx_d = \frac{2501}{14400}-\frac{3}{80}\pi-\frac{2}{45}\ln\left(2\right)$$
For case $3$: $$J_3 = \frac{1}{4}\int_0^1 \int_{y_d}^{2-y_d} \left(\int_{y_{d}}^{\sqrt{y_{d}y_{s}}}\int_{0}^{2-\frac{y_{d}y_{s}}{x_{d}}}A_{3}^{2}dx_{s}dx_{d}+\int_{\sqrt{y_{d}y_{s}}}^{1}\int_{0}^{2-x_{d}}A_{3}^{2}dx_{s}dx_{d}-\int_{y_{d}}^{\sqrt{y_{d}\left(2-y_{s}\right)}}\int_{0}^{\frac{y_{d}\left(2-y_{s}\right)}{x_{d}}}A_{3}^{2}dx_{s}dx_{d}-\int_{\sqrt{y_{d}\left(2-y_{s}\right)}}^{1}\int_{0}^{x_{d}}A_{3}^{2}dx_{s}dx_{d}\right) dy_s dy_d = \frac{2501}{14400}-\frac{3}{80}\pi-\frac{2}{45}\ln\left(2\right)$$
For case $4$: $$J_4 = \frac{1}{4}\int_{0}^{1}\int_{1-\sqrt{1-x_{d}^{2}}}^{x_{d}}\int_{2+\frac{y_{d}^{2}-2y_{d}}{x_{d}}}^{2-x_{d}}\int_{\frac{2x_{d}-x_{d}x_{s}}{y_{d}}}^{2-y_{d}}A_{4}^{2}dy_{s}dx_{s}dy_{d}dx_{d}+\frac{1}{4}\int_{0}^{1}\int_{x_{d}}^{\sqrt{2x_{d}-x_{d}^{2}}}\int_{\frac{y_{d}^{2}}{x_{d}}}^{2-x_{d}}\int_{2-\frac{x_{d}x_{s}}{y_{d}}}^{2-y_{d}}A_{4}^{2}dy_{s}dx_{s}dy_{d}dx_{d} = -\frac{95}{288}+\frac{1}{12}\pi+\frac{1}{8}\ln(2)$$
Adding these up yields that $J = \frac{39}{800}+\frac{\pi}{120}-\frac{\ln(2)}{180}$. Multiplying by $8$ gives the final answer as $$\frac{39}{100}+\frac{\pi}{15}-\frac{2\ln(2)}{45} \approx 0.569$$, which is in the simulated range Arthur mentioned in the comments.