Can we find the distribution and/or expected value of $$S=\sum_{n \geq 1} \frac{1}{(\text{lcm} (n,r_n))^2}$$
where $r_n$ is a uniformly distributed random integer, $r_n \in [1,n]$ and $\text{lcm}$ is the least common multiple?
Or maybe an estimate is possible?
Some boundaries are easy to see. For example:
$$1<S \leq \sum_{n \geq 1} \frac{1}{n^2}= \frac{\pi^2}{6} \approx 1.645$$
For two trials with $10000$ runs each and summing up to $n=20000$ I obtained the following distributions:
It seems as if there are two peaks.
For 20000 trials we obtain the following distribution:
For 100000 trials:
We can find an expression for the expected value of $S$ in terms of zeta function values.
We have $$ \mathbf{E}(S)=\sum_{n=1}^{\infty} \mathbf{E} \frac1{\mathrm{lcm}(n,r_n)^2}=\sum_{n=1}^{\infty} \frac1n\sum_{j=1}^n \frac1{\mathrm{lcm}(n,j)^2}. $$
We follow the proof of Theorem 7.1 in a paper by O. Bordelles here.
The sum over $j$ becomes
$$ \begin{align} \sum_{j=1}^n \frac1{\mathrm{lcm}(n,j)^2}&= \sum_{j=1}^n \frac{\mathrm{gcd}(n,j)^2}{n^2j^2}\\ &=\frac1{n^2}\sum_{d|n} d^2 \sum_{\substack{{j=1}\\{(j,n)=d}}}^n \frac1{j^2}=\frac1{n^2}\sum_{d|n} \sum_{\substack{{k\leq \frac nd}\\{(k,\frac nd)=1}}} \frac1{k^2}. \end{align} $$ In the second equality we use $\gcd(n,j)=d$. In the third equality, we use $j=dk$.
We use this into the following sum over $n$.
$$ \begin{align} \sum_{n\leq x} &\frac1n\sum_{j=1}^n \frac1{\mathrm{lcm}(n,j)^2} =\sum_{n\leq x} \frac1{n^3} \sum_{d|n} \sum_{\substack{{k\leq \frac nd}\\{(k,\frac nd)=1}}} \frac1{k^2}\\ &=\sum_{d\leq x} \frac1{d^3} \sum_{h\leq \frac xd} \frac1{h^3} \sum_{\substack{{k\leq h}\\{(k,h)=1}}} \frac1{k^2} \ (\textrm{ by }n=dh)\\ &=\sum_{d\leq x} \frac1{d^3}\sum_{h\leq \frac xd} \frac1{h^3}\sum_{\delta|h} \frac{\sum_{m\leq h/\delta} \mu(\delta)}{\delta^2 m^2} \ (\textrm{ by introducing }\delta|h, \delta|k, \mu(\delta), k=\delta m)\\ &=\sum_{d\leq x}\frac1{d^3}\sum_{\delta\leq \frac xd}\frac{\mu(\delta)}{\delta^5} \sum_{a\leq \frac x{d\delta}} \frac1{a^3}\sum_{m\leq a}\frac1{m^2} \ (\textrm{ by }h=\delta a). \end{align} $$
The inner double sum over $a$ and $m$, is a well-known Euler sum. The formula is given in Lemma 2.1 of this paper by Ce Xu, Yuhuan Yan, Zhijuan Shi. For the proof, see the citation [9] of this paper: here.
$$ C:=\sum_{a=1}^{\infty} \frac1{a^3}\sum_{m\leq a} \frac1{m^2} = 3\zeta(2)\zeta(3)-\frac{9\zeta(5)}2. $$
Then we have $$ \begin{align} \sum_{d\leq x}&\frac1{d^3}\sum_{\delta\leq \frac xd}\frac{\mu(\delta)}{\delta^5} \sum_{a\leq \frac x{d\delta}} \frac1{a^3}\sum_{m\leq a}\frac1{m^2} =\sum_{d\leq x}\frac1{d^3}\sum_{\delta\leq \frac xd} \frac{\mu(\delta)}{\delta^5} \left( C + O\left( \frac{d^2\delta^2}{x^2}\right)\right)\\ &=\sum_{d\leq x} \frac1{d^3}\left(\frac C{\zeta(5)} + O\left(\frac{d^4}{x^4}\right)\right)+O\left(\frac{\log x}{x^2}\right)\\ &=\frac{C\zeta(3)}{\zeta(5)} + O\left( \frac{\log x}{x^2}\right). \end{align} $$
Therefore, we obtain $$ \mathbf{E}(S)=\frac{\zeta(3)\left(3\zeta(2)\zeta(3)-\frac{9\zeta(5)}2\right)}{\zeta(5)}. $$
From this Wolfram Alpha input (here), the value of the above is approximately $$1.4673050041864341678704408312097934381417563912360385854969000128$$