Comenzando de la misma manera que Yanior Weg, asumimos $b_1$ y $b_2$ son fijos. Entonces $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$
En este gráfico de Desmos, se puede ver la representación de esto (mencionaré esto más tarde, por lo que puede ser útil tenerlo abierto). $P(r_1 \to b_1)$ es el área de la región sombreada dentro del cuadrado. Al limitar $x_2, y_2$ de manera tal que $x_1 < x_2 < 1$, $y_1 < y_2 < 1$, los cálculos posteriores pueden simplificarse. La respuesta final es entonces $$\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$$, donde $A$ es el área de la región sombreada dentro del cuadrado. $2A^2$ se multiplica por $4$ para tener en cuenta que $x_2 y $y_2.
La línea que separa la región sombreada tiene la ecuación $$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}$$
Esto se intersecta con $y = 1$ en $$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}$$
y con $y = 0$ en $$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}$$
Ahora hay cuatro casos: $(1)\ I_1 < 0, I_2 < 1, \ (2)\ I_1 < 0, I_2 > 1, \ (3)\ I_1 > 0, I_2 < 1$, y $(4)\ I_1 > 0, I_2 > 1$.
Tomando $$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$$ aquí $A_i$ representa el área de caso $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)$
Desde aquí, $$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}$$, donde $(1)$ es la región en $x_1, y_1, x_2, y_2$ tal que ocurre el caso 1 (restringido a $0 < x1 < x2 < 1$ y $0 < y1 < y2 < 1$), etc.
Para simplificar, hacer la sustitución $x_s = x_2 + x_1, x_d = x_2 - x_1$ y $y_s = y_2 + y_1, y_d = y_2 - y_1$ ayuda mucho. Luego, la integral necesitaría ser multiplicada por el Jacobiano de $\frac{1}{4}$.
Para el caso $1$, la integral se puede escribir como $$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)$$
Para el caso $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)$$
Para el caso $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)$$
Para el caso $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)$$
Sumando estos, obtenemos que $J = \frac{39}{800}+\frac{\pi}{120}-\frac{\ln(2)}{180}$. Multiplicando por $8$ da la respuesta final como $$\frac{39}{100}+\frac{\pi}{15}-\frac{2\ln(2)}{45} \approx 0.569$$, que está en el rango simulado que Arthur mencionó en los comentarios.
