By expanding the covariant derivative of the Scalar QED lagrangian one gets the following term, sometimes called "seagull" vertex.
$$\mathcal{L}_{seagull} = -q^2A_\mu A ^\mu \phi^\dagger \phi$$
Most of the literature give $2iq^2g_{\mu\nu}$ as the Feynman rule for this vertex, saying that the factor $2$ is coming from the "symmetry of the two $A$ fields".
However, even if it actually makes sense at a first look, I can't reproduce the same result by contracting the $S$-matrix with appropriate initial and final states, since I get $iq^2g_{\mu\nu}$. In other words, I can't understand how the symmetry in the $A^\mu A_\mu$ term actually affects the $\mathcal{T}$-product on the $S$-matrix. How should I reproduce this result?
Here are my calculations:
I choose the following initial $|i\rangle$ and final $|f\rangle$ states, in order to evaluate $S_{fi}$:
$$|i\rangle = |\gamma_\lambda(k_1),\gamma_\sigma(k_2)\rangle$$
$$|f\rangle = |s(p_1),\bar s(p_2)\rangle$$
Where $\gamma_\lambda(k_1)$ describes a photon of momentum $k_1$ and polarization $\lambda$, $s(p_1)$ a scalar particle of momentum $p_1$ and $\bar s(p_2)$ a scalar anti-particle of momentum $p_2$. Now I can write:
$$S_{fi} = iq^2 g^{\mu\nu}\langle s(p_1),\bar s(p_2)|\int d^4x \, \mathcal{N}[A_\mu(x) A_\nu(x) \phi^\dagger(x) \phi(x)]\, |\gamma_\lambda(k_1),\gamma_\sigma(k_2)\rangle$$
Expanding each field in terms of its "positive energy" and "negative energy" parts, and selecting the only contributing term one gets:
$$S_{fi} = iq^2 g^{\mu\nu}\langle s(p_1),\bar s(p_2)|\int d^4x \,\mathcal{N}[A_{\mu+} A_{\nu+} \phi^\dagger_- \phi_-]\, |\gamma_\lambda(k_1),\gamma_\sigma(k_2)\rangle$$
Making use of the $\mathcal{N}$-product, applying the fields to the states and finally integrating over $d^4x$ one gets the final amplitude:
$$S_{fi} = (2\pi)^4 \delta^4(p_1+p_2-k_1-k_2)iq^2 g^{\mu\nu} \epsilon_{\mu,\lambda}\epsilon_{\nu,\sigma}$$
From which one can deduce the vertex factor $iq^2g^{\mu\nu}$. The "symmetry" in the $A$ fields actually reflects in the fact that I could have obtained $\epsilon_{\nu,\lambda}\epsilon_{\mu,\sigma}$ instead of $\epsilon_{\mu,\lambda}\epsilon_{\nu,\sigma}$, but since the $\mu\nu$ indexes are contracted through $g^{\mu\nu}$ I would guess that it would not make any difference. Where am I failing?