Why does Quantum Phase Estimation first 'naturally' convert the phase into its Fourier basis states?

When you apply QFT on an state $|x \rangle$ the output can be shown as a tensor product form :
$$\text{QFT}(|x\rangle) = \frac{1}{\sqrt{N}}\bigotimes_{k=1}^n\bigg(|0\rangle+e^{\frac{2\pi ix}{2^k}}|1\rangle\bigg)\,, \tag{1}$$ where $$N=2^n\,.$$

Let $x$ be a binary string $x_1x_2..x_{n-1}x_n$, then (1) can be written as : $ \begin{align} \text{QFT}(|x_1x_2..x_{n-1}x_n\rangle) = \frac{1}{\sqrt{N}}&\bigg(|0\rangle + e^{\frac{2 \pi i x_1x_2..x_{n-1}x_n}{2^1}}|1\rangle\bigg)\\ &\otimes\bigg(|0\rangle + e^{\frac{2 \pi i x_1x_2..x_{n-1}x_n}{2^2}}|1\rangle\bigg)\\&\otimes\cdots\otimes\bigg(|0\rangle + e^{\frac{2 \pi i x_1x_2..x_{n-1}x_n}{2^n}}|1\rangle\bigg)\,. \tag{2} \end{align}$ Integer part of $\frac{2\pi i x}{2^k}$ in $e^{\frac{2\pi i x}{2^k}}$ can be ignored as it's 1. Thus, we just consider its fractional part, (2) can be written as :
$\begin{align} = \frac{1}{\sqrt{N}}&\bigg(|0\rangle + e^{2 \pi i 0.x_n}|1\rangle\bigg)\\&\otimes\bigg(|0\rangle + e^{2 \pi i 0.x_{n-1}x_n}|1\rangle\bigg)\\&\otimes\cdots\otimes\bigg(|0\rangle + e^{2 \pi i 0.x_1x_2...x_{n-1}x_n}|1\rangle\bigg)\,. \tag{3} \end{align}$

Next, let's see Fig.1 showing a diagram of phase estimation:

$|\psi_1 \rangle$ and $|\psi_2 \rangle$ can be written as follows :
$|\psi_1\rangle = H^{\otimes m}|0\rangle^{\otimes m}\otimes|\psi\rangle=\frac{1}{\sqrt{2^m}}(|0\rangle+|1\rangle)\otimes \cdots \otimes(|0\rangle+|1\rangle)\otimes|\psi\rangle\tag{4}$
When you apply controlled-U gate on $|j\rangle\otimes|\psi\rangle$ where $|j\rangle$ is control qubits it evolves to $|j\rangle U^j|\psi\rangle$. Hence $|\psi_2\rangle$ is represented as follows:

$|\psi_2\rangle =\frac{1}{\sqrt{2^m}}(|0\rangle+|1\rangle U^{2^{m-1}})\otimes(|0\rangle+|1\rangle U^{2^{m-2}})\otimes \cdots \otimes(|0\rangle+|1\rangle U^{2^{0}})\otimes|\psi\rangle\tag{5}$
Since $U^{2^{m-1}}=e^{2\pi i \theta 2^{m-1}}$ (5) can be written as (6):
$|\psi_2\rangle =\frac{1}{\sqrt{2^m}}(|0\rangle+e^{2\pi i \theta 2^{m-1}}|1\rangle )\otimes(|0\rangle+e^{2\pi i \theta 2^{m-2}}|1\rangle )\otimes \cdots \otimes(|0\rangle+e^{2\pi i \theta 2^{0}}|1\rangle )\otimes|\psi\rangle\tag{6}$
Therefore, first register can be represented as follows :
$\frac{1}{\sqrt{2^m}}(|0\rangle + e^{2\pi i \theta 2^{m-1}}|1\rangle ) \otimes (|0\rangle + e^{2\pi i \theta 2^{m-2}}|1\rangle ) \otimes \cdots \otimes (|0\rangle + e^{2\pi i \theta 2^{0}}|1\rangle ) \tag{7}$

A phase $\theta$ is an m-qubit real number. Let $\theta$ be :
$\theta = 0.\theta_1\theta_2...\theta_m \tag{8}$
Putting (8) into (7) gives :
$\frac{1}{\sqrt{2^m}}(|0\rangle + e^{2\pi i 0.\theta_m} |1\rangle ) \otimes (|0\rangle + e^{2\pi i 0.\theta_{m-1}\theta_{m}}|1\rangle ) \otimes \cdots \otimes (|0\rangle + e^{2\pi i 0.\theta_1\theta_2..\theta_{m-1}\theta_{m} }|1\rangle ) \tag{9}$

