Browning in this paper proves that if $f \in \mathbb{Z}[x]$ is an irreducible polynomial of degree $d \geq 3$ and $k$ is an integer $\geq 3d/4 + 1/4$, then we have $$\#\{n \in \mathbb{Z} \cap [1, X] : f(n) \text{ is $k$free}\}\sim C_{1}(k, f)X$$ as $X \rightarrow \infty$ where $$C_{1}(k, f) = \prod_{p}\left(1  \frac{\varrho(p^{k})}{p^{k}}\right)$$ and $\varrho_{f}(n) = \#\{a \pmod{n}: f(a) \equiv 0 \pmod{n}\}$. Is there anywhere that gives a similar asymptotic for the case when $f$ is reducible?

The reducible case is not much harder than the irreducible case, if one uses an argument in Greaves's paper (used originally by Gouvea and Mazur) "Powerfree values of binary forms". In particular, the main term can be defined just as easily for reducible polynomials $f(x) \in \mathbb{Z}[x]$. The error terms can be dealt with simply as follows. Write $f(x) = f_1(x) \cdots f_s(x)$, where $f_i(x)$ is irreducible over $\mathbb{Z}$ and define the error term $E_i(X) = \# \{x \in \mathbb{Z} : p^k  f(x) \text{ for some } x > \xi\}$, where $\xi = \frac{1}{k} \log X$. Then the overall error term is just $E(X) = \sum_{j=1}^s E_i(X)$. Thus one can use whatever techniques one would like (in this case, Browning's result is the strongest; although using a different determinant method I am able to reproduce the same result) to estimate each of the $E_i(X)$ and multiply the result by $s$, which is negligible since it is an absolute constant and the error term will be a power saving over the main term (assuming the constant in front of the main term is nonzero). I should remark that the above estimates do not cover the case when say $f_1(x) = uv^{kl}$ and $f_2(x) = u'v^l$, or when three factors of $f$ conspire to give us a non$k$free term. However, this situation is easily shown to give a negligible contribution. To see this, consider the equations $f_1(x) = uv^{kl}$, $f_2(x) = u'v^l$ and $f(x) = wv^k$. With $x$ and $v$ fixed, we see that $u, u'$ are divisors of $w$ and hence there are no more than $d(w) = O(w^\epsilon)$ many choices. A simple analysis on the size of $w$ yields that $w = O(X^d)$, and so $d(w) = O(X^\epsilon)$. We may then crudely estimate the number of points lying on the intersection of the varieties defined by $f_1(x) = uv^{kl}$ and $f_2(x) = u'v^l$ where we take $u, u'$ to be fixed as follows. For $u, u'$ fixed, each of the equations define a curve in $\mathbb{A}^2$, and they have at most finitely many intersection points. Thus the overall contribution is $O(X^\epsilon)$ and is negligible. 

