Feynman rule of $\bar{\phi} \phi F_{\mu \nu}F^{\mu \nu}$ and its corresponding 4-photon scattering amplitude Consider the Lagrangian:
$$\mathcal{L}~=~-\frac{1}{2}\bar{\phi} \square \phi - \frac{1}{4}F_{\mu \nu}F^{\mu \nu} + \lambda \bar{\phi} \phi F_{\mu \nu}F^{\mu \nu}$$$\hspace{200px}$
The vertex Feynman rule I obtained from this Lagrangian is
$$i\lambda \left[(p\cdot k)g^{\mu \nu} - p^{\mu}k^{\nu} \right]$$
where $p^{\mu}$ and $k^{\nu}$ are the photon momenta.
When I consider the 1-loop diagrams of 4-photon interaction with incoming momenta $p_1^{\mu}, p_2^{\nu}$ and outgoing momenta $p_3^{\rho}, p_4^{\sigma}$ , I obtained the corresponding Feynman amplitude:
$$
\begin{alignat}{4}
  & \Big(\frac{s}{2} & \Big)^2 (g^{\mu \nu}g^{\rho \sigma}) & \int \frac{\mathrm{d}^4k}{\left(2\pi\right)^4}\frac{1}{k^2\left(k+p_s\right)^2} \\
+~ & \Big(\frac{t}{2} & \Big)^2 (g^{\mu \rho}g^{\nu \sigma}) & \int \frac{\mathrm{d}^4k}{\left(2\pi\right)^4}\frac{1}{k^2\left(k+p_t\right)^2} \\
+~ & \Big(\frac{u}{2} & \Big)^2 (g^{\mu \sigma}g^{\rho \nu}) & \int \frac{\mathrm{d}^4k}{\left(2\pi\right)^4}\frac{1}{k^2\left(k+p_u\right)^2} 
\end{alignat}
$$
where $p_s = p_1 + p_2$,  $p_t = p_1 - p_3$, and  $p_u = p_1 - p_4$.
My problem is, I cannot verify if my answer satisfies the Ward identity, which states that the amplitude should vanish when I replace one of the polarization vector with its momenta.
Is my amplitude wrong or is my Feynman rule wrong or if I'm right and just need to work harder to check?
 A: First of all, your Lagrangian has a weird prefactor of the kinetic term, there is usually no $1/2$ when we are dealing with complex scalar fields. Then, I would actually put a $1/2$ in front of the interaction term to cancel the combinatorial prefactors. Then
$$\mathcal{L} = \partial_\mu \phi^\ast \partial^\mu \phi - \frac{1}{4}  F^{\mu \nu} F_{\mu \nu} + \lambda \phi^\ast \phi F^{\mu \nu} F_{\mu \nu} $$
The Feynman rule for the vertex (all momenta incoming) is then given by
$$ 2 i \lambda (p^\mu k^\nu - g^{\mu \nu} p \cdot k) $$
Considering the process $\gamma \gamma \to \gamma \gamma$ at 1-loop we have 3 diagrams

Their sum (where the photon polarization vectors are truncated) is given by
$$\begin{align}
i \mathcal{M}^{\mu \nu \rho \sigma} &= \lambda^2 (s g^{\mu \nu} - 2 p_2^\mu p_1^\nu) (s g^{\rho \sigma} - 2 p_4^\rho p_3^\sigma) \int \frac{d^d k}{(2 \pi)^d} \frac{1}{k^2 (k+ p_s)^2} \\ 
& + \lambda^2 (t g^{\mu \rho} + 2 p_3^\mu p_1^\rho) (t g^{\nu \sigma} + 2 p_4^\nu p_2^\sigma) \int \frac{d^d k}{(2 \pi)^d} \frac{1}{k^2 (k+ p_t)^2} \\
& + \lambda^2 (u g^{\mu \sigma} + 2 p_4^\mu p_1^\sigma) (u g^{\nu \rho} + 2 p_3^\nu p_2^\rho) \int \frac{d^d k}{(2 \pi)^d} \frac{1}{k^2 (k+ p_u)^2}
\end{align}
$$
Now you can easily recognize that contracting this with $p_1^\mu$, $p_2^\nu$, $p_3^\rho$ or $p_4^\sigma$ yields zero. For example 
$$
\begin{align}
i \mathcal{M}^{\mu \nu \rho \sigma} p_{1,\mu} &= \lambda^2 (s p_1^\nu - 2  (s/2) p_1^\nu) (s g^{\rho \sigma} - 2 p_4^\rho p_3^\sigma) \int \frac{d^d k}{(2 \pi)^d} \frac{1}{k^2 (k+ p_s)^2} \\ 
& + \lambda^2 (t p_1^\rho + 2 (-t/2) p_1^\rho) (t g^{\nu \sigma} + 2 p_4^\nu p_2^\sigma)  \int \frac{d^d k}{(2 \pi)^d} \frac{1}{k^2 (k+ p_t)^2} \\
& + \lambda^2 (u p_1^\sigma + 2 (-u/2) p_1^\sigma) (u g^{\nu \rho} + 2 p_3^\nu p_2^\rho)\int \frac{d^d k}{(2 \pi)^d} \frac{1}{k^2 (k+ p_u)^2} = 0
\end{align}
$$
Your result does not satisfy the Ward identity, because the amplitude is incorrect. For some reason you have neglected the second piece of the vertex $2 i \lambda (p^\mu k^\nu)$.
