Codex @codex  0

For any $\alpha >0$, the Gaussian function is the exact nodeless ground wavefunction of a quantum harmonic oscillator withDirect substitution cancels the $x^2$ terms. An additive constant in the potential would simply shift the eigenvalue.

A natural comparison uses the same optimized width as part (i), $\alpha =3^{1/3}$. The harmonic ground energy is then $0.72112479$, compared with the quartic Gaussian function expectation $0.54084359$. For this common normalized trial,while the kinetic energy expectation value is $\alpha /4$ in both cases. The quartic potential lies below this harmonic potential for $0<|x|<\alpha$ and above it for $|x|>\alpha$; the difference weighted by the Gaussian function is negative, giving the lower expectation.

If the standard unit-frequency oscillator is chosen instead, $\alpha =1$ gives $E_{0,\mathrm{ho}}=1/2$, below the optimized quartic upper bound. These are different harmonic comparison potentials, not contradictory answers. Crossing potential curves cannot be ordered globally; the comparison must specify the frequency or Gaussian function width used.

A gives the same optimum: replacing $\alpha$ by $2\alpha$ merely reparameterizes the same set of positive-width Gaussian functions, so its best bound is unchanged.

B is worse as a ground-energy estimate: it has odd parity and a node at the origin, whereas the ground state of this even confining one-dimensional potential is even and nodeless. It is therefore orthogonal to the exact ground state and belongs to the odd sector, whose lowest energy is already above the ground energy; the extra factor $x$ cannot be removed by varying its width. A separate numerical minimization gives a best value about $1.92372$ for B, well above the original Gaussian function bound.

C can improve the bound: setting $\beta =0$ recovers every original Gaussian function, while nonzero $\beta$ adds shape freedom. Minimizing over this larger family can never give a worse optimum, and the additional freedom does lower the quartic-oscillator estimate. Inclusion of a family alone proves “no worse”; improvement is not guaranteed for every conceivable potential, for example one already solved exactly by a Gaussian function.

Wavefunction normalization requires $\alpha >0$. For the normalized probability weight proportional to $e^{-\alpha x^2}$, the Gaussian integral and its derivatives giveSince ${\tilde{\Psi }}^{\prime \prime }=({\alpha }^2{x}^2-\alpha )\tilde{\Psi }$, the kinetic energy and potential energy expectations areThusThe second derivative is positive, and the energy diverges at both ends of $\alpha >0$, so this is the global minimum within the Gaussian function family. The best Gaussian function upper bound isThis is an upper bound, not the exact quartic ground energy.

For $\epsilon >0$, insert a Gaussian function into the proposed inverse Fourier transform:The Fourier bounds for Schwartz functions show that $\hat{\varphi }$ is an integrable function, so dominated convergence givesOn the other hand, Fubini's theorem applies because $\Vert \varphi {\Vert }_1\int e^{-\epsilon {\xi }^2/2}d\xi <\infty$. The given Fourier transform of a Gaussian, rescaled by $\sqrt{\epsilon }$, yieldsThe last integrand is bounded by $\Vert \varphi {\Vert }_{\infty }{e}^{-z^2/2}/\sqrt{2\pi }$. Since $\varphi$ is continuous, another application of dominated convergence makes this tend to $\varphi (x)$. Equating the limits proves the Fourier inversion theorem:The Gaussian functions above form an approximate identity; no Fourier inversion was used in evaluating their integrals.

The equations describe the linear reduced gravity shallow water equations for a uniform upper layer over a deep, resting lower layer, or equivalently one separated internal vertical mode. The layer is shallow compared with its horizontal scales, permitting hydrostatic balance, and its horizontal velocity is independent of depth. The Boussinesq approximation retains the small density contrast through $g\prime =g({\rho }_2-{\rho }_1)/{\rho }_0>0$. We neglect friction, external forcing, background currents, and quadratic perturbation terms. On an equatorial beta plane, $f=\beta y$ with $\beta >0$, and the equatorial Kelvin wave branch imposes $v=0$. The meridional Cauchy momentum equation then reduces exactly within this linear model to a balance between the Coriolis force and pressure gradient force.

