What is the advantages of minimizing Wasserstein divergence rather than Pearson divergence in a GAN? Generative Adversarial Networks (GANs) are generative models that jointly train two neural networks: a discriminator, that learns to say apart real data from generated data, and a generator, that learns to produce synthetic data that is realistic enough to fool the discriminator.
Many types of GANs have been proposed, generally minimizing slightly different cost functions. For example, Least Square GANs minimize the Pearson divergence between real and fake data, while Wasserstein GANs minimize Wasserstein divergence. 
I wonder what are the advantages of using one divergence instead of the other. This also raises questions as: can we use gradient penalty with Least Square GANs?
 A: For your last question, Yes!, You can use gradient penalty with LSGAN too and you can read this. Gradient penalty is a trick and independent of what divergence/ distance you use. 
First one is longer though. TO reiterate, WGAN is trained with Wasserstein distance, not divergence! This is important as divergence is a weaker notion of distance due to the fact that divergence is not symmetric. Ex. $KL(p || q) \neq KL(q || p)$and in fact exploits different properties!
I'm gonna explain it as in the chronological order to understand why WGAN is important.
Before WGAN, GANs(Vanilla GAN, DCGAN , LSGAN and many other GANs before WGAN) , were trained to minimize an f-divergence (KL, JSD, Pearson...). If we take derivative of JSD with respect to generator parameters while real and generator data distribuitons are far from each other, gradient converges to zero. Very bad generator! Pearson divergence provides gradient to generators even if distributions are far away.
Wasserstein  OR Kantorovich-Rubinstein metric OR Earth Mover's Distance is the distance between two continuous probability distributions defined as 
$$ W(p_r, p_g) = \inf_{\gamma \sim \Pi(p_r, p_g)} \mathbb{E}_{(x, y) \sim \gamma}[\| x-y \|] $$
where $\Pi(p_r, p_g)$ is the set of all possible joint probability distributions between real and generator data distribution. $\gamma \sim \Pi(p_r, p_g)$ defines
So, what makes Wasserstein different than others?
Referring to WGAN paper, say we have two distributions, $\textit{P}$ and $\textit{Q}$ :
$$ \forall (x, y) \in P, x = 0 \text{ and } y \sim U(0, 1)\\
\forall (x, y) \in Q, x = \theta, 0 \leq \theta \leq 1 \text{ and } y \sim U(0, 1)\\$$ 
When $\theta \neq 0$ , there's no overlap: 
]1)
$$
\begin{aligned}
D_{KL}(P \| Q) &= \sum_{x=0, y \sim U(0, 1)} 1 \cdot \log\frac{1}{0} = +\infty \\
D_{KL}(Q \| P) &= \sum_{x=\theta, y \sim U(0, 1)} 1 \cdot \log\frac{1}{0} = +\infty \\
D_{JS}(P, Q) &= \frac{1}{2}(\sum_{x=0, y \sim U(0, 1)} 1 \cdot \log\frac{1}{1/2} + \sum_{x=0, y \sim U(0, 1)} 1 \cdot \log\frac{1}{1/2}) = \log 2\\
W(P, Q) &= |\theta|
\end{aligned}
$$
Wasserstein provides a smooth measure even if distributions are far from each other. This helps for a stable learning procedure, eliminating mode collapse, and improving the learnable class of manifolds (check this).
However, people don't use the Wasserstein metric as it is due to intractability of computing infimum. Using Kantorovich-Rubinstein duality: 
$$
W(p_r, p_g) = \frac{1}{K} \sup_{\| f \|_L \leq K} \mathbb{E}_{x \sim p_r}[f(x)] - \mathbb{E}_{x \sim p_g}[f(x)]
$$
to measure least upper bound on the function. Function needs to be K-Lipschitz continuous(Strongly advising this to read.) . 
Skipping some formalities ,say our function is from a family of K-Lipschitz continuois functions , $\{ f_w \}_{w \in W}$ , parameterized by $w$, Wasserstein distance is measured by : 
$$
L(p_r, p_g) = W(p_r, p_g) = \max_{w \in W} \mathbb{E}_{x \sim p_r}[f_w(x)] - \mathbb{E}_{z \sim p_r(z)}[f_w(g_\theta(z))]
$$
