This is a reference request.
Assume that $f_1,\ldots,f_r \in \mathbb{Z}[t]$ are non-zero linear polynomial.
Letting $\mu$ be the M\"{o}bius function, is there any work on $$ \sum_{p\leq x} \prod_{i=1}^r \mu(f_i(p))^2 \ \ ?$$
(One needs some obvious local conditions which should be built within the leading constant of the asymptotic I guess.)
Edit: I neglected to see that the OP asked for square-free values at prime arguments, in which some of the arguments below need to be modified. The necessary adjustments are found in the three papers linked below, where in each they dealt with the case of prime arguments.
Let $f(x) = f_1 \cdots f_r(x)$ with $f_j(x) = a_j x + b_j$ pairwise non-proportional linear polynomials with integer coefficients. Define the multiplicative function
$$\displaystyle \rho_f(n) = \# \{m \pmod{n} : f(m) \equiv 0 \pmod{n}\}$$
Define
$$\displaystyle N_f(X) = \#\{n \leq X : f(n) \text{ is square-free}\}$$ and define
$$\displaystyle N_f^{(1)}(X) = \# \{n \leq X : p^2 | f(n) \text{ implies } p > (\log X)/r\}.$$
One can show without much trouble by using simple properties of the Mobius function that, for any $\varepsilon > 0$,
$$\displaystyle N_f^{(1)}(X) = \prod_{p \leq (\log X)/r} \left(1 - \frac{\rho_f(p^2)}{p^2}\right) X + O_{f, \varepsilon}\left(X^{1 - 1/d + \varepsilon}\right).$$
The main term in the above expression tends to the expected main term for $N_f(X)$, which is
$$\displaystyle \prod_p \left(1 - \frac{\rho_f(p^2)}{p^2}\right)X.$$ Note that this product converges absolutely. Thus it suffices to estimate the term
$$\displaystyle N_f^{(2)}(X) = \# \{n \leq X : \exists p \text{ s.t. } p^2 | f(n) \text{ and } p > (\log X)/r\}.$$
To do so, let us examine the possibilities. There can be two reasons for $p^2$ to divide $f(x)$: either $p^2 | f_j(x)$ for some $1 \leq j \leq r$ or $p | f_i(x), p | f_j(x)$ for distinct indices $i$ and $j$. In the second case, $p$ divides the resultant of $f$ and $g$, which is just a rational integer (in fact, a divisor of the discriminant of $f$). Since $\Delta(f)$ is assumed to be non-zero, this can only happen for finitely many primes; and by taking $X$ large enough, we can ignore this possibility.
Thus we assume that $p^2 | f_j(x)$ for some $j$. For each $j$, define
$$\displaystyle N_{f,j}(X) = \# \{n \leq X : \exists p^2 | f_j(n) \text{ with } p > (\log X)/r\}.$$
Observe that $f_j(n) \ll X$ since $n \leq X$, and thus $p \ll \sqrt{X}$. By using the same argument which allows one to show that the number of square-free numbers up to $X$ is asymptotic to $(6/\pi^2)X$, one sees that $N_{f,j}(X) = O_{f_j}(X^{1/2})$ for $j = 1,2, \cdots, r$. Moreover, we see that
$$\displaystyle N_f^{(2)}(X) \leq \sum_{j = 1}^r N_{f,j}(X) = O_f(X^{1/2}),$$
and so $N_f^{(1)}(X)$ approximates $N_f(X)$, obtaining the correct order of magnitude.
For a review of such results, see:
https://arxiv.org/abs/1103.2028
https://people.maths.bris.ac.uk/~matdb/preprints/powerfree.pdf