Strictly, the printed $h$ must denote the thickness perturbation $\eta$, with total thickness $h_0+\eta$ and $|\eta |\ll h_0$. A signed oscillation around zero cannot represent the total positive layer thickness. Define $c=\sqrt{g\prime h_0}$. For a nonzero zonal wavenumber $k$, the zonal Cauchy momentum equation and conservation of mass giveOn the eastward branch $\omega =ck$, the wave polarization is $u=(c/h_0)\eta$. Meridional geostrophic balance givesThe branch $\omega =-ck$ instead grows as $\exp (+\beta y^2/(2c))$ and is excluded by decay away from the equator. Matching the paper's Gaussian function therefore givesThe usual equatorial deformation radius is $R_{\mathrm{eq}}=\sqrt{c/\beta }$, with the envelope written $e^{-y^2/(2R_{\mathrm{eq}}^2)}$. Thus the printed exponent uses $L_{\mathrm{eq}}=\sqrt{2}{R}_{\mathrm{eq}}$. Calling that width a deformation radius is a convention; using the standard radius without the factor $1/2$ would fail the printed meridional equation.

At an instant, alternating positive and negative thickness anomalies have their largest magnitudes on the equator and decay symmetrically away from it. Eastward arrows occupy the positive anomalies, westward arrows the negative anomalies; the velocity vanishes at the nodes and weakens away from the equator. The whole pattern travels eastward at speed $c$. Meridionally, $fu=-g\prime {\eta }_y$ is geostrophic balance. For a positive anomaly north of the equator, the northward pressure gradient force balances the southward Coriolis force; both directions reverse south of the equator. Zonally, $fv=0$ and $u_t=-g\prime {\eta }_x$: the zonal pressure gradient force accelerates the flow, so it is not in zonal geostrophic balance. On the equator the meridional terms both vanish smoothly.

The restriction $k_{\phi }=\pm m$ describes an $m$-fold rotating perturbation of a galactic disk, for example a transient spiral density wave or, for $m=2$, a galactic bar. Its pattern speed is ${\Omega }_{\mathrm{p}}$, its amplitude peaks at $t_{\mathrm{peak}}$, and its Gaussian function envelope has duration of order $\tau$. The coefficients $u_{\mathbf{k}}(\mathbf{J})$ describe its strength along different stellar orbits. A real perturbing potential requires $u_{-\mathbf{k}}=u_{\mathbf{k}}^{\ast }$.

Write $g(t)=\exp [-(t-t_{\mathrm{peak}}{)}^2/(2{\tau }^2)]$. The product entering c isThe phase therefore measures detuning from the rotating pattern:Substitution givesThe strongest long-lived response occurs near an orbital resonance, ${\omega }_{\mathbf{k}}=0$, where successive force contributions do not rapidly cancel in phase.

The required time separation is $\tau \ll t_{\sec }$, where $t_{\sec }$ is the evolution timescale of $f_0$, together with $t_{\mathrm{peak}}\gg \tau$ and an observation time well after the pulse. The first condition permits freezing $f_0(\mathbf{J},t^{\prime })$ at its pre-pulse value; its neglected change would enter beyond the retained second order. The second makes the original initial time lie in a negligible tail of the Gaussian function, permitting replacement of the lower time limits by $-\infty$. Short duration alone would not justify that replacement for a pulse centered at the initial time.

Integrate the continuity equation in c over the whole pulse. For fixed $\mathbf{J}$, the time factor is the ordered Gaussian response integralA common translation by $t_{\mathrm{peak}}$ removes the peak time. Pair the positive and negative azimuthal harmonics: $(-m,-k_R)$ is paired with $(m,k_R)$. Both $k_{\phi }$ and $\mathbf{k}\cdot {\partial }_{\mathbf{J}}{f}_0$ change sign, so their product does not; ${\omega }_{-\mathbf{k}}=-{\omega }_{\mathbf{k}}$ makes the two time factors complex conjugates. Hence the pair contributes $2\mathrm{Re}\mathcal{I}$.

For completeness, the real part of the integrand is symmetric under $t\leftrightarrow t^{\prime }$. The triangular time domain therefore gives half of the full plane, and the Fourier transform of a Gaussian function givesThis derives the supplied identity rather than requiring a separate time-ordering calculation.

DefineUsing the paired harmonics and $\Delta F_0=-{\partial }_{J_{\phi }}\int Q_{0}dt$ now yieldsBy the stated assumption, ${\omega }_{m{k}_R}$ is independent of $J_R$, so its exponential can be taken outside the radial integral. Put $\Gamma =(\sqrt{2}\tau {)}^{-1}$ and useThereforeThe outer derivative acts on the Gaussian resonance broadening factor as well as on the radial integral. The result describes quasilinear diffusion in action space weighted by the squared pulse amplitude, with a detuning width $\Gamma$. A shorter pulse has broader frequency selectivity, while a longer pulse selects narrower orbital resonances. If the action-space boundary flux vanishes, integrating $\Delta F_0$ over $J_{\phi }$ gives zero, as required by conservation of mass.

