Se puede abordar esta cuestión formando una cuadrícula desplazada de cuadrados de radios $1/k$ de forma aproximadamente triangular que encaje en el sector superior izquierdo del rectángulo, con lo que los cuadrados se pueden rellenar con círculos y ya hemos terminado.
Para formalizar esta idea, empecemos por darnos suficiente margen para necesitar sólo el sector superior izquierdo. El cálculo del post señala que no necesitamos gran parte del área restante si rellenamos el primer $n = 20$ pero para tener más espacio (y puesto que lo necesitaremos más adelante de todos modos), rellenemos el primer campo $65$ bolas a mano para que podamos empezar con $n=66$ :
(Gráfico Desmos aquí modificado a partir del gráfico de la OP)
Ahora que tenemos esto, nuestro objetivo será construir un triángulo en el sector superior izquierdo a partir de cajas, con radios como tales:
Por supuesto, está la cuestión de cómo hacer que los cuadrados sean tangentes cuando los radios son decrecientes, así que estipulamos que en cada etapa primero ponemos una caja lo más a la izquierda y luego lo más arriba posible, y luego ponemos las cajas restantes a ras de la parte superior derecha de la caja asociada de la etapa anterior. Gráficamente, donde una caja está etiquetada $m$ si se añade en el $m$ paso, esto parece:
Si trazamos una línea con ángulo $45 \deg$ en la parte inferior derecha de la caja más inferior, obtenemos una caja delimitadora de un triángulo:
La longitud del lado corto de este triángulo isósceles se obtiene sumando primero los diámetros de la columna de casillas situadas más a la izquierda y añadiendo después el diámetro de la casilla situada más abajo a la izquierda. En concreto, si realizamos los números triangulares $T_n = \frac{n(n+1)}{2}$ podemos acotar la longitud lateral $S_m$ en el $m$ ª fase de la construcción como: $$S_m = \frac{2}{n + \frac{m(m-1)}{2}} + \sum_{k=0}^{m-1} \frac{2}{n + \frac{k(k+1)}{2}} \le \frac{2}{n} + \sum_{k=0}^{\infty} \frac{2}{n + \frac{k(k+1)}{2}}$$
De hecho, aquí utilizamos el hecho $n = 66$ para determinar:
$$S_m \le \frac{2}{66} + \sum_{k=0}^{\infty} \frac{2}{66 + \frac{k(k+1)}{2}} = \frac{2}{66} + \frac{4\pi \tanh \left(\frac{\sqrt{527} \pi}{2} \right)}{\sqrt{527}} = 0.577703\ldots$$
Ahora bien, esto es importante, porque resulta ser un límite suficientemente bueno. En efecto, si se traza la tangente al círculo unitario en el diagrama de la OP en el punto $(\cos(3\pi/4), \sin(3\pi/4))$ y ve dónde se cruza con el límite vertical de la izquierda, un cálculo rápido (utilizando el hecho de que la línea tiene ecuación $y = x + \sqrt{2}$ ) muestra el mayor triángulo isósceles posible que puede encajar tiene lado corto de longitud $2 - \sqrt{2} = 0.585786\ldots$ por lo que nuestro triángulo construido arriba encaja.
Como comprobación de cordura, ahora podemos poner los círculos restantes de radio $1/k$ en las cajas de radio $1/k$ y determinar el área resultante. En efecto, el área de las cajas puede calcularse como $\sum_{k = 66}^\infty \left(\frac{2}{k}\right)^2 \approx 0.061$ mientras que el área del sector superior izquierdo viene dada por $1 - \frac{\pi}{4} \approx 0.21$ así que nos sobra mucha superficie.