Equations 29–30: The Nyman–Beurling Criterion and Báez-Duarte's Strengthening — A Deep Dive

1. What the Criterion Says

The Nyman–Beurling criterion is one of the most remarkable reformulations of the Riemann Hypothesis (RH): it converts a statement about zeros of an analytic function into a pure approximation problem in a Hilbert space — no complex analysis appears in the statement at all.

Theorem (Nyman 1950, Beurling 1955). Let ρ_θ(x) = {θ/x} − θ{1/x} for 0 < θ ≤ 1, where {·} denotes the fractional part. Then:

RH is true ⟺ the linear span of { ρ_θ : 0 < θ ≤ 1 } is dense in L²(0,1).

Beurling proved the more general L^p version: the span is dense in L^p(0,1) ⟺ ζ(s) has no zeros in the half-plane Re s > 1/p. Taking p = 2 gives RH.

Theorem (Báez-Duarte 2002). It suffices to take θ = 1/k with k = 1, 2, 3, …:

RH is true ⟺ the indicator function χ_(0,1] lies in the L²(0,∞)-closure of span{ f_k(x) = {1/(kx)} : k ∈ ℕ }.

This is a genuine strengthening because a continuum of dilation parameters is replaced by a countable, arithmetic set — the criterion becomes a statement about Dirichlet-polynomial approximation.

2. Why It Works (the Mellin Mechanism)

The Mellin transform is the bridge. Since ρ_θ vanishes identically on (1,∞) (for x > 1 both fractional parts are the arguments themselves and cancel), a direct computation using the classical identity

$$\int_0^\infty \{1/x\}\,x^{s-1}\,dx = -\frac{\zeta(s)}{s}, \qquad 0 < \operatorname{Re} s < 1,$$

(substituting u = θ/x for the first term) gives

$$\int_0^1 \rho_\theta(x)\,x^{s-1}\,dx = \frac{\theta - \theta^{s}}{s}\,\zeta(s).$$

So every finite linear combination Σ c_j ρ_{θ_j} has Mellin transform ζ(s) · A(s)/s where A is entire and bounded appropriately. By the Mellin–Plancherel theorem, L²(0,1) is isometric (up to constants) to the Hardy space H² of the half-plane Re s > ½ restricted to the critical line with measure |ds|/|s|². Density of the span means one can make

$$\left\| \frac{1}{s} - \frac{\zeta(s)A(s)}{s} \right\|_{L^2(\text{critical line},\, dt/(\frac14+t^2))} \to 0.$$

• (⇐) If ζ(ρ) = 0 with Re ρ > ½, then every approximant ζ·A vanishes at ρ. Point evaluation is a continuous functional on H² (reproducing-kernel property), so the limit function 1 would also vanish at ρ — contradiction. Hence density ⟹ no zeros off the line.
• (⇒) The converse (RH ⟹ density) is the deep direction; Beurling's proof uses his theory of invariant subspaces for the semigroup of multiplications/dilations on H².

In Báez-Duarte's version, A(s) = Σ_{k≤N} c_k k^{−s} is an honest Dirichlet polynomial, so RH becomes: Dirichlet polynomials times ζ(s) can approximate 1 in a weighted L² sense on the critical line.

3. History and the People Behind It