Put $s=z-z_0$. In the source's complex beam parameter convention,Substitution separates the transverse wave amplitude and phase:Thus the Gaussian beam is circularly symmetric, has a Gaussian function transverse light intensity, and narrows to a beam waist at $z_0$ before spreading by diffraction. Its wavefront is planar at the beam waist, converging before it and diverging after it. The quadratic transverse phase identifies the signed radius of curvature:Define spot size as the wave amplitude $1/e$ radius, equivalently the light intensity $1/e^2$ radius. ThenThe Rayleigh range is where $w=\sqrt{2}{w}_0$; the far-field beam divergence half-angle is $\lambda /(\pi w_0)$. The factor $1/q$ also produces a longitudinal phase shift relative to a plane wave, changing by $\pi$ across the beam waist; its sign here follows the given $e^{i(kz-\omega t)}$ convention. The on-axis light intensity decreases as $1/(s^2+z_R^2)$ while transverse integration leaves the power constant. These are scalar paraxial approximation results.

Put $b=\Delta {\omega }_b>0$ and $s=t-t_m$. At the peak time, $e^{-ip\Delta \omega t_m}=1$, so the exact comb sum becomesFor many occupied normal modes, $b\gg \Delta \omega$, replace the sum by an integral with density $1/\Delta \omega$. Setting $x=\omega -{\omega }_0$ and evaluating the Fourier transform of a Gaussian givesThus the isolated-optical pulse approximation hasThe optical pulse has a Gaussian function envelope and no chirp. Its carrier angular frequency is ${\omega }_0$, and its light intensity full width at half maximum is $2\sqrt{2\ln 2}/b$. Successive envelope maxima have different carrier-envelope phases, consistent with b.

Approximation in the printed assertion. A discrete optical frequency comb gives a train, not one exact isolated Gaussian pulse. Applying the Poisson summation formula to the ideal sum over all integers instead givesNear $t_m$, the $j=m$ term dominates when $bT\gg 1$. The other terms are exponentially small but not identically zero. A single Gaussian function used globally would decay between $0$ and $T$, whereas the exact light intensity has equal maxima there; this supplies a concrete reason the isolated form cannot be an exact global identity. The extension to negative frequencies is the usual narrow-positive-band analytic-field approximation when ${\omega }_0\gg b$.

Rearrange as $y^{\prime }+(x-1)y=x-1$ and put $w=y-1$. The resulting linear differential equation is $w^{\prime }+(x-1)w=0$, soThe initial condition gives $C=1$. Its graph is a Gaussian function above the horizontal asymptote $y=1$, symmetric about $x=1$, with its maximum $y=2$ there. It rises for $x<1$ and falls for $x>1$, since $y^{\prime }=-(x-1)e^{-(x-1{)}^2/2}$.

Put $\xi =x+vt$ and $F(\xi )=A{e}^{-{\xi }^2/x_0^2}$, with $x_0\ne 0$. The chain rule givesFor a nonzero pulse the wave equation on a string therefore requiresIt is not a solution for arbitrary $v$; the zero pulse $A=0$ is the trivial exception. Here $\mu$ is linear mass density, and $T>0$ is tension. With $v>0$, the maximum lies at $x=-vt$, so the Gaussian function translates towards decreasing $x$ by $v(t_2-t_1)$ between the two times. Its amplitude and width stay constant.

First alternative. The Euler product of the Riemann zeta function is absolutely convergent for $\sigma >1$. Its logarithm and the nonnegative trigonometric polynomialgiveConsequently $\zeta (\sigma {)}^3|\zeta (\sigma +it){|}^4|\zeta (\sigma +2it)|\ge 1$. Suppose $t\ne 0$ and the Riemann zeta function had a zero of order $m\ge 1$ at $1+it$. As $\sigma \downarrow 1$, its simple pole at one gives $\zeta (\sigma )=O((\sigma -1{)}^{-1})$, the alleged zero gives $\zeta (\sigma +it)=O((\sigma -1{)}^m)$, and $\zeta (\sigma +2it)$ remains bounded. The product would be $O((\sigma -1{)}^{4m-3})\to 0$, contradicting its lower bound. Thus there are no zeta zeros on $\Re s=1$; the point $s=1$ itself is a pole, not a zero.

