We first show that $f(f(m)) = m$ for all $m$.
We start by substituting $n = 1$. This yields $f(f(m)f(1)) = m$. Then on the one hand $f$ is injective, since if $f(m) = f(n)$ then it immediately follows that $m = f(f(m)f(1)) = f(f(n)f(1)) = n$. On the other hand, $f$ is surjective, since for every $m$ there is an $a$ such that $f(a) = m$, namely $a = f(m)f(1)$. Hence $f$ is bijective. Now we substitute $m = n$, which yields $f(f(m)^2) = m^2 = f(f(m^2)f(1))$. Because of the injectivity of $f$, it follows
$$f(m)^2 = f(m^2)f(1).$$
Now choose $m_1$ for which $f(m_1) = 1$. (It follows from our proof of the surjectivity of $f$ that we can choose $m_1 = f(1)^2$ for this, but that is not important.) If we substitute $m = m_1$ in the above we obtain $1 = f(m_1^2)f(1)$. In particular, $f(1)$ is a divisor of $1$, and since it must be a natural number we conclude that $f(1) = 1$.
Substituting $n = 1$, now yields $f(f(m)) = m$.
Now we can substitute $f(m)$ and $f(n)$ for $m$ and $n$ respectively, so that we obtain the more classical equation
$$f(mn) = f(f(f(m))f(f(n))) = f(m)f(n).$$
It follows that $f$ is determined by the values $f(p)$ with $p$ prime.
Let $p$ be a prime number, and let $m$ and $n$ be natural numbers such that $f(p) = mn$. Then it follows that $f(m)f(n) = f(mn) = f(f(p)) = p$. Since $f(m)$ and $f(n)$ are natural numbers, it must hold that $f(m) = 1$ or $f(n) = 1$. Because of injectivity, the only $m$ for which $f(m) = 1$ is the number $m = 1$. So either $m = 1$ or $n = 1$. We conclude that $f(p)$ is a prime number.
Now suppose that $p \nmid 2024$ holds. In particular, we have $\gcd(p, 2024) = 1$, so $p$ has a multiplicative inverse modulo $2024$. Let $b$ be a natural number in the residue class of this inverse, i.e. $bp \equiv 1 \pmod{2024}$. Since $bp \ge p > 1$, it follows from the second equation that $f(bp) = bp$, but also that $f(bp+2024p) = bp+2024p$. However, it follows from the previous paragraph that $f(bp) = f(b)f(p)$ and $f((b+2024)p) = f(b+2024)f(p)$. Therefore $f(p)$ is a divisor of both $bp$ and $bp+2024p$. Now we calculate
$$\gcd(bp, bp + 2024p) = \gcd(bp, 2024p) = \gcd(b, 2024)p = p,$$
since $b$ is the multiplicative inverse of $p$ modulo $2024$. Because of injectivity, we know that $f(p) \ne 1$, so $f(p) \mid p$ implies that $f(p) = p$.
Suppose for a prime divisor $p$ of $2024$ that $f(p) = q$, with $q \ne p$ a different prime number. Then it follows that $f(q) = f(f(p)) = p$. So all prime numbers not mapped to themselves by $f$ form pairs $(p, q)$ with $f(p) = q$ and $f(q) = p$. Since $2024 = 2^3 \cdot 11 \cdot 23$, we have four possible cases:
* $f(p) = p$ for all primes $p$. So $f$ is the identity,
* $f(2) = 11$, $f(11) = 2$ and $f(p) = p$ for all $p \nmid 2, 11$,
* $f(2) = 23$, $f(23) = 2$ and $f(p) = p$ for all $p \nmid 2, 23$,
* $f(11) = 23$, $f(23) = 11$ and $f(p) = p$ for all $p \nmid 11, 23$.
It is easy to check that all four of these cases actually lead to a function that satisfies the given conditions. $\square$