Proof. At either endpoint \(x=X,2X\), the range \(q\le X^\sigma\) is contained in the range of (15) for every fixed logarithmic saving once \(X\) is sufficiently large. Removing proper prime powers from a \(\psi\)-sum costs at most \[\sum_{p^j\le3X,\ j\ge2}\log p\ll X^{1/2}L^2\] in each residue class. Indeed there are \(O(X^{1/2}L)\) pairs \(p,j\ge2\), and each summand is \(O(L)\). Translating the two endpoints through any \(|a|\le H_X\) costs \(O((H_X+1)L)\), uniformly in the residue class. Thus, writing \(Q=\lfloor X^\sigma\rfloor\), for every \(A_0>0\), \[ \sum_{q\le Q}E_X(q) \ll_{A_0}X L^{-A_0} +Q\{X^{1/2}L^2+(H_X+1)L\} \ll_{A_0,\sigma,D,B}X L^{-A_0}. \tag{17}\] The last absorption uses the fixed strict inequality \(\sigma<1/2\). The same bound controls the maximum over all shifts at once; no union bound over shifts is used.

An interval of length \(X\) contains at most \(X/q+1\) integers in one residue class. Also \[\frac q{\varphi(q)}=\prod_{p\mid q}\frac p{p-1} \le2^{\omega(q)}\le d(q).\] Since \(q\le Q<X\) for large \(X\), these two facts give \[ E_X(q)\ll \frac Xq\{L+d(q)\}. \tag{18}\] For every fixed real \(T\ge1\), \[ \sum_{q\le Q}\frac{\mu^2(q)T^{\omega(q)}}q \le\prod_{p\le Q}(1+T/p)\ll_T L^{c_T} \tag{19}\] for a finite constant \(c_T\). The elementary bound \(\sum_{p\le y}1/p=O(\log\log(3y))\) suffices: for \(y\ge2\), compare \(1/p\le e p^{-1-1/\log y}\) with \(\log\zeta(1+1/\log y)\), and bound \(\zeta(1+u)\le1+1/u\). No prime-number theorem is needed for this estimate.

As squarefree \(q\) satisfy \(d(q)=2^{\omega(q)}\), (18)–(19) imply, for some finite \(c_K\), \[\sum_{q\le Q}\mu^2(q)K^{2\omega(q)}E_X(q) \ll_K X L^{c_K}.\] Cauchy–Schwarz, applied to the nonnegative numbers \(E_X(q)\), therefore gives \[\begin{split} \sum_{q\le Q}\mu^2(q)K^{\omega(q)}E_X(q) &\le \left(\sum_{q\le Q}E_X(q)\right)^{1/2} \left(\sum_{q\le Q}\mu^2(q)K^{2\omega(q)}E_X(q)\right)^{1/2}\\ &\ll X L^{(c_K-A_0)/2}. \end{split}\] Choose \(A_0\ge c_K+2A\). This proves (16). ◻

Proof. Expand the product of divisor sums. Möbius factors kill all nonsquarefree divisors. Whenever the product of the \(F_r\)’s is nonzero, the support hypothesis gives \[ \prod_{\ell=1}^{j_r}d_{r,\ell}\le X^{\rho_r} \quad\text{for every }r, \qquad q\le\prod_\gamma d_\gamma\le X^\sigma. \tag{23}\] This uses the budget of the sum of the coordinates in each factor, and does not impose a budget on the dimension.

For a squarefree tuple, the conditions \(d_\gamma\mid m+b_\gamma\) prescribe at prime \(p\mid q\) the residues \(m\equiv-b_\gamma\pmod p\) for all selected coordinates. These congruences are inconsistent unless the selected shifts agree modulo every \(p\mid q\). When they agree, the Chinese remainder theorem prescribes exactly one class \(m\pmod q\). In the unmarked case this is precisely compatibility, and the count is \(X/q+O(1)\), with an absolute constant in the \(O(1)\).