For the first integral of the Chebyshev function, defineThe logarithmic derivative identity $-{\zeta }^{\prime }/\zeta (s)={\sum }_n\Lambda (n)n^{-s}$, valid for $\Re s>1$, leads toTo verify the kernel rather than assume an inversion formula, integrate $y^s/(s(s+1))$ on the vertical line $\Re s=c$. Closing to the left when $y>1$ gives the residues at zero and minus one, namely $1-y^{-1}$; closing to the right when $0<y<1$ gives zero. At $y=1$ the integral is zero as well. The horizontal integrals tend to zero because of the quadratic denominator, and the distant vertical integrals vanish in the indicated direction. Multiplying this kernel by $x$ with $y=x/n$ gives $(x-n{)}^{+}$. The absolute convergence of the Dirichlet series on $\Re s=c$, together with the integrable denominator, permits interchange with the integral and proves the formula.

Here is a quantitative contour route from this relation to the Prime number theorem. Question 3 establishes a zero-free region of the Riemann zeta function and the accompanying bounds for its logarithmic derivative, independently of the Prime number theorem. In a sufficiently narrow half-width region they implyon the left and horizontal sides of the rectangle with heights $\pm T$ and left edge ${\sigma }_L=1-b/(2\log (T+3))$, after making $b>0$ small enough also for bounded heights. The explanation of this bound is included at the end of Question 3. Take $c=1+1/\log x$ and $T=\exp (\sqrt{\log x})$. On the original line, the omitted tails aresince the Euler product bounds $|{\zeta }^{\prime }/\zeta (c+it)|$ by $-{\zeta }^{\prime }/\zeta (c)=O(\log x)$. Shift the truncated integral to the left edge. The only enclosed pole is at $s=1$, with residue $x^2/2$; the kernel poles at zero and minus one lie outside the rectangle. The left edge contributes $O(x^{2-b/(2\log (T+3))}{\log }^2(T+3))$, since $|s(s+1){|}^{-1}\ll (1+t^2{)}^{-1}$ there. The two horizontal edges contribute $O(x^2{\log }^2(T+3)/T^2)$. In particular, for some $a>0$,

It remains to remove the smoothing; differentiating an asymptotic formula would not justify this step. The monotonicity of the Second Chebyshev function gives, for $h=\epsilon x$ with $0<\epsilon <1$,First let $x\to \infty$ with $\epsilon$ fixed. After division by $x$, the lower and upper bounds tend to $1-\epsilon /2$ and $1+\epsilon /2$. Let $\epsilon \downarrow 0$ to obtain $\psi (x)\sim x$. Higher prime powers contribute $O(\sqrt{x}\log x)$ by the Chebyshev estimate, so $\theta (x)\sim x$. Finally partial summation givesThe integral is $O(x/{\log }^{2}x)$, by splitting at $\sqrt{x}$ and using $\theta (t)=O(t)$. This proves the Prime number theorem with $\pi (x)\sim x/\log x$.

Second alternative. We prove the functional equation of the Riemann zeta function from a theta function. Set $\Theta (u)={\sum }_{n\in \mathbb{Z}}{e}^{-\pi n^{2}u}$ for $u>0$. With Fourier transform convention $\hat{f}(\xi )={\int }_{\mathbb{R}}f(v)e^{-2\pi iv\xi }dv$, the Gaussian function $f(v)=e^{-\pi u{v}^2}$ hasIndeed, differentiation under the integral and integration by parts give ${\hat{f}}^{\prime }(\xi )=-2\pi \xi \hat{f}(\xi )/u$, while the Gaussian integral gives $\hat{f}(0)=u^{-1/2}$. Periodize $f$: the smooth period-one function ${\sum }_{n}f(v+n)$ has Fourier coefficients $\hat{f}(k)$, by integrating on $[0,1]$ and interchanging its absolutely convergent sum with the integral. Its absolutely convergent Fourier series, evaluated at zero, proves the Poisson summation formula in this case. HenceFor $\Re s>1$, termwise integration and the defining integral of the Gamma function give the Mellin representation of the completed Riemann zeta function:Split at one and substitute $u=1/v$ in the lower integral. The elementary part isand the remaining part becomes an integral over $[1,\infty )$. ThusThe integral is an entire function of $s$, because $\Theta (u)-1$ decays exponentially and uniform convergence holds on every compact subset of the complex plane. This proves meromorphic continuation of the completed Riemann zeta function, with only the displayed simple poles, and the right side is invariant under $s\mapsto 1-s$. ThereforeThe Gamma reflection formula and Gamma duplication formula give the equivalent form $\zeta (s)=2^s{\pi }^{s-1}\sin (\pi s/2)\Gamma (1-s)\zeta (1-s)$, understood by meromorphic continuation. Multiplication by $s(s-1)/2$ gives the entire Riemann xi function $\xi (s)=\xi (1-s)$.

