Relationship between genetic diversity within and between species
Metrics of interest
The two metrics you are interested in are

$\pi$ - the mean number of differences between two randomly sampled (with replacement) alleles in a population
$d$ - the mean number of differences between two randomly sampled (with replacement) alleles coming from two different species

Consider two sequences

ATCGTCAAT
ATAGTTAAT

There are 2 pairwise differences between these two sequences (positions 3 and 6).
The whole point here is to understand that two individuals in the same population coalesce at a given time in the past just like two individuals coming from two different species. The number of pairwise differences is just equal to the rate at which mutations accumulate multiplied by the coalescence time. Let me develop this idea with a few equations below.
Neutral Expectations
Let's do the math! We will do two important assumptions below.

Every mutation makes a new allele (it is an infinite allele model)
All mutations are substitutions (no indels, no gene duplication, etc...)

Intro to Coalescence Theory
Let's first make sure you understand the concept of coalescence. Imagine you were to look at an evolving population backward in time. A coalescent event, is an event by which two lineages become one (looking backward in time). For simplicity consider a case of an asexual population (but future calculations don't do this assumption). In an asexual population, two siblings coalesced the previous generation, as in the previous generation the genes they are carrying were in a single individual (their parent).
We will see below how to calculate coalescence time between two randomly drawn individual from a population and how we can make inference about genetic diversity from this coalescent time.
Rate of accumulation of neutral mutations
The rate of accumulation of neutral mutation is also called the fixation rate of neutral mutations. Let $fix$ be this fixation rate. Consider a panmictic diploid population of constant size $N$. At a given locus of interest, the mutation rate is $\mu$. The number of mutations occurring every generation in this population is $2N\mu$ (The $2$ comes from the fact that in a diploid population, there are two homologous copies of every gene in each cell). If the mutation is neutral, then every single individual in the population has the same probability to end up being the parent of the whole population in a long time in the future. In other words, starting to look at this problem from the future and looking backward in time, any individual in a distance far enough in the future descends from a single individual in the present. As a consequence, the probability of a new mutation to reach fixation is $\frac{1}{2N}$. Multiplying the probability of fixation by the rate at which mutations occur in the whole population we obtain $\frac{1}{2N} 2N\mu = \mu$. In other words the rate of fixation of new neutral mutations is simply equal to the rate at which mutation occur.
Key equation:
$$fix=\frac{1}{2N} 2N\mu = \mu$$
$d$ - number of pairwise differences between randomly sampled (with replacement) alleles coming from different species
Consider two species which common ancestor lived $T$ generations ago. Each lineage accumulated mutations at rate $\mu$ (and explained above) and therefore the expected (average) total number of pairwise differences is $\bar d=2\mu T$. The probability of having exactly $d$ pairwise differences comes from a Poisson distribution with rate $2\mu T$ ($P(D=d) = Poisson(2\mu T)$).
Key equation: $$\bar d=2\mu T$$
$\pi$ - number of pairwise differences between randomly sampled (with replacement) alleles coming from the same population
Here the calculations follows the same logic as above $\pi = 2 \mu T$. The whole issue is that $T$ (the coalescence time) is unknown for the moment and we need to calculate it.
Let $P(T)$ be the probability that two randomly sampled (with replacement) individuals from a diploid panmictic population coalesce exactly $T$ generations ago. The probability of coalescing in a given generation is simply the probability for drawing the same individual, that is $\frac{1}{2N}$ and the probability of not coalescing in any subsequent generation is $1-\frac{1}{2N}$. Iterating over the generations we obtain
$$P(T) = \frac{1}{2N} \left( 1-\frac{1}{2N} \right)^T \approx \frac{1}{2N}e^{-\frac{T}{2N}}$$
, where $e \approx 2.7$ is Euler's number. The above approximation is accurate for large $N$ (say larger than 100). In case you are interest, you will find on this Math.SE post, an explanation for this approximation.
The expected value (average) $\bar T$ of this distribution is $\bar T=2N$ (and the variance is $var(T) = 4N^2$). As a consequence the expected number of pairwise differences between these two individuals is $\pi = 2 \mu 2 N = 4N\mu$
Quickly speaking, these calculations can be extended to calculate the coalescent time between $k$ randomly sampled alleles. The expected coalescent time is then $\bar T = \frac{4N}{k(k-1)}$ and therefore some algebra shows that the total expected time (along all branches) is $4 N \sum_{i=1}^{k-1}\frac{1}{i}$ and the expected number of segregating sites is $4 N \mu \sum_{i=1}^{k-1}\frac{1}{i}$.
Key equation: $$\pi = 4N\mu$$
Conclusion
It is clear from the above that $\pi = 4N\mu$ and $d=2\mu T$ are correlated as they are both linearly related to the mutation rate. As it is said in your quote

For neutral mutations, a positive association between $\pi$ and $d$ should exist, because the neutral theory predicts that both quantities are linearly proportional to the rate at which neutral mutations arise

The quote goes one saying

Recent genome-scale data shows instead that this association is in fact negative.

Such negative relationship between $d$ and $\pi$ cannot be caused by neutral processes. Therefore selection must be involved somehow.
Imagine for example that two species are both selected at different optima (opposite selection pressures). In such case, you would see purifying selection within species reducing the probability of a new mutation to reach fixation to a value lower than $\frac{1}{2N}$ and therefore reducing the rate of fixation $fix$ so that $fix&lt;\mu$. However, such opposite selection pressure among species will yield the species to diverge more than expected by random processes that is $\pi &gt; 4N\mu$. Such selective pattern could explain the observed negative relationship.
More information
This post offers book recommendation on the subject (population genetics).