How to Renormalize The Schrödinger Equation

How to Renormalize The Schrödinger Equation#

Here we work through the example discussed in [Lepage, 1997], which uses the example of bound states in a spherically symmetric potential, with Coulomb-like behaviour for large radii \(r\rightarrow \infty\), but with “unknown” corrections at short distances \(r\rightarrow 0\):

\[\begin{gather*} \left(\frac{-\hbar^2\nabla^2}{2m} + V(r) - E\right)\Psi(r, \Omega) = 0, \qquad \lim_{r\rightarrow \infty} V(r) \rightarrow \frac{-\alpha}{r}. \end{gather*}\]

How to Solve the Schrödinger Equation#

To easily work through the discussion in [Lepage, 1997], we must be able to formulate and solve the Schrödinger equation fr spherically symmetric potentials.

Radial Equation#

Spherical symmetry allows use to express this as a simple 1D boundary value problem (BVP) for the radial equation:

With the usual spherical coordinates \(\phi\) and \(\theta\), and angular momentum quantum numbers \(l\) and \(m\), we have

\(d\)	\(\Omega\)	\(\lambda\)
\(2\)	\(\phi\)	\(\abs{m}\)
\(3\)	\((\phi, \theta)\)	\(l\)

for \(d=2\) dimensions (i.e. cylindrical coordinates, but no \(z\) dependence) and \(d=3\) dimensions respectively. The total orbital angular momentum operator is related to the Laplace-Beltrami operator by

\[\begin{gather*} \op{L}^2 = \hbar^2\op{\Delta}_{S^{d-1}} \end{gather*}\]