Let $N(T)$ count positive ordinates of Nontrivial zeros of the Riemann zeta function up to height $T$, including multiplicity. The Riemann–von Mangoldt formula isMore precisely, away from ordinates of zeros, the constant term is $7/8$, the additional term is $S(T)={\pi }^{-1}\arg \zeta (1/2+iT)$, and the remaining error is $O(T^{-1})$. Here the argument is continued along $2\to 2+iT\to 1/2+iT$, and $S(T)=O(\log T)$. The coarse formula displayed above holds at all large heights with the usual right-continuous counting convention. Subtract it at $T+1$ and $T$: the main term changes by $O(\log T)$ and each error is $O(\log T)$. Hence

For infinitely many reciprocal-logarithmic gaps between zeta zeros, fix any $0<c<2\pi$, for example $c=\pi$. Suppose all sufficiently late gaps satisfied ${\gamma }_{n+1}-{\gamma }_n\le c/\log n$. Summing and usinggives ${\gamma }_n\le (c+o(1))n/\log n$. The last elementary estimate follows by integral comparison with $1/\log t$ and one integration by parts. The Riemann–von Mangoldt formula would implya contradiction. This reasoning is valid whether repeated ordinates are listed according to multiplicity or only distinct ordinates are listed, since in either case $n\le N({\gamma }_n)$. Consequently ${\gamma }_{n+1}-{\gamma }_n>c/\log n$ infinitely often, for any fixed $0<c<2\pi$.

Repeated differentiation of the Gaussian function gives $D^n{e}^{-x^2}=P_n(x)e^{-x^2}$, where $P_0=1$ and $P_{n+1}=P_n^{\prime }-2x{P}_n$. Induction shows $P_n$ has degree of a polynomial $n$ and leading coefficient $(-2{)}^n$: the differentiated term has lower degree, while multiplication by $-2x$ raises it by one. Thus the Hermite polynomial $H_n=(-1{)}^n{P}_n$ has degree $n$ and leading coefficient $2^n$.

For $m<n$, apply integration by parts $n$ times to the Rodrigues expression:All boundary terms vanish because each derivative of $e^{-x^2}$ is a polynomial times that exponential, and such products tend to zero at either infinity. The same result for $m>n$ follows by symmetry of the inner product. The diagonal value is positive; indeed $H_n^{(n)}=2^{n}n!$ givesThis proves the requested orthogonality as well as its normalization.

Use the Fourier transform convention $\hat{h}(\xi )={\int }_{\mathbb{R}}h(x)e^{-2\pi ix\xi }dx$. For a Schwartz function $h$, the Poisson summation formula isTo prove it, periodize $h$ as $P(x)={\sum }_{n\in \mathbb{Z}}h(x+n)$. Rapid decay of all derivatives makes this a smooth periodic function, with uniform termwise differentiability. Its $m$th Fourier coefficient iswhere absolute convergence justifies splitting the real-line integral into unit intervals. Integration by parts gives rapid decay of the coefficients of $P$, so its Fourier series converges absolutely and uniformly to $P$. Evaluate it at zero to obtain the formula.

For $\mathrm{Im}\tau >0$ and $z\in \mathbb{C}$, take the complex Gaussian function $h(x)=e^{\pi i\tau x^2+2\pi izx}$. It is a Schwartz function of the real variable $x$. Its Fourier transform isThe real Gaussian integral proves this first for $\tau =iy$, $y>0$, and real $z$. Normal convergence of the defining integrals makes both sides entire in $z$ and holomorphic in $\tau \in \mathbb{H}$, so the identity theorem extends it to the stated domain. Choose the square-root branch holomorphic for $-i\tau$ in the right half-plane and positive when $\tau =iy$.

Applying Poisson summation gives the transformation of the Jacobi theta function:The series on the right is ${\vartheta }_{00}(z/\tau ,-1/\tau )$. Rearranging therefore givesThe branch choice is part of the formula; at $\tau =i$ its square-root factor equals one.