• Bertil Nyman (Swedish mathematician, doctoral student of Arne Beurling at Uppsala) proved the L² criterion in his 1950 thesis, "On the One-Dimensional Translation Group and Semi-Group in Certain Function Spaces" (Uppsala University, 1950). The thesis is famously hard to find; the result circulated mostly through Beurling.
• Arne Beurling (1905–1986, Uppsala, later IAS Princeton — also famous for breaking the German Geheimschreiber cipher in WWII) published the L^p generalization in 1955 in a 3-page PNAS note.
• Vasilii Vasyunin (Steklov Institute, St. Petersburg, 1995) showed rational θ suffice and derived beautiful closed-form cotangent-sum formulas for the inner products ⟨ρ_{1/p}, ρ_{1/q}⟩ — these "Vasyunin sums" are now a research topic in their own right (reciprocity formulas by Bettin–Conrey).
• Luis Báez-Duarte (1938–2018, IVIC, Venezuela) proved in 2002 that θ ∈ {1, ½, ⅓, …} suffices.
• Báez-Duarte, Balazard, Landreau & Saias (2000) studied the distance $$d_N^2 = \inf_{A_N} \frac{1}{2\pi}\int_{-\infty}^{\infty} \left|1 - \zeta A_N\left(\tfrac12+it\right)\right|^2 \frac{dt}{\tfrac14+t^2}$$ and conjectured the asymptotic $$d_N^2 \sim \frac{C}{\log N}, \qquad C = \sum_{\rho} \frac{1}{|\rho|^2} = \sum_\rho \frac{1}{\frac14+\gamma^2} = 2 + \gamma_{\text{Euler}} - \log 4\pi \approx 0.046191.$$ Note C = 2λ₁, twice the first Li coefficient — a pretty cross-connection between the NB and Li criteria.
• Jean-François Burnol (2002) proved the sharp lower bound: lim inf d_N² · log N ≥ Σ_{Re ρ = ½} m_ρ²/|ρ|² (m_ρ = multiplicity). So convergence, if it happens, is provably no faster than C/log N — this is exactly why numerical NB experiments plateau and cannot "verify" RH empirically.
• Bettin, Conrey & Farmer (2013) showed that, conditionally on RH plus a mild zero-spacing hypothesis, the choice c_k = μ(k)(1 − log k/log N) achieves the conjectured asymptotic.

4. Key Papers (with links)

(arXiv links are stable; please double-check journal DOIs when citing.)

5. Computational Examples in Python

Example 1 — Verify the Mellin identity, and watch it vanish at a zeta zero

For θ = 1/k the substitution t = 1/x makes the integral an exactly summable series: on each interval (n, n+1) the integrand is constant, and summation by parts gives Σ({n/k}−{(n−1)/k})n^{−s} = (θ − θ^s)ζ(s).

import mpmath as mp
mp.mp.dps = 25

def mellin_lhs(s, k, N=200000):
    """ integral_0^1 rho_{1/k}(x) x^{s-1} dx  computed as an exact series:
        sum_{n>=1} {n/k} * (n^{-s} - (n+1)^{-s}) / s,  plus mean-tail correction. """
    total = mp.mpf(0)
    for n in range(1, N):
        total += mp.mpf(n % k) / k * (mp.power(n, -s) - mp.power(n + 1, -s))
    mu = mp.mpf(k - 1) / (2 * k)              # mean of {n/k}
    total += mu * mp.power(N, -s)             # tail estimate
    return total / s

def mellin_rhs(s, k):
    th = mp.mpf(1) / k
    return (th - th**s) / s * mp.zeta(s)

s = mp.mpc('0.6', '3.0'); k = 3
print("LHS =", mp.nstr(mellin_lhs(s, k), 8))
print("RHS =", mp.nstr(mellin_rhs(s, k), 8))

# Now evaluate AT the first nontrivial zero: transform must vanish
rho1 = mp.zetazero(1)                          # 0.5 + 14.134725...i
print("At first zero:", mp.nstr(mellin_rhs(rho1, k), 8))   # exactly 0
print("Series there :", mp.nstr(mellin_lhs(rho1, k), 8))   # ~ 0 numerically

Example 2 — Compute the Báez-Duarte distance d_N² (exact piecewise Gram matrix)

We compute d_N² = 1 − bᵀG⁻¹b with exact closed-form integration between breakpoints (the integrand is rational on each piece), where G_{kl} = ∫₀^∞ {t/k}{t/l} t⁻² dt and b_k = (log k + 1 − γ)/k (exact).

import numpy as np
from math import gcd, log, floor, pi

EULER = 0.57721566490153286