When \(M\ge1\), the number of positive integer tuples whose product is at most \(Y\ge1\) is at most \[ Y\left(\sum_{d\le Y}\frac1d\right)^{M-1} \le Y(1+\log Y)^{M-1}. \tag{24}\] Indeed, fix the first \(M-1\) entries and bound the number of possible last entries by \(Y\) divided by their product, then enlarge each preceding summation to \(d\le Y\). Apply (24) with \(Y=X^\sigma\), multiply the \(O(1)\) count error by \(\prod_r\|F_r\|_\infty\), and divide by \(X\). This proves (22), hence (21) in the unmarked case. If \(M=0\), the unmarked average and (20) are both exactly \(\prod_rF_r\).

For the marked case, put \(n=m+a_0\). On a residue class satisfying these divisor congruences, the residue of \(n\pmod p\) is \(a_0-b_\gamma\) whenever \(p\mid d_\gamma\). Consequently this class is reduced modulo \(q\) exactly when the additional marked compatibility condition holds. If it is not reduced, a prime \(n\) in that class would have to be a prime divisor of \(q\), hence satisfy \(n\le q\). But (23) and \(\sigma<1/2\) give \[n\ge X-H_X> X^\sigma\ge q\] for large \(X\). The marked sum on a nonreduced class is therefore exactly zero, not an unestimated exceptional term. On a reduced class, the marked sum is \[\frac X{\varphi(q)}+O(E_X(q)).\] Since \(q\) is squarefree, \(\varphi(q)=\prod_{p\mid q}(p-1)\), proving the main density factor in (20).

For fixed squarefree \(q\), each prime dividing \(q\) must be assigned to a nonempty subset of the \(M\) divisor coordinates. There are at most \((2^M-1)^{\omega(q)}\le (2^M)^{\omega(q)}\) choices if \(M\ge1\); support and compatibility only reduce this count. The empty tuple, if \(M=0\), has \(q=1\) and is handled directly by the same bound with \(K=1\). Thus the total normalized error is at most \[\frac{\prod_r\|F_r\|_\infty}{X} \sum_{q\le X^\sigma}\mu^2(q)(2^M)^{\omega(q)}E_X(q),\] which is \(O_A(L^{-A})\) by Lemma 9. This proves the marked case, including all scalar factors and the wholly zero-dimensional case. ◻

Proof of Proposition 8. Let \(M=|\Gamma|\). By Lemma 10, it remains to evaluate \(\mathcal A_\delta\) from (20): the arithmetic error is \(O_N(L^{-N})\) for every fixed \(N\). If \(M=0\), that sum is exactly \(\prod_rF_r\), and the assertion follows. We henceforth assume \(M>0\).

Choose a smooth compact extension \(\widetilde F_r\) of each positive-dimensional function and define its Fourier transform by \[\Phi_r(\mathbf u)=(2\pi)^{-j_r}\int_{\mathbb R^{j_r}} e^{\sum_\ell v_\ell}\widetilde F_r(\mathbf v) e^{\mathrm i\mathbf u\cdot\mathbf v}\,d\mathbf v.\] The function \(\Phi_r\) is Schwartz, and Fourier inversion gives exactly \[ F_r(\mathbf v)=\int_{\mathbb R^{j_r}}\Phi_r(\mathbf u) e^{-\sum_\ell(1+\mathrm i u_\ell)v_\ell}\,d\mathbf u \qquad(\mathbf v\in[0,\infty)^{j_r}). \tag{25}\] For scalar factors use the scalar itself. Put \(\Phi=\prod_r\Phi_r\), \(z_\gamma=(1+\mathrm i u_\gamma)/L\), and \(z_S=\sum_{\gamma\in S}z_\gamma\).

Insert (25) into the density sum (20), retaining the sum over all squarefree tuples. The following bound justifies moving that sum inside the integrals: \[ \sum_{\mathbf d}\rho_\delta(\mathbf d)\prod_\gamma d_\gamma^{-1/L} \le\prod_p(1+C_Mp^{-1-1/L})\ll_M L^{C_M}. \tag{26}\] Indeed \(1/(p-\delta)\le2/p\), every nonempty selected subset contributes at most \(p^{-1/L}\), and there are at most \(2^M-1\) subsets. Moreover \(\sum_pp^{-1-\varepsilon}\le\log\zeta(1+\varepsilon) \le\log(1+1/\varepsilon)\). Multiplication by the finite integral \(\int|\Phi|\) proves absolute convergence of the sum-integral interchange.

