  Let $U^N$ be a family of $N\times N$ independent Haar unitary random matrices
and their adjoints, $Z^N$ a family of deterministic matrices, $P$ a
self-adjoint noncommutative polynomial, i.e. such that for any $N$,
$P(U^N,Z^N)$ is self-adjoint, $f$ a smooth function. We prove that for any $k$,
if $f$ is smooth enough, there exist deterministic constants
$\alpha_i^P(f,Z^N)$ such that $$ \mathbb{E}\left[\frac{1}{N}\text{Tr}\left(
f(P(U^N,Z^N)) \right)\right]\ =\ \sum_{i=0}^k \frac{\alpha_i^P(f,Z^N)}{N^{2i}}\
+\ \mathcal{O}(N^{-2k-2}) .$$ Besides the constants $\alpha_i^P(f,Z^N)$ are
built explicitly with the help of free probability. As a corollary, we prove
that given $\alpha<1/2$, for $N$ large enough, every eigenvalue of $P(U^N,Z^N)$
is $N^{-\alpha}$-close to the spectrum of $P(u,Z^N)$ where $u$ is a $d$-tuple
of free Haar unitaries. We also prove the convergence of the norm of any
polynomial $P(U^N\otimes I_M, I_N\otimes Y^M)$ as long as the family $Y^M$
converges strongly and that $M\ll N \ln^{-5/2}(N)$.
