Question. Let $N_k(n)$ count how many numbers $1\le x\le n$ for which $x$ has exactly $k$ unitary prime divisors, or equivalently $x$'s prime factorization has exactly $k$ primes with exponent $1$. Can we find an asymptotic for $N_k(n)$?
Motivation. I've been playing around to see what kinds of things I can (heuristically) calculate with Euler products. Particularly, the "probability" a number's prime factorization has certain features. A typical example is the probability a number is squarefree, which is $\prod (1-1/p^2)=1/\zeta(2)=6/\pi^2$.
For $v>1$, the probability a number $x$ has exactly $k$ prime factors $p$ with $v_p(x)=v$ is the sum of all Euler products
$$\prod_p \left(1-\frac{1}{p}\right)\left(\frac{1}{1-1/p}-\frac{1}{p^v}\right) $$
with exactly $k$ of their factors replaced with $\displaystyle\left(1-\frac{1}{p}\right)\frac{1}{p^v}$, or in other words
$$ e_k\left(\frac{\left(1-\frac{1}{p}\right)\frac{1}{p^v}}{\left(1-\frac{1}{p}\right)\left(\frac{1}{1-1/p}-\frac{1}{p^v}\right)}\right) \prod_p \left(1-\frac{1}{p}\right)\left(\frac{1}{1-1/p}-\frac{1}{p^v}\right) $$
where $e_k$ is the $k$th elementary symmetric "polynomial" of the infinitely many values in parentheses parametrized by primes $p$. The Newton-Girard identities allow us to write this as a finite polynomial in a number of infinite sums. We can specify the multiplicities of multiple valuations simultaneously if we replace elementary symmetric polynomials with more complicated monomial symmetric polynomials. Either way, the key is that everything converges.
This doesn't work for $v=1$ though, since the formula gives an indeterminate $0\cdot\infty$ expression. I think this case is broaching a barrier in the heuristic's applicability: the "probability measure" (which isn't one) on the naturals isn't countably additive. But the probability a number $x$ has zero primes $p$ with $v_p(x)=1$ (aka no unitary prime divisors) I think is given by the Euler product $\prod_p\left(1-\frac{1}{p}\right)(1+\frac{1}{p^2}+\cdots)=0$. I suspect the probability is $0$ for any $k>0$ too, but the Euler product approach isn't helping.
Basically, all you need is to evaluate the sum (the reason we put $<$ instead of $\le$ is merely to simplify the calculation):
$$ N_k(x)=\sum_{\substack{p_1<p_2<\cdots<p_k\\p_1p_2\cdots p_k<x}}1. $$
In particular, it follows from the prime number theorem that
$$ N_1(x)=\pi(x)+O(1)\sim{x\over\log x}. $$
Similar to what Chebyshev did to the case of $k=1$, we study the generalized sum to make our task easier:
$$ \vartheta_k(x)=\sum_{\substack{p_1<p_2<\cdots<p_k\\p_1p_2\cdots p_k<x}}\log(p_1p_2\cdots p_k), $$
where $\vartheta_1(x)=\vartheta(x)\sim x$ follows from the prime number theorem. Since the condition $p_1<p_2<\cdots <p_k$ is annoying, we remove them by multiplying both sides by $k!$:
\begin{align} k!\vartheta_k(x) &=\sum_{\substack{p_1,p_2,\dots,p_k\text{ distinct}\\p_1p_2\cdots p_k\le x}}\log(p_1p_2\cdots p_k) =\sum_{p_k<x}\sum_{\substack{p_1,p_2,\cdots,p_{k-1}\text{ distinct}\\p_1p_2\cdots p_{k-1}<x/p_k}}\log p_1p_2\cdots p_k \\ &=\sum_{p_k<x}\left[(k-1)!\vartheta_{k-1}\left(x\over p_k\right)+(k-1)!N_{k-1}\left(x\over p_k\right)\log p_k\right]. \end{align}
This suggests that
$$ \vartheta_k(x)=\frac1k\sum_{p<x}\left[\vartheta_{k-1}\left(\frac xp\right)+N_{k-1}\left(\frac xp\right)\log p\right].\tag1 $$
Having obtained the recursion formula, we see that
$$ \vartheta_2(x)=\frac12\sum_{p<x}\vartheta\left(\frac xp\right)+\frac12\sum_{p<x}\pi\left(\frac xp\right)\log p+O(x) $$
Suppose we choose $0<\theta<1$ and $y=x^\theta$. Then:
\begin{align} \sum_{p<x}\vartheta\left(\frac xp\right) &=\sum_{p<y}\vartheta\left(\frac xp\right)+\sum_{y\le p<x}\vartheta\left(\frac xp\right) \\ &=x\sum_{p<y}\frac1p+O\left\{\sum_{y\le p<x}\frac xp\right\} \\ &=x\log\log y+o(x\log\log y)+O\left(x\log{\log x\over\log y}\right) \\ &=x\log\log x+o(x\log\log x). \end{align}
Similarly, we can pick some $2<z<x$ such that
\begin{align} \sum_{p<x}\pi\left(\frac xp\right)\log p &=\sum_{p<z}\pi\left(\frac xp\right)\log p+O\left\{x\sum_{z\le p<x}{\log p\over p}\right\} \\ &=\sum_{p<z}\frac xp{\log p\over\log(x/p)}+O\left(x\log\frac xz\right) \\ &\ll x\int_2^z{\log t\over\log x-\log t}{\mathrm dt\over t\log t}+x\log\frac xz \\ &=x\int_{\log x/\log z}^{\log x/\log2}{\mathrm du\over u(u-1)}+x\log\frac xz \\ &\ll{x\log z\over\log x}+x\log\frac xz \end{align}
Setting $z=x/2$ gives
$$ \sum_{p<x}\pi\left(\frac xp\right)\log p=O(x) $$
Conclusively, we have
$$ \vartheta_2(x)\sim\frac12x\log\log x $$
Iterating (1) with $\vartheta_2(x)$ should yield a formula for general $\vartheta_k(x)$. Then, by partial summation one would obtain a formula for $N_k(x)$ via
$$ N_k(x)\sim{\vartheta_k(x)\over\log x} $$