[proofplan]
We prove the identity by comparing the coefficient of each formal monomial $n^{-s}$. First we make the coefficientwise meaning of the infinite Euler product precise by taking finite products over finite sets of primes. For a fixed integer $n$, the coefficient stabilizes once the finite set contains all prime divisors of $n$. Unique prime factorization identifies exactly one choice of local exponents with $n$, and multiplicativity identifies the corresponding coefficient with $f(n)$.
[/proofplan]
[step:Define the finite Euler products and their coefficients]
Let $\mathcal P$ denote the set of prime numbers. For a finite subset $P\subset \mathcal P$, define the finite formal Dirichlet series
\begin{align*}
E_P
:=
\prod_{p\in P}\left(\sum_{a=0}^{\infty}\frac{f(p^a)}{p^{as}}\right).
\end{align*}
The product is a formal product, so the coefficient of $n^{-s}$ in $E_P$ is obtained by choosing, for each $p\in P$, an exponent $a_p\in \mathbb N\cup\{0\}$ and multiplying the resulting local terms.
Thus a choice of exponents $(a_p)_{p\in P}$ contributes to the coefficient of $n^{-s}$ precisely when
\begin{align*}
n=\prod_{p\in P}p^{a_p}.
\end{align*}
For such a choice, the contributed coefficient is
\begin{align*}
\prod_{p\in P} f(p^{a_p}).
\end{align*}
[guided]
We first replace the infinite product by finite products, because formal infinite products are understood coefficient by coefficient. Let $\mathcal P$ be the set of all prime numbers, and let $P\subset \mathcal P$ be finite. We define
\begin{align*}
E_P
:=
\prod_{p\in P}\left(\sum_{a=0}^{\infty}\frac{f(p^a)}{p^{as}}\right).
\end{align*}
Expanding this finite product means making one choice from each factor. For each prime $p\in P$, choose an exponent $a_p\in \mathbb N\cup\{0\}$, contributing the local term
\begin{align*}
\frac{f(p^{a_p})}{p^{a_p s}}.
\end{align*}
Multiplying these selected local terms gives
\begin{align*}
\prod_{p\in P}\frac{f(p^{a_p})}{p^{a_p s}}
=
\frac{\prod_{p\in P}f(p^{a_p})}{\left(\prod_{p\in P}p^{a_p}\right)^s}.
\end{align*}
Therefore this selected term contributes to the coefficient of $n^{-s}$ exactly when
\begin{align*}
n=\prod_{p\in P}p^{a_p}.
\end{align*}
In that case, its coefficient is
\begin{align*}
\prod_{p\in P}f(p^{a_p}).
\end{align*}
[/guided]
[/step]
[step:Show that each coefficient stabilizes in the directed finite products]
Fix $n\in\mathbb N$. Let $P(n)\subset \mathcal P$ denote the finite set of primes dividing $n$. If $P$ is a finite subset of $\mathcal P$ with $P(n)\subset P$, then the coefficient of $n^{-s}$ in $E_P$ is independent of $P$.
Indeed, by unique prime factorization, there is a unique family of exponents $(\nu_p(n))_{p\in P}$ with $\nu_p(n)\in \mathbb N\cup\{0\}$ such that
\begin{align*}
n=\prod_{p\in P}p^{\nu_p(n)},
\end{align*}
where $\nu_p(n)=0$ for $p\notin P(n)$. Hence the coefficient of $n^{-s}$ in every such $E_P$ is
\begin{align*}
\prod_{p\in P}f(p^{\nu_p(n)}).
\end{align*}
For $p\notin P(n)$, one has $\nu_p(n)=0$, so $p^{\nu_p(n)}=1$ and $f(p^{\nu_p(n)})=f(1)=1$. Therefore
\begin{align*}
\prod_{p\in P}f(p^{\nu_p(n)})
=
\prod_{p\in P(n)}f(p^{\nu_p(n)}),
\end{align*}
which depends only on $n$.
[guided]
Now fix a particular integer $n\in\mathbb N$. The infinite product is meaningful if the coefficient of $n^{-s}$ eventually stops changing as we enlarge the finite set of primes. Let $P(n)\subset\mathcal P$ be the finite set of primes that divide $n$.
Take any finite set $P\subset\mathcal P$ with $P(n)\subset P$. By unique prime factorization, there is exactly one exponent $\nu_p(n)\in\mathbb N\cup\{0\}$ for each $p\in P$ such that
\begin{align*}
n=\prod_{p\in P}p^{\nu_p(n)}.
\end{align*}
For primes $p\in P$ that do not divide $n$, the exponent is $\nu_p(n)=0$.
From the coefficient computation in the previous step, the coefficient of $n^{-s}$ in $E_P$ is therefore
\begin{align*}
\prod_{p\in P}f(p^{\nu_p(n)}).
\end{align*}
The extra primes in $P\setminus P(n)$ do not change this product, because for each such prime $p$ we have $\nu_p(n)=0$, hence
\begin{align*}
f(p^{\nu_p(n)})=f(p^0)=f(1)=1.
\end{align*}
Thus
\begin{align*}
\prod_{p\in P}f(p^{\nu_p(n)})
=
\prod_{p\in P(n)}f(p^{\nu_p(n)}).
\end{align*}
So the coefficient of $n^{-s}$ stabilizes once all prime divisors of $n$ have been included.
[/guided]
[/step]
[step:Use multiplicativity to identify the stabilized coefficient with $f(n)$]
Write the prime factorization of $n$ as
\begin{align*}
n=\prod_{p\in P(n)}p^{\nu_p(n)}.
\end{align*}
The factors $p^{\nu_p(n)}$ are pairwise coprime as $p$ ranges over $P(n)$. Since $f$ is multiplicative, repeated application of the defining identity gives
\begin{align*}
f(n)
=
f\left(\prod_{p\in P(n)}p^{\nu_p(n)}\right)
=
\prod_{p\in P(n)}f(p^{\nu_p(n)}).
\end{align*}
Therefore the stabilized coefficient of $n^{-s}$ in the Euler product is exactly $f(n)$.
[guided]
We have reduced the problem to identifying the stabilized coefficient
\begin{align*}
\prod_{p\in P(n)}f(p^{\nu_p(n)}).
\end{align*}
The prime factorization of $n$ is
\begin{align*}
n=\prod_{p\in P(n)}p^{\nu_p(n)}.
\end{align*}
If $p$ and $q$ are distinct primes in $P(n)$, then $p^{\nu_p(n)}$ and $q^{\nu_q(n)}$ are coprime. Therefore all the factors in this product are pairwise coprime.
The multiplicativity hypothesis says that $f(xy)=f(x)f(y)$ whenever $\gcd(x,y)=1$. Applying this repeatedly to the pairwise coprime prime-power factors gives
\begin{align*}
f(n)
=
f\left(\prod_{p\in P(n)}p^{\nu_p(n)}\right)
=
\prod_{p\in P(n)}f(p^{\nu_p(n)}).
\end{align*}
This is exactly the coefficient obtained from the formal Euler product after stabilization.
[/guided]
[/step]
[step:Conclude equality of the formal Dirichlet series]
For every $n\in\mathbb N$, the coefficient of $n^{-s}$ in the coefficientwise infinite product
\begin{align*}
\prod_{p}\left(\sum_{a=0}^{\infty}\frac{f(p^a)}{p^{as}}\right)
\end{align*}
is $f(n)$. The coefficient of $n^{-s}$ in
\begin{align*}
\sum_{m=1}^{\infty}\frac{f(m)}{m^s}
\end{align*}
is also $f(n)$. Since formal Dirichlet series are equal exactly when all their coefficients are equal, the desired identity follows:
\begin{align*}
\sum_{n=1}^{\infty}\frac{f(n)}{n^s}
=
\prod_{p}\left(\sum_{a=0}^{\infty}\frac{f(p^a)}{p^{as}}\right).
\end{align*}
[/step]
