In Deformation quantization of Poisson manifolds, Kontsevich gives the quantization formula
$$f \star g = \sum_{n=0}^\infty \hbar^n \sum_{\Gamma \in G_n} w_\Gamma B_{\Gamma,\alpha}(f,g).$$
He gives an explicit formula up to $O(\hbar^3)$:
$$f \star g = fg + \hbar \sum_{i,j} \alpha^{ij} \partial_i(f)\partial_j(g) + \frac{\hbar^2}{2} \sum_{i,j,k,l} \alpha^{ij}\alpha^{kl} \partial_i\partial_k(f)\partial_j\partial_l(g) + \frac{\hbar^2}{3}\left(\sum_{i,j,k,l} \alpha^{ij}\partial_j(\alpha^{kl})(\partial_i\partial_k(f)\partial_l(g) - \partial_k(f)\partial_i\partial_l(g))\right) + O(\hbar^3).$$
This is an equality up to gauge equivalence. Other than that, the equality $w_\Gamma B_{\Gamma,\alpha} = w_{\Gamma'} B_{\Gamma',\alpha}$ for graphs that differ only in the labeling (by virtue of skew symmetry in the weight and the Poisson bivector) has been used to collect like terms, and one term vanishes (also explained in the preceding link).
Anyhow,
We are assured that this star product is associative up to second order, but I don't see it:
We can write the above star product up to $O(\hbar^3)$ in terms of graphs without losing any information:
Here it is understood that the order of the ground vertices is fixed (first is $f$, second is $g$), and all arrows correspond to independent indices; the labels $L$ and $R$ just indicate whether the contraction is with the first or second upper index of $\alpha$.
Computing $f \star (g \star h)$ and $(f \star g) \star h$ is an exercise in the Leibniz rule, which can also be done graphically:
Of course, many terms do match identically (struck out with green). Some terms match because $1/2 + 1/2 = 1$ (red underline 1 and 2). Some other terms match because swapping $L$ and $R$ once gives a sign change (black underline). But what happens with the boxed terms?
(I've verified using a computer program that the expressions above are correct.)
The desired equality is the following:
In terms of bidifferential operators, this is: $$\sum_{i,j,k,l} \alpha^{ij}\partial_i(\alpha^{kl})\partial_k(f)\partial_l(g)\partial_j(h) = \sum_{i,j,k,l} \alpha^{ij}\partial_j(\alpha^{kl})\partial_i(f)\partial_k(g)\partial_l(h).$$
We can relabel the indices on the right hand side (such that the last 3 match) and use skew symmetry: $$\sum_{i,j,k,l} \alpha^{ij}\partial_i(\alpha^{kl})\partial_k(f)\partial_l(g)\partial_j(h) = \sum_{i,j,k,l} \alpha^{ik}\partial_i(\alpha^{jl})\partial_k(f)\partial_l(g)\partial_j(h).$$
The only difference is the position of $j$ and $k$ in the first two factors.
If I'm not mistaken, the Jacobi identity for the Poisson bivector is $$0 = [\alpha,\alpha]^{ijk} = \alpha^{li}\partial_l(\alpha^{jk}) + \alpha^{lj}\partial_l(\alpha^{ki}) + \alpha^{lk}\partial_l(\alpha^{ij}).$$
Using this on the factor in the right hand side gives $$\alpha^{ik}\partial_i(\alpha^{jl}) = -\alpha^{ij}\partial_i(\alpha^{lk}) - \alpha^{il}\partial_i(\alpha^{kj}) = \alpha^{ij}\partial_i(\alpha^{kl}) + \alpha^{il}\partial_i(\alpha^{jk}),$$ but this forces $\alpha^{il}\partial_i(\alpha^{jk}) = 0$, which is nonsense. What am I doing wrong?