Background
We start by defining the following function $S(x)$:
$$ S(x) + S(x^2) + S(x^3) + \dots = x $$
Using the mobius inversion formula:
$$ S(x) = \sum_{r=1}^\infty \mu(r) x^r $$
We now define a nil-potent $(n+1) \times (n+1)$ matrix $\epsilon$ such that:
$$ \epsilon^{n+1} = \hat 0$$
Hence,
$$ S(\epsilon) = \sum_{r=1}^n \mu(r)\epsilon^r $$
and,
$$ S(\epsilon)+ S(\epsilon^2)+ \dots + S(\epsilon^n)= \epsilon $$
Let, us call the above equation the epsilon-n equation.
$$ \implies \frac{\partial}{\partial \epsilon} (S(\epsilon)+ S(\epsilon^2)+ \dots + S(\epsilon^n))= 1 $$
$$ \implies S'(\epsilon)+ 2 \epsilon S'(\epsilon^2)+ \dots + n \epsilon^{n-1 }S'(\epsilon^n)= 1 $$
We know $ (\epsilon^k)_{1,n+1} = 0 $ where $(\epsilon^k)_{1,n+1} $ is the n'th row and 1'st row and n'th column. Hence,
$$ \implies S'(\epsilon_{1,n+1})+ 2 \epsilon S'((\epsilon^2)_{1,n+1})+ \dots + n \epsilon^{n-1 }S'((\epsilon^n)_{1,n+1})= 1 $$
$$ \implies S'(0) = 1 $$
Now, taking the second derivative of the epsilon-n equation:
$$ \frac{\partial^2}{\partial \epsilon^2} (S(\epsilon)+ S(\epsilon^2)+ \dots + S(\epsilon^n))= 0 $$
$$ S''(\epsilon) + 2S'(\epsilon^2) + 4 \epsilon S''(\epsilon^2) + \dots = 0 $$
Again, comparing the 1'st row n'th column of each martix:
$$ S''(0) + 2 S'(0) = 0 $$
Repeating the above procedure $n$ times, we obtain a $n \times n$ matrix, $A$ :
\begin{equation*} \underbrace{\begin{pmatrix} 1&0&0&\dots&0\\ 2&1&0&\dots&0 \\ \vdots& & &\dots& \\ \end{pmatrix}}_{A} \underbrace{\begin{pmatrix} S'(0) \\ \vdots \\ S'^n (0) \\ \end{pmatrix}}_{B} = \underbrace{\begin{pmatrix} 1\\ \vdots\\ 0\\ \end{pmatrix}}_{C} \end{equation*}
Futher, we know, $S'^k(0) = \mu(k) k!$
\begin{equation*} \implies \underbrace{\begin{pmatrix} 1&0&0&\dots&0\\ 2&1&0&\dots&0 \\ \vdots& & &\dots& \\ \end{pmatrix}}_{A} \underbrace{\begin{pmatrix} \mu(1) \\ \vdots \\ \mu(n) n! \\ \end{pmatrix}}_{B} = \underbrace{\begin{pmatrix} 1\\ \vdots\\ 0\\ \end{pmatrix}}_{C} \end{equation*}
As the system of equations given by $AB=C$ is solvable, $A^{-1}$ must exist. Multiplying from the right a column matrix $D= \begin{pmatrix} 1/1! & 1/2! & 1/3! & \dots & 1/n!\end{pmatrix}$:
$$\implies BD = A^{-1} CD $$.
Taking the trace:
$$ M(n) = \sum_{r=1}^n \mu(r) = \text{Tr} (BD) = \text{Tr} (A^{-1} CD) $$
Questions
Is this correct? Is there any general formula for the general $ k \times k$ matrix of $A$ or $A^{-1}$ matrix? Is this method viable to calculate Merten's function? Does this method already exist?