Harmonic Oscillator in QM
Let's start with quantum mechanics. The Lagrangian of an harmonic oscillator is
\begin{equation}
L=\frac{m}{2}\dot{x}^2-\frac{k}{2}x^2
\end{equation}
The equations of motion are
\begin{equation}
\frac{d^2}{dt^2}x+\omega^2 x=0
\end{equation}
where $\omega=\sqrt{k/m}$. We know that when we quantize this theory, we get that the position and momentum operators are
\begin{aligned}
\hat{X}&=\frac{1}{\sqrt{2\omega}}\Big[\hat{a}+\hat{a}^\dagger\Big]\\
\hat{P}&=-i\sqrt{\frac{\omega}{2}}\Big[\hat{a}-\hat{a}^\dagger\Big]
\end{aligned}
The Hamiltonian, which classically was $H=\frac{p^2}{2m}+\frac{k}{2}x^2$ becomes
\begin{equation}
\hat{H}=\hbar \omega \big(\hat{N}+\frac{1}{2}\big)
\end{equation}
where $\hat{N}=\hat{a}^\dagger \hat{a}$. The interpretation is that $\hat{a}^\dagger$ creates a quantum of energy $\Delta E=\hbar \omega$. You can apply the operators $\hat{X}$ and $\hat{P}$ to measure the position and momentum of this state, or any state with any superposition of energies.
When we go into the wavefunction picture of this problem, we see that instead of having a position $x(t)$ we have a probability density $\psi(x,t)$ for each possible position. General states are superpositions of eigenstates of position
\begin{equation}
|\psi\rangle = \int dx \psi(x) |x\rangle
\end{equation}
where $|x\rangle$ is an eigenstate of the operator $\hat{X}$ (a.k.a. the Dirac delta function).
Notice that no one ever says that $\hat{X}$ creates a particle, although if you apply it to the vacuum, that is true. What is true is that $\hat{X}$ creates a superposition of destroying and creating a particle. That is, if you are in a state with $n$ particles, applying $\hat{X}$ will leave you in a superposition of having $n-1$and $n+1$ particles. Also, notice that it is clear from this that an eigenstate of $\hat{X}$ necessarily is a superposition of all eigenstates of $\hat{N}$, since that's the only way that decreasing and increasing $n$ will still give you the same state.
Free Field in QFT
Let's use a real Klein-Gordon scalar theory. The Lagrangian is
\begin{equation}
L=\frac{1}{2}\partial_\mu \phi \partial^\mu \phi
\end{equation}
The equations of motion are simply
\begin{equation}
\Big[\frac{\partial^2}{\partial t^2}-c^2\frac{\partial^2}{\partial x^2}\Big]\phi(x,t)=0
\end{equation}
Let's do a Fourier transform of the field $\phi(x,t)$
\begin{equation}
\phi(x,t)=\int \frac{d^3 k}{\big(2\pi\big)^{3/2}}\tilde{\phi}_k(t)e^{ik\cdot x}
\end{equation}
Plugging this into the equations of motion we get
\begin{equation}
\Big[\frac{d^2}{dt^2}+k^2c^2\Big]\tilde{\phi}_k(t)=0
\end{equation}
If we identify $k^2c^2=\omega^2$ we see that the Fourier components of the field satisfy an harmonic oscillator equation, for each wavenumber $\vec{k}$. Therefore, we can quantize them in exactly the same way
\begin{equation}
\hat{\tilde{\phi}}_k=\frac{1}{\sqrt{2\omega}}\Big[\hat{a}_k+\hat{a}_{-k}^\dagger\Big]
\end{equation}
The minus sign on the second term might be weird, but it follows from the constraint that $\phi$ must be real. Indeed, a general derivation would give a superposition of $\hat{a}$ and $\hat{b}^\dagger$ and it's only the extra constraint of being a real variable (classically) what makes $b$ equal to $a$. In particular, real functions have a fourier transform that satisfies $\tilde{\phi}^*_k=\tilde{\phi}_{-k}$.
Plugging this back into the field we get
\begin{equation}
\hat{\phi}(x)=\int \frac{d^3k}{\big(2\pi\big)^{3/2}\sqrt{2\omega_k}}\Big\{\hat{a}_k e^{ik\cdot x}+\hat{a}^\dagger_k e^{-ik\cdot x}\Big\}
\end{equation}
Interpretation of $\hat{\phi}$
The field operator $\hat{\phi}$ creates a superposition of a state with one more quanta with wavenumber $k$ and one less quanta with wavenumber $-k$ for all values of $k$. It might be a bit simpler to think about $\hat{\tilde{\phi}}_k$, since this one does the same job, but only for one value of $k$. This is literally what motivates the Wheeler-Feynman interpretation of antiparticles as being particles going back in time.
If we want to have a notion of where the particles we need to build something similar to the wavefunction. This is something that I haven't seen in many textbooks but it gives some intuition about what's going on. In QM, we go from $x$ to a wavefunction $\psi(x)$ that gives a probability for each classically allowed configuration $x$. In quantum field theory, we get a wave functional $\Phi[\phi]$ that gives a probability for each allowed classical configuration $\phi$. It is this wavefunctional $\Phi[\phi]$ that describes the state you are in. In QM, the wavefunction is
\begin{equation}
\psi(x)=\langle \psi | x\rangle
\end{equation}
It is the inner product between your state and the basis that represents what you identify as a "classical configuration": An eigenstate of position. In QFT, we can build the equivalent of $|x\rangle$ like this.
First, notice that our $\hat{\phi}(x,t)$ is an operator-valued function. This means that it is a function of coordinates $(x,t)$ and when you give it a point $(x,t)$ it gives you an operator. However, the rigorous way to define operators in QFT is to define operator-valued distributions. A distribution is a linear map that takes a function and gives you a number. The Dirac delta function is a distribution, for example. In this sense, the field operator $\hat{\phi}$ is actually
\begin{equation}
\hat{\phi}[f]=\int \frac{d^3k}{\big(2\pi\big)^{3/2}\sqrt{2\omega_k}}\Big\{\tilde{f}_k \hat{a}_k+\tilde{f}^*_k \hat{a}^\dagger_k\Big\}
\end{equation}
This has the very clear interpretation of creating a quantum state whose spectrum corresponds to the spectrum of the classical state $f(x)$. Therefore, we can identify the state it creates as
\begin{equation}
\hat{\phi}[f]|0\rangle=|f\rangle
\end{equation}
Now, if you have any given state $|\psi\rangle$, you can build your wavefunctional $\Phi[f]$ as
\begin{equation}
\Psi[f]= \langle \psi |f\rangle
\end{equation}
In this way, you can kind of visualize what classical field configurations build up your quantum state $|\psi\rangle$. In particular, if you pick $f_k=e^{ik\cdot x}$, that corresponds to $f(x)=\delta(x)$ so there is a very real sense in which $\hat{\phi}$ as defined with plane waves, creates a quantum superposition of infinitely many wavenumber states whose spectrum corresponds with the classical state of having a field $\phi$ localized at $x=0$. Moving $x$ around gets you any other position.
To answer the questions
what exactly is position in QFT and why is $\hat{\phi}|0\rangle$ interpreted as a state with a definite position?
There is no position operator, but you can identify your quantum state with a superposition of classical configurations $\phi(x)$. If your state corresponds to a field configuration that is a delta function, then you can interpret your quantum field to be localized there. $\hat{\phi}$ is interpreted as a state with definite position because it is the operator-valued distribution that corresponds to $f_k=e^{ik\cdot x}$ which are the Fourier modes of a Dirac delta function.
what is the interpretation of the correlation functions $\langle 0| \phi(x)\phi(y)|0\rangle$?
If we understand it in the context of operator-valued distributions, it literally creates and then destroys a particle at positions $x$ and $y$. Given that the Fourier modes $\hat{\tilde{\phi}}_k$ also have the interpretation of creating $k$ and destroying $-k$ the particle antiparticle Wheeler-Feynman interpretation is also compatible with this.
I hope some of this helps!