Comparing (3) and (9), you see they have the same form. Hence, by applying an Inverse QFT on (9) we can get (8).

-EDIT-

QFT can be expressed in a tensor product form. I've added its derivation below.

Let input $x$ be a binary string $x=[x_1x_2..x_n] \, x_k \in \{0, 1\} \, k=\{1,2,..,n\}$

QFT can be expressed in a form (10): $$\text{QFT}(|x\rangle) = \frac{1}{\sqrt{N}}\sum_{j=0}^{N-1}\omega_N^{jx}|j\rangle_n \tag{10}$$

Let $j$ also be a binary string $j = [j_1j_2..j_n] \, j_k \in \{0,1\}, k = \{1,2,..,n\}$
Since $\sum_{j=0}^{N-1}$ is $\sum_{j_1=0}^{1}\sum_{j_2=0}^{1}\cdots \sum_{j_n=0}^{1}$, (10) can be written as (11):
$$\frac{1}{\sqrt{N}}\sum_{j_1=0}^{1}\sum_{j_2=0}^{1}\cdots \sum_{j_n=0}^{1}\omega_N^{jx}|j\rangle_n \tag{11}$$
Rewriting j as $j_1j_2...j_n$ in (11) gives (12): $$\frac{1}{\sqrt{N}}\sum_{j_1=0}^{1}\sum_{j_2=0}^{1}\cdots \sum_{j_n=0}^{1}\omega_N^{jx}|j_1j_2..j_n\rangle \tag{12}$$

Next, the exponent of $\omega$ in (12) can be expressed using $\sum$ form. As j can also be expressed using the sum of power of 2 times $j_k$, namely
$j = \sum_{k=1}^{n}2^{n-k} \cdot j_k \tag{13}$

Putting (13) to (12) gives (14):
$$\frac{1}{\sqrt{N}}\sum_{j_1=0}^{1}\sum_{j_2=0}^{1}\cdots \sum_{j_n=0}^{1}\omega_N^{x \cdot \sum_{k=1}^{n}2^{n-k} \cdot j_k}|j_1j_2..j_n\rangle \tag{14}$$

Now let's rewrite $\omega_N^{x \cdot \sum_{k=1}^{n}2^{n-k} \cdot j_k}|j_1j_2..j_n\rangle$ in (14) with tensor product notation:
$$\frac{1}{\sqrt{N}}\sum_{j_1=0}^{1}\sum_{j_2=0}^{1}\cdots \sum_{j_n=0}^{1}\bigotimes_{k=1}^n\omega_N^{x \cdot \sum_{k=1}^{n}2^{n-k} \cdot j_k}|j_k\rangle \tag{15}$$
Putting $\sum_{j_1=0}^{1}\sum_{j_2=0}^{1}\cdots \sum_{j_n=0}^{1}$ inside $\bigotimes_{k=1}^n$ in (15) gives (16):
$$\frac{1}{\sqrt{N}}\bigotimes_{k=1}^n \sum_{j_k=0}^{1}\omega_N^{x \cdot \sum_{k=1}^{n}2^{n-k} \cdot j_k}|j_k\rangle \tag{16}$$
Expanding $\sum_{j_k=0}^{1}$ in (16) gives : $$\frac{1}{\sqrt{N}}\bigotimes_{k=1}^n (\omega_N^{x \cdot 2^{n-k}\cdot 0}|0\rangle + \omega_N^{x \cdot 2^{n-k} \cdot 1}|1\rangle) \tag{17}$$
In (17), $\omega_N^{x \cdot 2^{n-k}\cdot 0} = \omega_N^{0} = 1$. Therefore, (17) can be written as (18):
$$\frac{1}{\sqrt{N}}\bigotimes_{k=1}^n (|0\rangle + \omega_N^{x \cdot 2^{n-k}}|1\rangle) \tag{18}$$
In (18), $\omega_N^{x \cdot 2^{n-k}}$ can be written as $e^{\frac{2\pi i x}{2^k}}$ because :
$$\omega_N^{x \cdot 2^{n-k}} = e^{\frac{2\pi i x \cdot 2^{n-k}}{N}} = e^{\frac{2\pi i x \cdot 2^{n-k}}{2^n}} = e^{\frac{2\pi i x}{2^k}} \tag{19}$$

Putting (19) to (18) gives (20):
$$\frac{1}{\sqrt{N}}\bigotimes_{k=1}^n (|0\rangle + e^{\frac{2\pi i x}{2^k}}|1\rangle) \tag{20}$$
As (20) is (1), (1) has been derived.