In a QFT with Euclidean signature, the correlation functions can only be well-defined in a time-ordered manner (This is Claim 1 on Page 2 of Simmons-Duffin's lecture note). For example, a scalar 2pt function $\langle O(x_1)O(x_2)\rangle$ $(x\equiv(\tau,\vec x))$should be regarded as
$$
\
\langle O_1( x_1)O_2( x_2)\rangle\equiv\langle\Omega|T_EO_1( x_1)O_2( x_2)|\Omega\rangle=
\begin{cases}
\langle\Omega|O_1(x_1)O_2(x_2)|\Omega\rangle,&\tau_1>\tau_2,\\
\langle\Omega|O_2(x_2)O_1(x_1)|\Omega\rangle,&\tau_1<\tau_2.
\end{cases}
$$
where $T_E$ is the Euclidean time-ordering, $|\Omega\rangle$ the vacuum.
On the other hand, correlation functions should be constrained by the global symmetries of the theory. For example, given a $d$-dimensional CFT and a conformal transformation $f(x)$, the correlation functions of primary operators should satisfy $$ \langle O_1(x_1)...O_n(x_n)\rangle=\Big|\frac{\partial f^{\mu}(x_1)}{\partial x_1^\nu}\Big|^{\Delta_1/d}...\Big|\frac{\partial f^{\mu}(x_n)}{\partial x_n^\nu}\Big|^{\Delta_n/d}\langle O_1(f(x_1))...O_n(f(x_n))\rangle.\tag I $$ where $\Delta_i$ is the scaling dimension of $O_i$.
My problem occurred when I was trying to derive (I) using the operator equation for primaries $$ U_fO(x)U_f^{-1}=\Big|\frac{\partial f^{\mu}}{\partial x^\nu}\Big|^{\Delta/d}O(f(x)),\tag{II} $$ where $U_f$ is the representation of the conformal transformation $f$ in Hilbert space. Here is my derivation (for simplicity, take $n=2$). Assuming $\tau_1>\tau_2$, one can remove the time-ordering $T_E$ and insert $U_f$ into the correlator $$ \begin{align} \langle O_1(x_1)O_2(x_2)\rangle=&\langle\Omega| O_1(x_1)O_2(x_2)|\Omega\rangle\\ =&\langle\Omega| U_fO_1(x_1)U_f^{-1} U_fO_2(x_2) U_f^{-1}|\Omega\rangle\\ =&\Big|\frac{\partial f^{\mu}(x_1)}{\partial x_1^\nu}\Big|^{\Delta_1/d}\Big|\frac{\partial f^{\mu}(x_2)}{\partial x_2^\nu}\Big|^{\Delta_2/d}\langle\Omega| O_1(f(x_1))O_2(f(x_2))|\Omega\rangle. \tag{III} \end{align} $$ From the second line to the third line, I use Eq. (II) and the fact that $U_f|\Omega\rangle=|\Omega\rangle$ (since $f$ is a global symmetry). Obviously, the time ordering is missing in the final result, and it should be wrong because even if we have $\tau_1>\tau_2$, it doesn't ensure $\tau_{f(x_1)}>\tau_{f(x_2)}$. So $\langle\Omega| O_1(f(x_1))O_2(f(x_2))|\Omega\rangle$ could be ill-defined for some choices of $f$ and $x_i$.
Where did I go wrong?
EDIT: In many textbooks (e.g., the Yellow Book (YB) by Di Francesco), Eq. (I) is derived by the Euclidean Path integral method. Assuming the conformal invariance of the action and the functional integration measure, Eq. (I) is merely a consequence of the conformal transformation law for primaries (II). (see the argument for (eq. 4.48) in YB). My main concern is to find an operator method to derive Eq. (I). However, what I derived was Eq. (III) when I attempted to do this.
The contradiction between Eq. (III) and Eq. (I) can somewhat be more clearly demonstrated using the Heaviside step function $\theta(x)$. Again, take $n=2$ for simplicity. The Euclidean 2pt functions can be written as follows $$ \begin{align} \langle O_1(x_1)O_2(x_2)\rangle=&\langle\Omega|T_EO_1(x_1)O_2(x_2)|\Omega\rangle\\=&\theta(\tau_1-\tau_2)\langle\Omega|O(x_1)O(x_2)|\Omega\rangle+\theta(\tau_2-\tau_1)\langle\Omega|O(x_2)O(x_1)|\Omega\rangle.\tag{IV} \end{align} $$ Then using the same trick (inserting $U_f^{-1}U_f$ into the correlator), one obtains $$ \begin{align} &\langle O_1(x_1)O_2(x_2)\rangle\\ =&\theta(\tau_1-\tau_2)\langle\Omega|U_f^{-1}U_fO_1(x_1)U_f^{-1}U_fO_2(x_2)U_f^{-1}U_f|\Omega\rangle\\ +&\theta(\tau_2-\tau_1)\langle\Omega|U_f^{-1}U_f O_2(x_2)U_f^{-1}U_f O_1(x_1)U_f^{-1}U_f|\Omega\rangle\\ =&\Big|\frac{\partial f^{\mu}(x_1)}{\partial x_1^\nu}\Big|^{\Delta_1/d}\Big|\frac{\partial f^{\mu}(x_2)}{\partial x_2^\nu}\Big|^{\Delta_2/d}\Big(\theta(\tau_1-\tau_2)\langle\Omega|O_1(f(x_1))O_2(f(x_2))|\Omega\rangle\\ +&\theta(\tau_2-\tau_1)\langle\Omega|O_2(f(x_2))O_1(f(x_1))|\Omega\rangle\Big).\tag{V} \end{align} $$ However, according to the definition (IV), we have $$ \begin{align} &\langle O_1(f(x_1))O_2(f(x_2))\rangle\\ =&\theta(\tau_{f(x_1)}-\tau_{f(x_2)})\langle\Omega|O_1(f(x_1))O_2(f(x_2))|\Omega\rangle\\ +&\theta(\tau_{f(x_2)}-\tau_{f(x_1)})\langle\Omega|O_2(f(x_2))O_1(f(x_1))|\Omega\rangle.\tag{VI} \end{align} $$
Comparing (V) and (VI), we may naively claim that in general $$ \langle O_1(x_1)O_2(x_2)\rangle\neq \Big|\frac{\partial f^{\mu}(x_1)}{\partial x_1^\nu}\Big|^{\Delta_1/d}\Big|\frac{\partial f^{\mu}(x_2)}{\partial x_2^\nu}\Big|^{\Delta_2/d}\langle O_1(f(x_1))O_2(f(x_2))\rangle, $$ since $f(x)$ in general does not need to be a chronological spacetime transformation, i.e, $\theta(\tau_{f(x_1)}-\tau_{f(x_2)})\neq \theta(\tau_{1}-\tau_{2})$.