Strehl identity for the sum of cubes of binomial coefficients

In 1993 Strehl showed that $$ \sum_k\binom nk^3=\sum_k\binom nk^2\binom{2k}n. $$ I’m interested in a combinatorial proof.

Upd (Jan '14). Maybe the original question was too restrictive — I'm now asking for any proof.

I've written a combinatorial proof (see answers), but it's somewhat convoluted. Maybe there is a nice proof using generating functions / residues (something like one for the Dixon's identity), for example?..

Comments and thoughts

1. There is of course an algebraic proof (any such identity can be proved using WZ-magic of hypergeometric functions AFAIK). Feel free to post any proof, but I’m mainly interested in a combinatorial proof.
2. It’s well known that $\sum_k\binom nk=2^n$ and $\sum_k\binom nk^2=\binom{2n}n$. Perhaps, one should think about Stehl identity as the next member of this family. (Is there any analogue for degree 4, I wonder, btw?)
3. Quadratic case suggests that it’s, perhaps, better to write LHS as $\binom nk^2\binom n{n-k}$. And maybe instead of one $n$ one should introduce two or three parameters (s.t. Strehl identity corresponds to n=m[=l]). A naive approach is $\sum\binom nk\binom mk\binom m{m-k}$ but looks like $\sum\binom nk\binom mk\binom n{m-k}$ might be a better choice — actually, it's equal to $\sum\binom nk\binom mk\binom{2k}m$.
4. arXiv:0712.3946 might be of some relevance. There LHS is interpreted as the number of ways to deal a deck of 3n cards (n denominations, 3 colors) to 3 players so that no player receives a card of her own color.
5. Let's write the identity in the form $$ \sum_{i+j=n}\binom Bi\binom Wj\binom ni=\sum_k\binom Bk\binom Wk\binom{2k}n, $$ where $B$ and $W$ are two $n$-element sets. What I want to say is that both sides count ways to choose an $n$-element subset of $B\sqcup W$ with some kind of... extension (cf. two formulas for a trinomial coefficient, perhaps).

Suppose we have $n$ black cards and $n$ white cards and we want to divide cards into «good» and «bad» — in such way that there is the same number of bad cards of each color — and burn $n$ of bad cards.

(RHS) Obviously there are exactly $\sum\binom nk^2\binom{2k}n$ ways to do it: we choose $k$ bad black cards, $k$ bad white cards and burn $n$ bad cards.

(LHS) If we just want to choose what $n$ cards we burn — it can be done in $\sum_k\binom nk\binom n{n-k}$ ways: we choose $k$ black cards to burn and $n-k$ white cards to burn.

Let’s call bad white cards and good black cards «perverse». After we’ve burned that $n$ cards we need to choose exactly $k$ perverse cards (from the cards left): if we have $l$ good black cards, we've had $n-l$ bad white cards of which $(n-l)-(n-k)=k-l$ left.

So $$\sum\binom nk^2\binom{2k}n=\sum\binom nk\binom n{n-k}\binom nk.$$

// This is of course just a rewrite of an algebraic proof (based on Vandermonde convolution and aforementioned two presentations of the same trinomial coefficient) — namely, of the proof on p. 16 of V. Strehl’s «Binomial Identities...» (I’ve tried to look into the paper before — but missed it).

Suppose we seek to show that $$\sum_{k=0}^n {n\choose k}^3 = \sum_{k=\lceil n/2 \rceil}^n {n\choose k}^2 {2k\choose n}.$$

With

$${n\choose k} {2k\choose n} = \frac{(2k)!}{k!\times (n-k)! \times (2k-n)!} = {2k\choose k} {k\choose n-k}$$

we find that the RHS is $$\sum_{k=\lceil n/2 \rceil}^n {n\choose k} {2k\choose k} {k\choose n-k}.$$

Introduce $${2k\choose k} = \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{(1+z)^{2k}}{z^{k+1}} \; dz$$

and (this integral is zero when $0\le k\lt \lceil n/2\rceil$) $${k\choose n-k} = \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{(1+w)^{k}}{w^{n-k+1}} \; dw$$

to get for the RHS

$$\frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z} \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{1}{w^{n+1}} \sum_{k=0}^n {n\choose k} \frac{w^k(1+w)^k (1+z)^{2k}}{z^k} \; dw\; dz \\ = \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z} \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{1}{w^{n+1}} \left(1+\frac{w(1+w)(1+z)^2}{z}\right)^n \; dw\; dz \\ = \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{1}{w^{n+1}} (z+w(1+w)(1+z)^2)^n \; dw\; dz \\ = \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \frac{1}{2\pi i} \int_{|w|=\gamma} \frac{1}{w^{n+1}} (z+w(z+1))^n (1+w(z+1))^n \; dw\; dz.$$

Extracting first the residue in $w$ in next the residue in $z$ we get

$$\frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{1}{z^{n+1}} \sum_{q=0}^n {n\choose q} z^{n-q} (1+z)^q {n\choose n-q} (1+z)^{n-q} \; dz \\ = \sum_{q=0}^n {n\choose q}^2 \frac{1}{2\pi i} \int_{|z|=\epsilon} \frac{(1+z)^n}{z^{q+1}} \; dz \\ = \sum_{q=0}^n {n\choose q}^3$$

QED.

Addendum May 27 2018. We compute this using formal power series as per request in comment. Start from

$${2k\choose k} = [z^k] (1+z)^{2k}$$

and

$${k\choose n-k} = [w^{n-k}] (1+w)^k.$$

Observe that this coefficient extractor is zero when $n-k\gt k$ or $k\lt\lceil n/2\rceil$ where $k\ge 0.$ Hence we are justified in lowering $k$ to zero when we substitute these into the sum and we find

$$\sum_{k=0}^n {n\choose k} [z^k] (1+z)^{2k} [w^{n-k}] (1+w)^k \\ = [z^0] [w^n] \sum_{k=0}^n {n\choose k} \frac{1}{z^k} (1+z)^{2k} w^k (1+w)^k \\ = [z^0] [w^n] \left(1+\frac{(1+z)^2 w (1+w)}{z}\right)^n \\ = [z^n] [w^n] (z + (1+z)^2 w (1+w))^n \\ = [z^n] [w^n] (1+w(1+z))^n (z+w(1+z))^n.$$

We extract the coefficient on $[w^n]$ then the one on $[z^n]$ and get

$$[z^n] \sum_{q=0}^n {n\choose q} (1+z)^q {n\choose n-q} (1+z)^{n-q} z^q \\ = \sum_{q=0}^n {n\choose q}^2 [z^{n-q}] (1+z)^n = \sum_{q=0}^n {n\choose q}^2 {n\choose n-q} = \sum_{q=0}^n {n\choose q}^3.$$

The claim is proved.