Use the angular-frequency Fourier transform conventionA precise sufficient interpretation of the regularity assumption is $f\ne 0$, $f\in H^1(\mathbb{R})$ in the first-order Sobolev space and $xf\in L^2(\mathbb{R})$; in particular, nonzero Schwartz functions are admissible. We use the Plancherel theorem, in the normalization $\Vert \hat{f}{\Vert }_2^2=2\pi \Vert f{\Vert }_2^2$, together with its derivative identity $\widehat{f^{\prime }}(\lambda )=i\lambda \hat{f}(\lambda )$. Both identities hold in $L^2$ for $f\in H^1$, by extending their identities for Schwartz functions. Thus the frequency factor isFor a Schwartz function, integration by parts gives $\Vert f{\Vert }_2^2=-2\mathrm{Re}\int x{f}^{\prime }(x)\overline{f(x)}dx$. This remains true under the stated assumptions: insert a smooth cutoff ${\chi }_R(x)=\chi (x/R)$ into the derivative of $x|f{|}^2$, where $\chi$ is one near zero and compactly supported. Integration givesThe cutoff-derivative term tends to zero, since $x{\chi }_R^{\prime }$ is bounded independently of $R$ and supported in a tail. The other terms converge by dominated convergence and the integrability of $|f^{\prime }||xf|$, furnished by the Cauchy-Schwarz inequality. HenceAfter squaring and using the Plancherel theorem, we obtain the uncentred Fourier uncertainty principle:If equality holds, both displayed inequalities must be equalities. The equality condition in the Cauchy-Schwarz inequality gives $f^{\prime }=cxf$ almost everywhere for a constant $c\in \mathbb{C}$; the real-part inequality and the identity for $\Vert f{\Vert }_2^2$ then force $c=-a$ with $a>0$ real. The weak derivative equation has the solution $f(x)=C{e}^{-a{x}^2/2}$: locally multiply by $e^{a{x}^2/2}$ to obtain a function whose weak derivative is zero, hence a constant. Conversely, for this Gaussian function the normalized second moments are $1/(2a)$ in position and $a/2$ in angular frequency, whose product is $1/4$. ThereforeThese are raw second moments, rather than variances about their respective means. Consequently translations and nonzero frequency modulations of the Gaussian function are not equality cases here; the equality case of the Heisenberg uncertainty relation for centred variances has those additional freedoms.

Retain $\Delta =\mathrm{div}\nabla$ and define the heat operator by $L={\partial }_t-{\Delta }_q$. A Riemannian heat kernel $H(p,q,t)$ satisfies $L_{q}H=0$ for $t>0$ and has the delta initial limitOn standard Euclidean space the kernel isOn the Riemannian manifold, the Gaussian function factor in the local construction uses $r(p,q)=\mathrm{dist}(p,q)$ in place of $|p-q|$.

Choose an open neighbourhood $\mathcal{U}$ of the diagonal such that every pair $(p,q)\in \mathcal{U}$ is joined by a unique short minimizing geodesic and $q={\exp }_{p}x$ with the inverse exponential map smooth in both variables. It can be obtained from normal coordinates tubes of radius smaller than the local injectivity radius, or from normal neighbourhoods with geodesic convexity. Then $r^2$ is smooth on $\mathcal{U}$. Write the Riemannian volume form in normal coordinates as $d{V}_q=J(p,x)dx$, where $J>0$ is smooth and $J(p,0)=1$.

For fixed $p$, radial divergence in normal coordinates givesFor $g=(4\pi t{)}^{-d/2}{e}^{-r^2/(4t)}$, use $|\nabla r|=1$ away from the diagonal and the product rule to obtainThe identities extend smoothly in the expressions involving $r{\partial }_r$, even though $r$ itself is not smooth on the diagonal. The singular powers cancel ifThe correct initial normalization is $w_0(p,p)=1$. These heat-kernel transport equations have the explicit solutionsand, recursively,To derive the integral, multiply the radial equation by $r^{i-1}{J}^{1/2}$ and integrate the derivative of $r^i{J}^{1/2}{w}_i$ from zero to $r$. The zero integration constant is forced by smoothness when $i\ge 1$. Changing variables from radial distance to $s\in [0,1]$ gives the displayed nonsingular expression. Its integrand is smooth in $(p,x,s)$, so differentiation under the integral proves inductively that every $w_i$ is smooth across the diagonal. In particular $w_i(p,p)=({\Delta }_q{w}_{i-1})(p,p)/i$.