Call a nonempty subset \(S\subset\Gamma\) allowed at \(p\) if its shifts satisfy the same compatibility and mark-exclusion conditions. The exact Euler factor is then \[P_p(\mathbf z)=1+\frac1{p-\delta} \sum_{S\text{ allowed at }p}(-1)^{|S|}p^{-z_S}, \qquad \mathcal A_\delta=\int\Phi(\mathbf u)\prod_pP_p(\mathbf z)\,d\mathbf u.\] By (26) and Schwartz decay, the part of this integral where \(\max_\gamma|u_\gamma|>\sqrt L\) is \(O_N(L^{-N})\) for every fixed \(N\). All subsequent comparisons are made on the complementary region.

The factors carrying the zeta poles are \[B_p(\mathbf z)=\prod_{i=1}^t\prod_{\varnothing\ne S\subset\mathcal M_i} (1-p^{-1-z_S})^{-(-1)^{|S|}},\qquad H_p(\mathbf z)=P_p(\mathbf z)/B_p(\mathbf z).\] For \(\Re z_\gamma\ge0\), all denominators in this expression have absolute value at least \(1/2\), and \(|H_p(\mathbf z)|\le\exp(C_M/p)\). If distinct shifts in \(\mathcal H\) have distinct residues modulo \(p\), the allowed subsets are exactly the nonempty subsets of the individual fibers \(\mathcal M_i\). Matching the terms of order \(p^{-1}\) gives \[ H_p(\mathbf z)=1+O_M(p^{-2}), \tag{27}\] uniformly for \(\Re z_\gamma\ge0\).

At \(\mathbf z=0\), summing signs over the nonempty subsets of each occupied allowed residue class gives \[P_p(0)=\frac{p-\nu_p(\mathcal H)}{p-\delta},\qquad B_p(0)=(1-1/p)^t.\] It follows that \(H_p(0)\) is the local factor of \(\mathfrak S(\mathcal H)\). This identity holds also when that factor is zero.

We now compare the product of the \(H_p\)’s additively with its value at zero. Take \(Y=\max(2,2DL^B)\). For \(p\le Y\), differentiation along the segment from zero to \(\mathbf z\) gives \[|H_p(\mathbf z)-H_p(0)|\ll_M L^{-1/2}\frac{\log p}{p}.\] For example, \(|p^{-z_S}-1|\le|z_S|\log p\) on that segment, while the denominators stay bounded away from zero. Each mixed partial product in the additive telescoping identity is at most \[\exp\left(C_M\sum_{p\le Y}1/p\right)\ll_M(\log Y)^{C'_M}.\] Using even the elementary bound \(\sum_{p\le Y}(\log p)/p\le\sum_{n\le Y}(\log n)/n=O((\log Y)^2)\), we obtain an error \(O(L^{-1/2}(\log L)^A)\) for the finite product. For \(p>Y\) there are no collisions of distinct shifts, so (27) shows that both remaining tail products are \(1+O_M(1/Y)\). Consequently \[ \prod_pH_p(\mathbf z)=\mathfrak S(\mathcal H) +O\!\left(L^{-1/2}(\log L)^A\right),\qquad \mathfrak S(\mathcal H)\ll(\log L)^A. \tag{28}\] This argument never divides by a possibly zero local singular-series factor.

It remains to evaluate the zeta factors. Absolute Euler products for \(\Re z_\gamma>0\) give \[\prod_pB_p(\mathbf z)=\prod_{i=1}^t \prod_{\varnothing\ne S\subset\mathcal M_i} \zeta(1+z_S)^{(-1)^{|S|}}.\] Write \(w_\gamma=1+\mathrm i u_\gamma\) and \[K_i(\mathbf u)=\prod_{\varnothing\ne S\subset\mathcal M_i} \left(\sum_{\gamma\in S}w_\gamma\right)^{(-1)^{|S|+1}}.\] The pole expansion \(\zeta(1+z)=z^{-1}(1+O(|z|))\) and the identity \(\sum_{\varnothing\ne S\subset\mathcal M_i}(-1)^{|S|}=-1\) imply \[ \prod_pB_p(\mathbf z) =L^{-t}\prod_iK_i(\mathbf u)\{1+O_M(L^{-1/2})\}. \tag{29}\] Every denominator in \(K_i\) has positive integer real part, so the product of the \(K_i\) has at most polynomial growth. Thus \[\mathcal C=\int_{\mathbb R^M}\Phi(\mathbf u)\prod_iK_i(\mathbf u)\,d\mathbf u\] is absolutely convergent. Integrating (28) and (29) against the Schwartz transform, then restoring the negligible frequency tails, proves (13). The constant is real by conjugation under \(\mathbf u\mapsto-\mathbf u\).

Finally, a singleton fiber contributes \(w_\gamma\), and a double fiber contributes \[\frac{w_\gamma w_{\gamma'}}{w_\gamma+w_{\gamma'}} =w_\gamma w_{\gamma'}\int_0^\infty e^{-(w_\gamma+w_{\gamma'})y_i}\,dy_i.\] After taking absolute values, the new variables contribute \(e^{-2\sum_i y_i}\); the remaining factor is a Schwartz function times a fixed polynomial. Fubini therefore applies. Recombining the Fourier integrals by differentiating (25) proves (14). ◻

