Kontsevich gives a construction that produces deformation quantization of $C^\infty(M)$ for general Poisson manifolds $M$. The resulting formula (on $\mathbb{R}^n$) is $$ f\star g = \sum_{n=0}^\infty \hbar^n \sum_{\Gamma\in G_n}w_\Gamma B_{\Gamma,\alpha}(f,g) $$ where the weights $w_\Gamma$ come from some integration over a configuration space of points. It appears that this expansion ultimately comes from a path integral for the Poisson sigma model on a disk, as described in A path integral approach to the Kontsevich quantization formula by Cattaneo and Felder. They give the formula $$ f\star g(x) = \int_{X(\infty) = x} f(X(1))g(X(0))e^{\frac{i}{\hbar}S[X,\eta]}dX d\eta. $$ My question is: what is the physical intuition behind this "2D Poisson sigma model", and why should one expect that the above path integral gives a deformation of the classical observables, in the first place? How does the path-integral formulation relate to Kontsevich's original formula, indexed by certain graphs $\Gamma$?
2 Answers
To address the significance of the PSM: we can think of the Poisson algebra $C^\infty(M)$ as the algebra of observables of a one-dimensional classical field theory, which we can call "Poisson topological mechanics". Unfortunately, Poisson topological mechanics cannot be described in terms of an action functional, a fact which has to do with the degeneracy of Poisson structures by comparison to symplectic structures. On the other hand, the Poisson sigma model can be described in terms of a Lagrangian, and Poisson topological mechanics can be understood as a boundary condition for the PSM. I think it is widely believed that the PSM is the "universal" 2D topological field theory for which Poisson topological mechanics is a boundary condition. I learned this perspective from Degenerate Classical Field Theories and Boundary Theories by Butson and Yoo. In the introduction to A Classical Bulk-Boundary Correspondence (a video presentation on which can be found here), I try to articulate what this "universality" condition might mean in terms of "$\mathbb{P}_0$-factorization algebras." I should also mention A holography theory of Poisson sigma model and deformation quantization by Cui and Zhu and Towards holography in the BV-BFV setting by Mnev, Schiavina, and Wernli, where similar ideas are discussed. I am undoubtedly missing references, and I invite other commenters to add their own suggestions.
To address the relationship between the path integral formula of Cattaneo–Felder and the sum over graphs appearing in Kontsevich's paper: I point to Section 3.7 of the Cattaneo and Felder paper, where the details are spelled out. I will try to highlight a few of the salient parts of the argument. First, in QFT, whenever one computes correlation functions via a path integral, one performs a sum over Feynman diagrams, where the vertices of the Feynman diagrams are labeled by terms in the interaction functional and the edges of the diagrams are labeled by propagators. One also performs a "position-space" integral over the positions of the vertices. Since there are two types of field in the Poisson Sigma Model ($X$ and $\eta$), we can put arrows on all edges: an arrow going into a vertex means that the corresponding term in the interaction "eats" an $X$-field, and an arrow going out corresponds to an $\eta$-field. The interaction in the Poisson Sigma Model is determined by the Poisson bivector on $M$; it is in particular quadratic in $\eta$, and can have any degree of dependence on the $X$ field. This explains why Kontsevich's admissible graphs have to have two edges going out of any vertex, but any number can arrive at a vertex. Moreover, the Feynman weights "factor" into a part that depends on the source two-manifold (e.g. a disk or the upper half plane) and a part that depends on the target (e.g. $\mathbb{R}^d$). The part of the weights which depends on the source manifold is $w_\Gamma$, while the part that depends on the target manifold is $B_{\Gamma, \alpha}(f,g)$. The position-space integral is generally divergent in QFT (because the propagator is distributional in nature); because the PSM is topological, one can actually use configuration space techniques to show that the desired integrals converge. This trick is due, I think, to Fulton–MacPherson and Axelrod–Singer. That's the reason the configuration space integrals appear in Kontsevich's formula.
I don't know how much this answers your questions, but it might be helpful to observe that the formula you mention does reproduce the deformation of classical observables in quantum mechanics, non-perturbatively. Consider the usual phase space path integral from quantum mechanics \begin{equation} \langle q_f,t_f| q_i,t_i\rangle=\int_{q(t_i)=q_i}^{q(t_f)=q_f}Dq(t)Dp(t)\,e^{\frac{i}{\hbar}\int_{t_i}^{t_f}p\frac{dq}{dt}-H\,dt}\,. \end{equation} This involves integrating over maps $X:[t_i,t_f]\to T^*\mathbb{R}\,.$ Instead, one can integrate over all maps $X:S^1\to T^*\mathbb{R},$ where we put $H=0$ and consider the same observable Cattaneo and Felder do, ie. $f(X(1))g(X(0))\delta_{(p,q)}(X(\infty))\,.$ The result is \begin{equation} (f\star g)(p,q)=\int_{X:S^1\to T^*\mathbb{R},\,X(\infty)=(p,q)} DX\,f(X(1))g(X(0))e^{\frac{i}{\hbar}\int_{S^1}X^*pdq}\,. \end{equation} Note that, this is equal to the integral you wrote since the symplectic form on $T^*\mathbb{R}$ is exact (see remark 3 in the paper of Catteneo-Felder).
You can compute this similarly to how the first path integral is computed, ie. partition $S^1$ into $(N+1)$ pieces such that the marked points $\infty, 0,1$ are nodes in the partition, and denote the images of these nodes under the map $X:S^1\to T^*\mathbb{R}$ by $(p_{0},q_{0}),(p',q'),(p'',q'')\,,$ respectively. Then you have to compute \begin{equation} \lim_{N\to +\infty} \int_{(p_0,q_0)=(p,q)} \left\{\prod_{n=1}^{N} \mathrm{d}q_n\right\} \left\{\prod_{n=1}^{N} \frac{\mathrm{d}p_n}{2\pi\hslash}\right\} f(p'',q'')g(p',q') \exp\left[{\frac{i}{\hslash} \sum_{n=1}^N \frac{(p_{n-1}+p_n)}{2} {(q_n - q_{n-1})}}\right]\,. \end{equation} Note that, $(p',q')=(p_i,q_i), (p'',q'')=(p_j,q_j)$ for some $i,j\in \{1,\ldots, N\}\,,$ so $(p',q'),(p'',q'')$ are being integrated over.
This (heuristically) results in \begin{align} &(f\star g)(p,q)=\frac{1}{(4\pi\hbar)^2}\int_{\mathbb{R}^4}f(p'',q'')g(p',q')e^{\frac{i}{2\hbar}[(p''-p)(q-q')-(q''-q)(p-p')]}\,dp''\,dq''\,dp'\,dq'\,, \end{align} which is the non-perturbative form of the Moyal product (all of these computations are up to some constant factor). The exponent in the integral is there area of a triangle connecting $(p,q),(p',q'),(p'',q'')\,.$ This gives a strict deformation quantization, in the sense of Rieffel, and it can be used to determine the quantizations of the observables $x,p,$ along with any observable in $L^2(T^*\mathbb{R})\,.$
There are (secretly) interesting things involving $\infty$-groupoids that are relevant here.