Fix integers $r,s\geq1$, not both 1.
Definition. We say that integers $a_{1},\ldots,a_{s}$ are relatively $r$-prime if their greatest common divisor has no perfect $r$th power factors > 1.
(When $r=1$, we say that the integers are relatively prime. When $s=1$, we say that the integer is $r$-free.)
Let $A_{n}\left(d\right)$ be the set of integers in $I_{n}=\left\{0,\ldots,n\right\}$ which are divisible by $d$ and let $S^{r,s}_{n}$ be the set of $s$-tuples of integers in $I_{n}$ which are relatively $r$-prime. Then, the following result holds.
Proposition. Let $P_{n}$ be any probability distribution on $I_{n}$. Then, for $n>1$, we have $$P_{n}^{s}\left(S^{r,s}_{n}\right)=\sum_{1\leq d\leq n^{1/r}}\mu\left(d\right)\left(P_{n}\left(A_{n}\left(d^{r}\right)\right)^{s}-P_{n}\left(\left\{0\right\}\right)^{s}\right),$$ where $\mu$ is the Möbius function.
For the uniform and binomial distributions, it is known that $$\lim_{n\to\infty}P_{n}^{s}\left(S^{r,s}_{n}\right)=\frac{1}{\zeta\left(rs\right)}.$$ As a curiosity, I want to evaluate the probability on other distributions. For now, I want to use the proposition on the following distribution (ideal soliton distribution):
Define $P_{n}$ on $I_{n}$ by $$P_{n}\left(\left\{k\right\}\right)=\begin{cases} 0, & k=0\\ 1/n, & k=1\\ \frac{1}{k\left(k-1\right)}, & \mathrm{otherwise} \end{cases}.$$
My analysis so far:
Since $P_{n}\left(\left\{0\right\}\right) = 0$, we may omit that part from the sum, leaving $$P_{n}^{s}\left(S^{r,s}_{n}\right)=\sum_{1\leq d\leq n^{1/r}}\mu\left(d\right)\left(P_{n}\left(A_{n}\left(d^{r}\right)\right)^{s}\right).$$
Also, since $A_{n}\left(1\right)=I_{n}$, we may pull that term out as 1, leaving us with $$P_{n}^{s}\left(S^{r,s}_{n}\right)=1+\sum_{2\leq d\leq n^{1/r}}\mu\left(d\right)\left(P_{n}\left(A_{n}\left(d^{r}\right)\right)^{s}\right).$$
Since, for $d>1$, $1\not\in A_{n}\left(d^{r}\right)$, we have $$P_{n}\left(A_{n}\left(d^{r}\right)\right)=\sum_{\substack{2\leq k\leq n \\ k\equiv0\,\left(d^{r}\right)}}\frac{1}{k\left(k-1\right)}=\sum_{k=1}^{n/d^{r}}\frac{1}{d^{r}k\left(d^{r}k-1\right)}.$$
According to Wolfram Alpha, t (See Edit) Taking $n\to\infty$, we get (*) $$\lim_{n\to\infty}P_{n}\left(A_{n}\left(d^{r}\right)\right)=\sum_{k=1}^{\infty}\frac{1}{d^{r}k\left(d^{r}k-1\right)}=-d^{-r}\left(\psi\left(1-d^{-r}\right)+\gamma\right),$$ where $\psi$ is the digamma function and $\gamma$ is the Euler-Mascheroni constant.
From this, we get $$\lim_{n\to\infty}P_{n}^{s}\left(S^{r,s}_{n}\right)=1+\sum_{d=2}^{\infty}\mu\left(d\right)\left(-d^{-r}\left(\psi\left(1-d^{-r}\right)+\gamma\right)\right)^{s}=$$
$$1+\left(-1\right)^{s}\sum_{d=2}^{\infty}d^{-rs}\mu\left(d\right)\left(\psi\left(1-d^{-r}\right)+\gamma\right)^{s}.$$
With some manipulation, we have $$\sum_{d=2}^{\infty}d^{-rs}\mu\left(d\right)\left(\psi\left(1-d^{-r}\right)+\gamma\right)^{s}=$$
$$\gamma^{s}\sum_{d=2}^{\infty}d^{-rs}\mu\left(d\right) + \sum_{j=1}^{s}\binom{s}{j}\gamma^{s-j}\sum_{d=2}^{\infty}\psi\left(1-d^{-r}\right)^{j}d^{-rs}\mu\left(d\right)=$$
$$\gamma^{s}\left(\frac{1}{\zeta\left(rs\right)}-1\right)+\sum_{j=1}^{s}\binom{s}{j}\gamma^{s-j}\sum_{d=2}^{\infty}\psi\left(1-d^{-r}\right)^{j}d^{-rs}\mu\left(d\right).$$
Therefore, $$\lim_{n\to\infty}P_{n}^{s}\left(S^{r,s}_{n}\right) = 1 + \left(-1\right)^{s}\left(\gamma^{s}\left(\frac{1}{\zeta\left(rs\right)}-1\right)+\sum_{j=1}^{s}\binom{s}{j}\gamma^{s-j}\sum_{d=2}^{\infty}\psi\left(1-d^{-r}\right)^{j}d^{-rs}\mu\left(d\right)\right).$$
My problem is getting some workable form for $$\sum_{d=2}^{\infty}\psi\left(1-d^{-r}\right)^{j}d^{-rs}\mu\left(d\right)$$ and using that to evaluate $$\sum_{j=1}^{s}\binom{s}{j}\gamma^{s-j}\sum_{d=2}^{\infty}\psi\left(1-d^{-r}\right)^{j}d^{-rs}\mu\left(d\right)$$ in terms of $r$ and $s$.
In case it's easier to evaluate, here is the probability if we don't pull the $\gamma^{s}\left(\frac{1}{\zeta\left(rs\right)} - 1\right)$ term out in the manipulation: $$\lim_{n\to\infty}P_{n}^{s}\left(S^{r,s}_{n}\right)=1+\left(-1\right)^{s}\sum_{j=0}^{s}\binom{s}{j}\gamma^{s-j}\sum_{d=2}^{\infty}\psi\left(1-d^{-r}\right)^{j}d^{-rs}\mu\left(d\right).$$
Two questions:
How would I prove the Wolfram Alpha result (*) above, that$$\sum_{k=1}^{\infty}\frac{1}{d^{r}k\left(d^{r}k-1\right)}=-d^{-r}\left(\psi\left(1-d^{-r}\right)+\gamma\right)?$$
SOLVED Q1: See edit.
- How would I even begin to analyze the sums that I ended up getting?
It's been a couple of years since I've really worked on a mathematics problem beyond basic calculus / linear algebra, so I'm a bit rusty. If there are some references in the literature to sums of these forms, or similar, they would be very much appreciated. Thank you.
EDIT: I've proven the Wolfram Alpha identity (*) (by assuming another series formula without proof).
On the Wikipedia page for the digamma function, there's a section with a series formula that fits, namely $$\psi\left(z+1\right)=-\gamma + \sum_{k=1}^{\infty}\frac{z}{k\left(k+z\right)}$$ for $z\neq-1,-2,-3\ldots.$ Setting $z=-d^{-r}$ and adding $\gamma$ to both sides, and putting $d^{-r}$ into the denominator as $d^{r}$, we get $$\psi\left(1-d^{-r}\right)+\gamma=-\sum_{k=1}^{\infty}\frac{d^{-r}}{k\left(k-d^{-r}\right)}=-\sum_{k=1}^{\infty}\frac{1}{k\left(d^{r}k-1\right)}.$$
Finally, multiplying by $-d^{-r}$ on both sides gives $$-d^{-r}\left(\psi\left(1-d^{-r}\right)+\gamma\right)=\sum_{k=1}^{\infty}\frac{1}{d^{r}k\left(d^{r}k-1\right)},$$ as required.
I suppose, since the source that Wikipedia provides is just a handbook of formulas (6.3.16) listed without proof, I could amend question 1 to ask for a proof of the series formula itself, but I'll postpone that until I've taken the time to try it myself.