Proof. The cases \(s=0,1\) are immediate. Otherwise let \(y=(1/4)\log h\) and \(Q=\prod_{p\le y}p\). The prime number theorem gives \(Q\le h^{1/2}\) for large \(h\). For all tuples, including those with equal entries, define \[\mathfrak S_y(\mathbf a)=\prod_{p\le y} (1-\nu_p(\{a_1,\ldots,a_s\})/p)(1-1/p)^{-s}.\] The exponent remains \(-s\) even on tuples with collisions. For independent uniform residues \(a_1,\ldots,a_s\) modulo a prime \(p\), averaging over an additional uniform residue \(x\) gives \[\mathbb E(1-\nu_p/p) =\mathbb P(a_i\ne x\text{ for every }i)=(1-1/p)^s.\] The Chinese remainder theorem therefore shows that \(\mathfrak S_y\) has exactly mean one on complete residue tuples modulo \(Q\).

A uniform integer in a length-\(h\) interval has residue distribution modulo \(Q\) within total variation \(O(Q/h)\) of uniform. The elementary bound \(\sum_{p\le y}1/p=O(\log\log(3y))\) used above gives \[0\le\mathfrak S_y(\mathbf a) \le\prod_{p\le y}(1-1/p)^{-s}\ll_s(\log y)^{A_s}.\] It follows that \[h^{-s}\sum_{a_i\in I_i}\mathfrak S_y(\mathbf a) =1+O_s\big((Q/h)(\log y)^{A_s}\big).\] There are \(O_s(h^{s-1})\) tuples with equal entries, so deleting them changes this normalized mean by at most \(O_s(h^{-1}(\log y)^{A_s})\).

For a distinct tuple put \(\Delta=\prod_{i<j}|a_i-a_j|\), a positive integer. For \(y>2s\) every factor of the tail \[T_y(\mathbf a)=\prod_{p>y}(1-\nu_p/p)(1-1/p)^{-s}\] is positive. A prime not dividing \(\Delta\) has \(\nu_p=s\) and a local factor \(1+O_s(p^{-2})\); the logarithm of a remaining factor is \(O_s(1/p)\). Hence \[|\log T_y|\ll_s \frac1y+\sum_{\substack{p>y\\p\mid\Delta}}\frac1p \le C_s\left(\frac1y+\frac{\log\Delta}{y\log y}\right).\] Since \(\log\Delta\le\binom s2\log(Kh)\), this gives, uniformly in the distinct tuples, \[T_y(\mathbf a)=1+O_{s,K}(1/\log\log h).\] Now use the identity \(\mathfrak S(\{a_1,\ldots,a_s\})=\mathfrak S_y(\mathbf a)T_y(\mathbf a)\) and the nonnegative truncated mean. This identity is valid also if a small-prime factor vanishes; no quotient by \(\mathfrak S_y\) is taken. The stated error follows from \(Q\le h^{1/2}\) and \(\log y\asymp\log\log h\). ◻