which has eigenvalues \(\hbar^2 m^2\) in \(d=2\) (i.e. \(\op{L}^2 \equiv \op{L}_z^2\) and \(\hbar^2 l (l - 1)\) in \(d=3\).

Our aim is to satisfy the Schrödinger equation for central potentials in \(d\)-dimensions, which we can do in the usual way by expressing the wavefunction \(\Psi(r, \Omega) = \psi(r)Y_{\lambda}(\Omega)\) in terms of the radial wavefunction \(u(r)\) and the appropriate generalized spherical harmonics \(Y_{\lambda}(\Omega)\):

\[\begin{gather*} \left(\frac{-\hbar^2\nabla^2}{2m} + V(r) - E\right)\Psi(r, \Omega) = 0, \\ \left[ \frac{\hbar^2}{2m} \underbrace{ \left(-\diff[2]{}{r} + \frac{\nu^2 - 1/4}{r^2}\right) }_{\op{K}} + V(r)\right]u(r) = E u(r),\\ u(r) = r^{(d-1)/2}\psi(r), \qquad \nu = \lambda + \frac{d}{2} - 1. \end{gather*}\]

Details

Here \(\Omega\) is the generalized solid angle. Rotational invariance implies that angular momentum is a good quantum number, so all eigenstates can be factored \(\Psi(r, \Omega) = \psi(r)Y_{\lambda}(\Omega)\) where \(Y_{\lambda}(\Omega)\) is the generalized spherical harmonic on the (\(d-1\))-dimensional sphere and \(\lambda \in \{0, 1, \dots\}\) is the generalized angular momentum:

\[\begin{gather*} \nabla^2 = \frac{1}{r^{d-1}}\diff{}{r}\left(r^{d-1} \diff{}{r}\right) + \frac{1}{r^2}\Delta_{S^{d-1}}, \\ \Delta_{S^{d-1}}Y_{\lambda}(\Omega) = \lambda (\lambda + d - 2)Y_{\lambda}(\Omega), \end{gather*}\]

where \(\Delta_{S^{d-1}}\) is the Laplace-Beltrami operator. Introducing the radial wavefunction \(u(r)\), this becomes quadratic:

\[\begin{gather*} u(r) = r^{(d-1)/2}\psi(r), \qquad \nu = \lambda + \frac{d}{2} - 1,\\ r^{(d-1)/2}\nabla^2 \psi(r) = \left(\diff[2]{}{r} - \frac{\nu^2 - 1/4}{r^2}\right)u(r). \end{gather*}\]

This follows after a little algebra from

\[\begin{gather*} \frac{1}{r^{d-1}}\diff{}{r}\left(r^{d-1} \diff{}{r}\right)r^{(1-d)/2}u(r)\\ = r^{(1-d)/2} \left( \diff[2]{}{r} - \frac{(d-3)(d-1)}{4} \right)u(r). \end{gather*}\]

This can be solved quite simply – but not very accurately – with finite differences. Highly accurate solutions can be obtained by shooting, and we present details in How to Solve the (Radial) Schrödinger Equation about how to do this.

The Essence#

The essential idea is low-energy properties, like the bound state energies \(E_n\), should not be highly sensitive to details about the nature of the potential \(V(r)\) at short distances \(r \ll \lambda\)

Numerical Approach#

The technical problem is to compute properties such as the bound state energies \(E_n\) of the potential \(V(r)\). [Lepage, 1997] presents a nice derivation about how to do this using perturbation theory with various contact interaction terms, and I urge you to follow and reproduce this discussion.

Here we will take an alternative numerical approach, defining our own effective potential as follows:

\[\begin{gather*} V(r) = P(r/z)f_a(r) - \frac{\alpha}{r}\bigl(1-f_a(r)\bigr), \qquad f_a(r) = e^{-r^2/2a^2}, \qquad P(z) = \sum_{n} c_n \frac{z^n}{n!}. \end{gather*}\]

This satisfies the criteria laid out in [Lepage, 1997]:

We incorporate the correct long-range behavior through the cutoff function \(1-f_a(r)\) which goes to zero exponentially fast for \(r>a\).
We have introduced an ultraviolet cutoff \(a\) into our theory which softens the potential at short distances.
We have added “local” corrections via the parameters \(c_n\). The locality is provided by the cutoff factor \(f_a(r)\).

import warnings;warnings.simplefilter("error")
from functools import partial
from scipy.optimize import root
from tqdm import tqdm
from phys_552 import seq

class SEQ(seq.CoulombSEQ):
    c = [-2.0]
    a = 1.0
    
    E_tol = 1e-4
    
    def f(self, r):
        return np.exp(-(r/self.a)**2/2)
        
    def V(self, r):
        f = self.f(r)
        return np.polyval(self.c, r/self.a)*f + (1-f)*super().V(r)
    
    def fit(self, Es, cs=None, **kw):
        """Fit the specified energies."""
        if cs is None:
            cs = self.c
        cs = list(cs)[:len(Es)]
        cs = cs + [0.0]*(len(Es) - len(cs))
        self.c = cs
        
        def objective(cs, Es):
            self.c[:len(Es)] = cs
            Es_ = np.array([self.compute_E(_E, tol=self.E_tol, lam=0.999) for _E in Es])
            return Es_/Es - 1
        
        for n in tqdm(range(1, len(Es)+1)):
            f = partial(objective, Es=Es[:n])
            res = root(f, self.c[:n], **kw)
        self.res = res
        if not res.success:
            raise Exception(res.message)
        self.c[:len(Es)] = res.x

        
s = SEQ()

Es = np.array([
    -1.28711542, 
    -0.183325753,
    -0.0703755485,
    -0.0371495726,
    -0.0229268241,
    -0.0155492598,
    -0.00534541931,
    -0.00129205010])

---------------------------------------------------------------------------
ModuleNotFoundError                       Traceback (most recent call last)
Cell In[2], line 4
      2 from functools import partial
      3 from scipy.optimize import root
----> 4 from tqdm import tqdm
      5 from phys_552 import seq
      7 class SEQ(seq.CoulombSEQ):

ModuleNotFoundError: No module named 'tqdm'

#s.compute_E(-1)
s.fit(Es[:2])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[3], line 2
      1 #s.compute_E(-1)
----> 2 s.fit(Es[:2])

NameError: name 's' is not defined

%time s.fit(Es[:2])

CPU times: user 4 μs, sys: 1 μs, total: 5 μs
Wall time: 8.34 μs

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[4], line 1
----> 1 get_ipython().run_line_magic('time', 's.fit(Es[:2])')

File ~/checkouts/readthedocs.org/user_builds/physics-552-quantum-iii/conda/latest/lib/python3.13/site-packages/IPython/core/interactiveshell.py:2532, in InteractiveShell.run_line_magic(self, magic_name, line, _stack_depth)
   2530     kwargs['local_ns'] = self.get_local_scope(stack_depth)
   2531 with self.builtin_trap:
-> 2532     result = fn(*args, **kwargs)
   2534 # The code below prevents the output from being displayed
   2535 # when using magics with decorator @output_can_be_silenced
   2536 # when the last Python token in the expression is a ';'.
   2537 if getattr(fn, magic.MAGIC_OUTPUT_CAN_BE_SILENCED, False):

File ~/checkouts/readthedocs.org/user_builds/physics-552-quantum-iii/conda/latest/lib/python3.13/site-packages/IPython/core/magics/execution.py:1448, in ExecutionMagics.time(self, line, cell, local_ns)
   1446 if interrupt_occured:
   1447     if exit_on_interrupt and captured_exception:
-> 1448         raise captured_exception
   1449     return
   1450 return out

File ~/checkouts/readthedocs.org/user_builds/physics-552-quantum-iii/conda/latest/lib/python3.13/site-packages/IPython/core/magics/execution.py:1398, in ExecutionMagics.time(self, line, cell, local_ns)
   1396 st = clock2()
   1397 try:
-> 1398     out = eval(code, glob, local_ns)
   1399 except KeyboardInterrupt as e:
   1400     captured_exception = e

File <timed eval>:1

NameError: name 's' is not defined

Figure 1#

Here we reproduce Fig. 1 from [Lepage, 1997] using our potential:

s2 = SEQ(a=1, E_tol=1e-4)
s2.fit(Es[:2])
Es2 = np.array([s2.compute_E(_E, lam=0.999) for _E in tqdm(Es)])

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[6], line 1
----> 1 s2 = SEQ(a=1, E_tol=1e-4)
      2 s2.fit(Es[:2])
      3 Es2 = np.array([s2.compute_E(_E, lam=0.999) for _E in tqdm(Es)])

NameError: name 'SEQ' is not defined

s2.compute_E(Es[0]) - Es[0]

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[7], line 1
----> 1 s2.compute_E(Es[0]) - Es[0]

NameError: name 's2' is not defined

plt.loglog(abs(Es), abs(Es2/Es - 1))

---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
Cell In[8], line 1
----> 1 plt.loglog(abs(Es), abs(Es2/Es - 1))

NameError: name 'Es' is not defined