  We establish convergence in norm and pointwise almost everywhere for the
non-conventional (in the sense of Furstenberg) bilinear polynomial ergodic
averages \[ A_N(f,g)(x) := \frac{1}{N} \sum_{n =1}^N f(T^nx) g(T^{P(n)}x)\] as
$N \to \infty$, where $T \colon X \to X$ is a measure-preserving transformation
of a $\sigma$-finite measure space $(X,\mu)$, $P(\mathrm{n}) \in \mathbb
Z[\mathrm{n}]$ is a polynomial of degree $d \geq 2$, and $f \in L^{p_1}(X), \ g
\in L^{p_2}(X)$ for some $p_1,p_2 > 1$ with $\frac{1}{p_1} + \frac{1}{p_2} \leq
1$. We also establish an $r$-variational inequality for these averages (at
lacunary scales) in the optimal range $r > 2$. We are also able to "break
duality" by handling some ranges of exponents $p_1,p_2$ with
$\frac{1}{p_1}+\frac{1}{p_2} > 1$, at the cost of increasing $r$ slightly. This
gives an affirmative answer to Problem 11 from Frantzikinakis' open problems
survey for the Furstenberg--Weiss averages (with $P(\mathrm{n})=\mathrm{n}^2$),
which is a bilinear variant of Question 9 considered by Bergelson in his survey
on Ergodic Ramsey Theory from 1996. This also gives a contribution to the
Furstenberg-Bergelson-Leibman conjecture. Our methods combine techniques from
harmonic analysis with the recent inverse theorems of Peluse and Prendiville in
additive combinatorics. At large scales, the harmonic analysis of the adelic
integers $\mathbb A_{\mathbb Z}$ also plays a role.