Proof of Proposition 5. First we check the smoothness and support needed by the moment formula. Extend \(f_j\) by zero to \(\mathbb R^j\). Its cumulative integral is smooth on \(\mathbb R^j\), and multiplication by smooth coordinate cutoffs that equal one near \([0,\infty)\) and vanish below \(-1\) makes it compactly supported without changing its orthant values or derivatives. Indeed, the cumulative integral already vanishes when any coordinate exceeds an upper support bound of \(f_j\). On the nonnegative orthant, \(F_j\) vanishes if \(\sum_i v_i\ge\tau\), and its signed complete mixed derivative is given by (4). This derivative is zero on every coordinate face. For \(j\ge1\), \(Tf_{j+1}\) has the same smoothness, interior-support, and face-vanishing properties.

Expand \(Z^2\) over pairs of subsets \(S,T\subset I\). Any shift in \(S\mathbin\triangle T\) occurs once and sets a coordinate of one of the functions in (4) to zero. Its main coefficient is therefore zero. If \(S=T\) has size \(j\), the coefficient is \(\|f_j\|_2^2\). Averaging the singular series and using \[\frac{\binom hj}{L^j}\longrightarrow\alpha_j\] proves the first limit in (5). With a prime mark in \(J\), the same matching rule applies and \[\frac{h\binom hj}{L^{j+1}}\longrightarrow\lambda\alpha_j.\] The two divisor factors have total support \(2\tau<1/2\), so the marked moment formula proves the first limit in (6).

For an internal mark \(a\in I\), the exact deletion identity (8) replaces \(F_j\) by \(\widetilde F_j=F_j+F_{j+1}(\cdot,0)\) on the prime event. Deletion is a bijection, so no extra dimension factor occurs. The support radius of \(\widetilde F_j\) is still \(\tau\), and its signed complete mixed derivative is \(f_j+Tf_{j+1}\). It vanishes on coordinate faces when \(j\ge1\). Apply the external marked formula to these functions. Only equal subsets \(U\) survive, and \[\frac{h\binom{h-1}{j}}{L^{j+1}}\longrightarrow\lambda\alpha_j.\] This proves the formula for \(v\), including the scalar term \(\lambda|f_0+\int f_1|^2\).

The equality-pattern estimate preceding this subsection justifies summing the errors in these limits. Applied to \(Z^4\), the same argument uses four divisor factors, total support \(4\tau<1\), and only finiteness of their pattern constants. Bounded averages of the singular series then give the second assertion in (5).

It remains to prove a detector estimate whose constant is independent of the chosen family. Write \(q_X=\log(2X+2h)/L=1+o(1)\). At a relevant prime \(m+b\), we have \(D_G(m;b)=1\). Consequently, pointwise, \[\frac{\vartheta(m+b)}L\le q_XD_G(m;b)^2, \qquad \left(\frac{\vartheta(m+b)}L\right)^2 \le q_X\frac{\vartheta(m+b)}L.\] The first inequality holds at composites as well, since its left side is zero. Because \(Z^2\ge0\), these inequalities give \[\begin{align*} \mathbb E_X(Z^2V_J^2) &\le q_X\mathbb E_X(Z^2V_J)\\ &\quad+q_X^2\sum_{\substack{b,c\in J\\b\ne c}} \mathbb E_X\bigl[Z^2 D_G(m;b)^2D_G(m;c)^2\bigr]. \end{align*}\] An off-diagonal summand contains six divisor factors, with total support \(6\tau<1\). The shifts \(b,c\) are paired and lie outside \(I\). The two subsets from \(Z^2\) must again be identical. For their common size \(j\), the coefficient is \(c_G^2\|f_j\|_2^2\). The pair \((b,c)\) is ordered, and \[\frac{h(h-1)\binom hj}{L^{j+2}} \longrightarrow\lambda^2\alpha_j.\] Thus the off-diagonal sum tends to \(\lambda^2c_G^2w\). Together with (6), this proves (7). Although the convergence rates may depend on the fixed family, the resulting \(C_*\) does not. ◻