def gram_entry(k, l, Tmin=20000):
    L = k * l // gcd(k, l)
    T = L * max(1, -(-Tmin // L))                      # T = multiple of period L
    bp = np.union1d(np.arange(0, T + 1, k), np.arange(0, T + 1, l)).astype(float)
    total = 0.0
    for t1, t2 in zip(bp[:-1], bp[1:]):
        a, b = floor((t1 + t2) / (2 * k)), floor((t1 + t2) / (2 * l))
        H = lambda t: (t/(k*l) - (b/k + a/l)*log(t) - a*b/t) if t > 0 else 0.0
        total += H(t2) - H(t1)
    # tail: periodic mean mu over one period, tail ~ mu / T
    mu, bpp = 0.0, np.union1d(np.arange(0, L + 1, k), np.arange(0, L + 1, l)).astype(float)
    for t1, t2 in zip(bpp[:-1], bpp[1:]):
        a, b = floor((t1 + t2) / (2 * k)), floor((t1 + t2) / (2 * l))
        P = lambda t: t**3/(3*k*l) - (b/k + a/l)*t**2/2 + a*b*t
        mu += P(t2) - P(t1)
    return total + (mu / L) / T

# sanity check: G_11 = log(2*pi) - gamma = 1.2606...
print("G_11 =", gram_entry(1, 1), " expected:", log(2*pi) - EULER)

C = 2 + EULER - log(4 * pi)                     # = 0.0461914... = 2*lambda_1
for N in [1, 2, 4, 8, 12]:
    G = np.array([[gram_entry(k, l) for l in range(1, N+1)] for k in range(1, N+1)])
    b = np.array([(log(k) + 1 - EULER) / k for k in range(1, N+1)])
    c = np.linalg.lstsq(G, b, rcond=None)[0]    # lstsq: G is ill-conditioned
    dN2 = 1 - b @ c
    print(f"N={N:3d}   d_N^2 = {dN2:.6f}   conjectured C/log N = "
          f"{C/log(N) if N > 1 else float('inf'):.6f}")

You will see d_N² decrease very slowly — Burnol's theorem guarantees it can never beat C/log N. This is worth emphasizing in your paper: NB numerics cannot converge fast, by theorem.

Example 3 — The constant C = Σ_ρ 1/(¼+γ²) = 2 + γ − log 4π

import mpmath as mp
mp.mp.dps = 20

target = 2 + mp.euler - mp.log(4 * mp.pi)
print("target C =", mp.nstr(target, 12))          # 0.0461914179...

s = mp.mpf(0)
for n in range(1, 101):
    g = mp.im(mp.zetazero(n))
    s += 2 / (mp.mpf(1)/4 + g**2)                 # factor 2: pair rho, conj(rho)
# density-based tail estimate: (log(g/2pi) + 1)/(pi*g)
tail = (mp.log(g / (2*mp.pi)) + 1) / (mp.pi * g)
print("100 zeros:", mp.nstr(s, 8), " + tail est:", mp.nstr(s + tail, 8))

The sum over the first 100 zeros plus a zero-density tail estimate reproduces C ≈ 0.04619, confirming the identity Σ_ρ 1/(ρ(1−ρ)) = 2 + γ − log 4π (which follows from the Hadamard product for ξ(s); note C = 2λ₁, Li's first coefficient).

6. Points Worth Making in Your Paper

1. Conceptual payoff: NB turns RH into functional analysis — density of dilations of the fractional part in L². It connects RH to Beurling's invariant subspace theory and to the "Hilbert–Pólya" circle of ideas.
2. Arithmetic re-enters through the coefficients: the optimal c_k are mollified Möbius values, c_k ≈ μ(k)(1 − log k/log N) (Bettin–Conrey–Farmer) — the Möbius function is hiding inside the "pure analysis" problem.
3. The log N wall: Burnol's lower bound is the rigorous explanation for why every numerical NB experiment stalls near d_N² ≈ 0.05–0.1 for feasible N — an honest limitation to report.
4. Related directions: de Roton extended the criterion to the Selberg class; Vasyunin/cotangent sums connect to the Estermann function and quantum modular forms (Bettin–Conrey).