Algebraic DVR Approaches Applied to Describe the Stark Effect

Abstract

Two algebraic approaches based on a discrete variable representation are introduced and applied to describe the Stark effect in the non-relativistic Hydrogen atom. One approach consists of considering an algebraic representation of a cutoff 3D harmonic oscillator where the matrix representation of the operators $$r^2$$ and $$p^2$$ are diagonalized to define useful bases to obtain the matrix representation of the Hamiltonian in a simple form in terms of diagonal matrices. The second approach is based on the $$U\left(4\right)$$ dynamical algebra which consists of the addition of a scalar boson to the 3D harmonic oscillator space keeping constant the total number of bosons. This allows the kets associated with the different subgroup chains to be linked to energy, coordinate and momentum representations, whose involved branching rules define the discrete variable representation. Both methods, although originating from the harmonic oscillator basis, provide different convergence tests due to the fact that the associated discrete bases turn out to be different. These approaches provide powerful tools to obtain the matrix representation of 3D general Hamiltonians in a simple form. In particular, the Hydrogen atom interacting with a static electric field is described. To accomplish this task, the diagonalization of the exact matrix representation of the Hamiltonian is carried out. Particular attention is paid to the subspaces associated with the quantum numbers $$n=2,3$$ with $$m=0$$.

1. Introduction

A discrete variable representation approach is based on the search of a discrete basis in terms of which any function of the coordinates is diagonal. The discrete variable representation methods (DVR) in configuration space were developed with some variants since the 1960s, but systematically widely used during the 1980s with different names: discrete-variable representation method [1,2,3,4], quadrature discretization method [5,6], configuration localized states (CLS) approach [7], and Lagrange mesh method (LMM) [8,9,10,11], whose similarities and differences are discussed in Ref. [11]. The basic ingredient of these methods is the use of orthogonal polynomials and their associated grids from Gaussian quadratures.

In this contribution, we introduce two algebraic DVR methods where the discrete bases are obtained using purely algebraic means without any explicit reference to polynomials. This is accomplished through the diagonalization of the matrix representations associated with the coordinates and momenta. From this perspective, the recently-proposed unitary group approach ($U\left(4\right)$-UGA) belongs to an algebraic DVR method [12,13,14,15,16]. In this approach, the discretization provided by the zeros of the polynomials in configuration space is substituted by the branching rules involved in the subgroup chains embedded in the dynamical group. An alternative approach to establish an algebraic DVR method consists of taking a cutoff 3D harmonic oscillator in Fock space and diagonalize the matrix representation of the squares of the coordinates and momenta to obtain two discrete bases (HO-DVR approach), which in turn are used to diagonalize the matrix representation of the Hamiltonian in a simple manner. Both methods are based on a harmonic oscillator basis, albeit they provide different convergence behavior due to the fact that the provided discrete bases are distinct.

In this contribution we present the HO-DVR approach as well as the salient features of the $U\left(4\right)$-UGA, both applied to the Coulombic potential. We first show the confidence of our approaches to reproduce the analytical wave functions for the non relativistic Hydrogen atom. Both approaches are able to arrive at the same level of accuracy with the appropriate basis dimension. In order to explore the capabilities and differences between these approaches, the Stark effect is studied by adding the potential energy for the electron in an homogeneous electric field. Our treatment involves a complete dipole matrix representation contemplating the possibility of the action of fields beyond the perturbation region. Since the conservation of the angular momentum is broken, the calculations become particularly heavy. This fact forces us to constraint the basis to a limited angular momentum. With the proposal of these DVR algebraic approaches, we do not intend to improve the variational approaches used to describe the Stark effect, but rather for showing a simple way to deal with this problem as well as to evaluate the differences between them.

Our contribution is organized as follows. In Section 2 the salient features of both the algebraic HO-DVR method and the $U\left(4\right)$-UGA are presented. Section 3 is devoted to applying the methods to describe the isolated non-relativistic Hydrogen atom. To this end, both energies and wave functions are tested. In Section 4 the Stark effect is described, paying attention to the bound states originated from the subspaces with $n=2,3$. Finally, in Section 5 the summary and our conclusions are presented.

2. Algebraic DVR Methods in 1D Systems

In this section, we present the salient features of the HO-DVR method as well as of the $U\left(4\right)$-UGA. Since both approaches are based on the 3D harmonic oscillator we start presenting the main features associated with its algebraic solutions.

2.1. 3D Harmonic Oscillator

We start considering the isotropic 3D harmonic oscillator with reduced mass m and frequency $\omega$. The Hamiltonian is given by

(1)${\widehat{H}}_{\mathrm{CS}}^{\mathrm{HO}}=\frac{1}{2m}{p}^2+\frac{m{\omega }^2}{2}{r}^2.$

The solutions in configuration space are well known, but this is not the case in Fock space [17,18]. The Hamiltonian, as well as the eigenfunctions can be translated into the Fock space through the introduction of the bosonic creation $a_{\mu }^{†}$ and annihilation $a_{\mu }$ operators with commutation relations

(2)$[a_{\mu },a_{\nu }^{†}]={\delta }_{\mu \nu };[a_{\mu }^{†},a_{\nu }^{†}]=[a_{\mu },a_{\nu }]=0.$

The reflection and rotation properties of $a_{\mu }^{†}$ are summarized by their transformation under a rotation $\widehat{R}(\alpha ,\beta ,\gamma )$:

(3)${\widehat{R}}^{-1}{a}_{\mu }^{†}\widehat{R}=\sum _{{\mu }^{\prime }}{D}_{{\mu }^{\prime }\mu }^{\left(1\right)}(\alpha ,\beta ,\gamma )a_{{\mu }^{\prime }}^{†},$

characterized by the three Euler angles $(\alpha ,\beta ,\gamma )$. The annihilation operators $a_{\mu }$, however, do not transform as the components of a spherical tensor operator, but the modified operator

(4)${\tilde{a}}_{\mu }={(-1)}^{1-\mu }{a}_{-\mu }$

do satisfy the transformation property (3). In this way the the operators $a_{\mu }^{†}$ and ${\tilde{a}}_{\mu }$ can be coupled to definite angular momentum $\lambda$ and projection m

(5)${[a^{†}\times \tilde{a}]}_m^{\left(\lambda \right)}=\sum _{{\mu }^{\prime },\mu }⟨1\mu 1{\mu }^{\prime }|\lambda m⟩a_{\mu }^{†}{\tilde{a}}_{{\mu }^{\prime }},$

using the standard notations × for the angular momentum coupling and the Clebsh–Gordan coefficients. These operators in terms of the spherical coordinate and momentum components are given by

(6)$\begin{array}{cc}\hfill a_{\mu }^{†}=& \frac{1}{\sqrt{2}}\left(\frac{1}{d}{q}_{\mu }-i\frac{d}{\hslash }{(-1)}^{\mu }{p}_{-\mu }\right);\hfill \end{array}$

(7)$\begin{array}{cc}\hfill {\tilde{a}}_{\mu }=& -\frac{1}{\sqrt{2}}\left(\frac{1}{d}{q}_{\mu }+i\frac{d}{\hslash }{(-1)}^{\mu }{p}_{-\mu }\right),\hfill \end{array}$

with commutation relations

(8)$[{\tilde{a}}_{\mu },a_{\nu }^{†}]={(-1)}^{1-\mu }{\delta }_{-\mu \nu };[a_{\mu }^{†},a_{\nu }^{†}]=[{\tilde{a}}_{\mu },{\tilde{a}}_{\nu }]=0,$

where $d=\sqrt{\hslash /m\omega }$ has been introduced as unit length, $q_{\mu }(\mu =1,0,-1)$ are the coordinate components with $r^2={\sum }_{\mu }{(-1)}^{\mu }{q}_{\mu }{q}_{-\mu }$ and corresponding momenta $p_{\mu }=-i\hslash \frac{\partial }{\partial q_{\mu }}$. The solutions of the eigensystem associated with (1) in Fock space are given by the kets [17,18]

(9)$|nlm⟩=B_{nl}{({\mathbf{a}}^{†}·{\mathbf{a}}^{†})}^{(n-l)/2}{\mathcal{Y}}_m^l\left({\mathbf{a}}^{†}\right)|0⟩,$

characterized by the eigensystem

(10)$\begin{array}{cc}\hfill \widehat{n}|nlm⟩=& n|nlm⟩,\hfill \end{array}$

(11)$\begin{array}{cc}\hfill {\widehat{L}}^2|nlm⟩=& l(l+1)|nlm⟩,\hfill \end{array}$

(12)$\begin{array}{cc}\hfill {\widehat{L}}_z|nlm⟩=& m|nlm⟩,\hfill \end{array}$

where ${\mathcal{Y}}_m^l\left({\mathbf{a}}^{†}\right)$ are the solid spherical harmonics defined by

(13)${\mathcal{Y}}_m^l\left({\mathbf{a}}^{†}\right)=2^{-l/2}{({\mathbf{a}}^{†}·{\mathbf{a}}^{†})}^{l/2}{Y}_m^l\left({\mathbf{a}}^{†}\right)$

in terms of the spherical harmonics $Y_m^l\left({\mathbf{a}}^{†}\right)$ with normalization [17,18]

(14)$B_{nl}={(-1)}^{(n-l)/2}\sqrt{\frac{4\pi }{(n-l)!!(n+l+1)!!}},$

and eigenvalues

(15)$E_n^{\mathrm{HO}}=\hslash \omega (n+3/2),$

with frequency $\omega =\sqrt{k/m}$ in terms of the force constant k and number operator

(16)$\widehat{n}=\sum _{\mu }{a}_{\mu }^{†}{a}_{\mu }=\sqrt{3}[{\mathbf{a}}^{†}\times \tilde{\mathbf{a}}{]}_{\mathbf{0}}^{\left(\mathbf{0}\right)},$

while the notation ${\mathbf{a}}^{†}·{\mathbf{a}}^{†}$ stands for the dot product defined by

(17)${\mathbf{a}}^{†}·{\mathbf{a}}^{†}=\sqrt{3}{[{\mathbf{a}}^{†}\times {\mathbf{a}}^{†}]}_0^{\left(0\right)}=\sum _{\mu }{(-1)}^{\mu }{a}_{\mu }^{†}{a}_{-\mu }^{†},$

with the realization of the angular momentum components

(18)${\widehat{L}}_{\mu }=\sqrt{2}{[a^{†}\times \tilde{a}]}_{\mu }^{\left(1\right)}.$

This summary for the harmonic oscillator paves the way to introduce the algebraic DVR methods.

2.2. HO-DVR Approach

The HO-DVR approach consists of establishing discrete bases where functions of the radial coordinate and the corresponding momentum are diagonal. To accomplish this task we should consider the algebraic realization of the variables involved in Equation (1). Based on the relations (6) and (7), the squared of the radius and momenta in Fock space takes the form

(19)$\begin{array}{cc}\hfill r^2& =(\widehat{n}+3/2)+\frac{1}{2}({\mathbf{a}}^{†}·{\mathbf{a}}^{†}+\mathbf{a}·\mathbf{a}),\hfill \end{array}$

(20)$\begin{array}{cc}\hfill p^2& =(\widehat{n}+3/2)-\frac{1}{2}({\mathbf{a}}^{†}·{\mathbf{a}}^{†}+\mathbf{a}·\mathbf{a}),\hfill \end{array}$

Since both operators (19) and (20) preserve the angular momentum we can obtain their matrix representation for a given l and m in the N-dimensional space defined by

(21)${\mathcal{L}}_N^l=\left\{\right|nlm⟩,n=0,1,\dots ,N-1\},$

with the restrictions $l=n,n-2,\dots ,1\mathrm{or}0$, and the usual reduction $SO\left(3\right)\supset SO\left(2\right)$ given by $-l\le m\le l$. It is clear that the N-dimensional space (21) is not complete, but it establishes a space to carry our the calculations in an approximate way. Explicitly the matrix elements of the operators (19) and (20) in dimensionless units turn out to be

(22)$\begin{array}{cc}\hfill ⟨nlm\left|r^2\right|nlm⟩& =⟨nlm\left|p^2\right|nlm⟩,\hfill \\ \hfill & =n+\frac{3}{2},\hfill \end{array}$

(23)$\begin{array}{cc}\hfill ⟨n+2,lm\left|r^2\right|nlm⟩& =-⟨n+2,lm\left|p^2\right|nlm⟩\hfill \\ \hfill & =\frac{1}{2}\sqrt{(n+l+2)(n-l+2)}.\hfill \end{array}$

The numerical diagonalizations of the representation matrix of these operators provide a discrete set of eigenvectors

(24)$\begin{array}{cc}\hfill |lm,r_i^2⟩=& \sum _{n=0}^{N-1}⟨nlm|lm,r_i^2⟩|nlm⟩,\hfill \end{array}$

(25)$\begin{array}{cc}\hfill |lm,p_i^2⟩=& \sum _{n=0}^{N-1}⟨nlm|lm,p_i^2⟩|nlm⟩,\hfill \end{array}$

which provide the radial coordinate and momentum representations characterized by

(26)$\begin{array}{cc}\hfill r|lm,r_i^2⟩=& \sqrt{r_i^2}|lm,r_i^2⟩,\hfill \end{array}$

(27)$\begin{array}{cc}\hfill p|lm,p_i^2⟩=& \sqrt{p_i^2}|lm,p_i^2⟩.\hfill \end{array}$

Because of these properties, the basis $|lm,r_i^2⟩$ is identified with the position representation, while $|lm,p_i^2⟩$ stands for the momentum representation. In addition since the Hamiltonian (15) is diagonal in the basis (9) we refer the $U\left(3\right)$ basis as the energy representation. As a consequence, for any functions $V\left(r\right)$ and $G\left(p\right)$ we have for their matrix elements

(28)$\begin{array}{cc}\hfill ⟨lm,r_j^2\left|V\left(r\right)\right|lm,r_i^2⟩& =V\left(\sqrt{r_i^2}\right){\delta }_{ij},\hfill \end{array}$

(29)$\begin{array}{cc}\hfill ⟨lm,p_j^2\left|G\left(p\right)\right|lm,p_i^2⟩& =G\left(\sqrt{p_i^2}\right){\delta }_{ij}.\hfill \end{array}$

The eigenvectors (24) and (25) define the transformation coefficients

(30)$\begin{array}{cc}\hfill \mathbf{W}& =\left|\right|⟨nlm|lm,p_i^2⟩||,\hfill \end{array}$

(31)$\begin{array}{cc}\hfill \mathbf{T}& =\left|\right|⟨nlm|lm,r_i^2⟩||,\hfill \end{array}$

while the properties (28) and (29) lead to a DVR method where the matrix representation of the general 3D Hamiltonian depending on the radial coordinate

(32)$\widehat{H}=\frac{p^2}{2m}+V\left(r\right),$

takes the following matrix form in the energy representation:

(33)$\mathbf{H}={\mathbf{W}}^{†}{\mathsf{\Lambda }}^{\left(p\right)}\mathbf{W}+{\mathbf{T}}^{†}{\mathsf{\Lambda }}^{\left(r\right)}\mathbf{T},$

where ${\mathsf{\Lambda }}^{\left(p\right)}$ is the diagonal contribution of the kinetic energy in the momentum representation (25):

(34)$\left|\right|{\mathsf{\Lambda }}^{\left(p\right)}\left|\right|=\left|\right|⟨lm,p_j^2|\frac{p^2}{2m}|lm,p_i^2⟩\left|\right|=\frac{p_i^2}{2m}{\delta }_{ij},$

where the notation for the mass m should not be confused with the angular momentum projection m. In addition the matrix ${\mathsf{\Lambda }}^{\left(r\right)}$ corresponds to the diagonal matrix of the potential $V\left(r\right)$ in the position representation:

(35)$\left|\right|{\mathsf{\Lambda }}^{\left(r\right)}\left|\right|=\left|\right|⟨lm,r_j^2|\left[V\left(\sqrt{r^2}\right)\right]|lm,r_i^2⟩\left|\right|=\left[V\left(\sqrt{r_i^2}\right)\right]{\delta }_{ij}.$

Both transformation coefficients $\mathbf{W}$ and $\mathbf{T}$ are numerically calculated only once for a given l and m and space dimension N, providing the methodology for any potential in the framework of the HO-DVR method. However since the basis corresponds to the harmonic oscillator functions, adding and subtracting the quadratic potential makes possible to reformulate the Equation (33) in terms of only the transformation coefficients $\mathbf{T}$ in the form:

(36)$\mathbf{H}={\mathsf{\Lambda }}^{\left(E\right)}+\kappa {\mathbf{T}}^{†}{{\mathsf{\Lambda }}^{\prime }}^{\left(r\right)}\mathbf{T},$

where for a given l and m, ${\mathsf{\Lambda }}^{\left(E\right)}$ is the diagonal matrix in the energy representation

(37)$\left|\right|{\mathsf{\Lambda }}^{\left(n\right)}\left|\right|=\left|\right|⟨n^{\prime }lm|\hslash \omega (\widehat{n}+3/2)|nlm⟩||=\hslash \omega (n+3/2){\delta }_{n^{\prime }n},$

while in the position representation

(38)$\left|\right|{{\mathsf{\Lambda }}^{\prime }}^{\left(r\right)}\left|\right|=-\frac{m{\omega }^2}{2}{r}_i^2{\delta }_{ij}+V\left(\sqrt{r_i^2}\right){\delta }_{ij}.$

In Equation (36), $\kappa$ stands for a control parameter in the interval $[0,1]$, whose variation allows the potential to be transformed from a harmonic oscillator ($\kappa =0$) to the arbitrary potential $V\left(\sqrt{r^2}\right)$ ($\kappa =1$).

The methodology of this approach is simple. We diagonalize both matrix representations (22) and (23) to obtain the coefficients (30) and (31). Then we proceed to diagonalize in the energy representation the Hamiltonian (36) associated with the potential we are interested in for a given angular momentum and projection.

2.3. $SU\left(4\right)$-UGA

This approach starts by establishing the $U\left(4\right)$ dynamical group for 3D systems by adding a scalar boson $s^{†}\left(s\right)$ to the space of the same bosonic operators $a_{\mu }^{†}\left({\tilde{a}}_{\mu }\right)$ previously introduced carrying angular momentum $l=1$ in such a way that the bilinear products satisfy the commutation relations associated with the generators of the $U\left(4\right)$ group [19]:

(39)$\begin{array}{cc}\hfill {\widehat{n}}_a& =\sqrt{3}{[a\times \tilde{a}]}_0^{\left(0\right)};{\widehat{n}}_s=s^{†}s,\hfill \end{array}$

(40)$\begin{array}{cc}\hfill {\widehat{L}}_{\mu }& =\sqrt{2}{[a\times \tilde{a}]}_{\mu }^{\left(1\right)};{\widehat{D}}_{\mu }=a_{\mu }^{†}s-s^{†}{\tilde{a}}_{\mu },\hfill \end{array}$

(41)$\begin{array}{cc}\hfill {\widehat{Q}}_{\mu }& ={[a\times \tilde{a}]}_{\mu }^{\left(2\right)};{\widehat{R}}_{\mu }=i(a_{\mu }^{†}s+s^{†}{\tilde{a}}_{\mu }),\hfill \end{array}$

with the constraint that the total number of bosons $\widehat{N}={\widehat{n}}_a+s^{†}s$ is constant. Hence, the associated quantum number N establishes the dimension of the dynamical space. From the dynamical group $U\left(4\right)$, three subgroup chains can be obtained, preserving the angular momentum [18,19]:

(42)$\begin{array}{cc}\hfill U\left(4\right)& \supset U\left(3\right)\supset SO\left(3\right)\supset SO\left(2\right),\hfill \end{array}$

(43)$\begin{array}{cc}\hfill U\left(4\right)& \supset SO\left(4\right)\supset SO\left(3\right)\supset SO\left(2\right),\hfill \end{array}$

(44)$\begin{array}{cc}\hfill U\left(4\right)& \supset \overline{SO}\left(4\right)\supset SO\left(3\right)\supset SO\left(2\right).\hfill \end{array}$

Each chain provides a basis with a definitive meaning, as will become clear shortly. A particular chain, (42), leads to the basis $|\left[N\right]n_{a}LM⟩$ defined by the operator equations

(45)$\begin{array}{cc}\hfill \widehat{N}|\left[N\right]n_{a}LM⟩=& N|\left[N\right]n_{a}LM⟩,\hfill \end{array}$

(46)$\begin{array}{cc}\hfill {\widehat{n}}_a|\left[N\right]n_{a}LM⟩=& n_a|\left[N\right]n_{a}LM⟩,\hfill \end{array}$

(47)$\begin{array}{cc}\hfill {\widehat{L}}^2|\left[N\right]n_{a}LM⟩=& L(L+1)|\left[N\right]n_{a}LM⟩,\hfill \end{array}$

(48)$\begin{array}{cc}\hfill {\widehat{L}}_z|\left[N\right]n_{a}LM⟩=& M|\left[N\right]n_{a}LM⟩.\hfill \end{array}$

We now turn our attention to the connection with configuration space. Using the approach recently proposed [20] with a mapping to the 3D harmonic oscillator basis

(49)$|\left[N\right]n_{a}LM⟩\to |nlm⟩.$

It is possible to show that the $U\left(4\right)$ algebraic realization ${\widehat{\mathcal{Q}}}_{\mu }$ and ${\widehat{\mathcal{P}}}_{-\mu }$ for the coordinates and momenta, respectively, take the form [16]

(50)$\begin{array}{cc}\hfill {\widehat{\mathcal{Q}}}_{\mu }=& \frac{d}{\sqrt{2N}}{\widehat{D}}_{\mu },\hfill \end{array}$

(51)$\begin{array}{cc}\hfill {\widehat{\mathcal{P}}}_{-\mu }=& -{(-)}^{1-\mu }\frac{\hslash }{d}\frac{1}{\sqrt{2N}}{\widehat{R}}_{\mu }.\hfill \end{array}$

The mapping

(52)$q_{\mu }\to {\widehat{\mathcal{Q}}}_{\mu }{p}_{\mu }\to {\widehat{\mathcal{P}}}_{\mu }$

allows the algebraic Hamiltonian (1) to be obtained in the $U\left(4\right)$ space:

(53)${\widehat{H}}_{alg}^{U\left(4\right)}=\frac{1}{2m}{\widehat{\mathcal{P}}}^2+\frac{m{\omega }^2}{2}{\widehat{\mathcal{Q}}}^2.$

Notice that the operators ${\widehat{\mathcal{P}}}^2$ and ${\widehat{\mathcal{Q}}}^2$ are not related in a straightforward way with the physical momentum $p^2$ and coordinate $q^2$, but their matrix elements coincide in the N large limit:

(54)$\underset{N\to \infty }{lim}\left|\right|⟨N{n}_a^{\prime }LM|{\widehat{\mathcal{P}}}^2|N{n}_{a}LM⟩\left|\right|=\left|\right|⟨n^{\prime }lm|p^2|nlm⟩||,$

with similar equation for the coordinate.

Taking into account (50) and (51), the Hamiltonian (53) is recast in the form:

(55)${\widehat{H}}_{alg}^{U\left(4\right)}=\frac{\hslash \omega }{4}\left[\frac{{\widehat{R}}^2+{\widehat{D}}^2}{N}\right].$

The second order Casimir operators of $U\left(4\right)$ and $U\left(3\right)$ are given by [18]

(56)$\begin{array}{cc}\hfill {\widehat{C}}_2\left[U\left(4\right)\right]=& \frac{1}{3}{\widehat{n}}_a^2+\frac{1}{2}({\widehat{L}}^2+{\widehat{D}}^2+{\widehat{R}}^2)+{\widehat{Q}}^2+{(\widehat{N}-{\widehat{n}}_a)}^2,\hfill \end{array}$

(57)$\begin{array}{cc}\hfill {\widehat{C}}_2\left[U\left(3\right)\right]=& \frac{1}{3}{\widehat{n}}_a^2+\frac{1}{2}{\widehat{L}}^2+{\widehat{Q}}^2,\hfill \end{array}$

and in the symmetrical space of identical bosons the following replacements are valid [18]

(58)$\begin{array}{cc}\hfill C_2\left[U\left(4\right)\right]\to & \widehat{N}(\widehat{N}+3)\hfill \end{array}$

(59)$\begin{array}{cc}\hfill C_2\left[U\left(3\right)\right]\to & {\widehat{n}}_a({\widehat{n}}_a+2),\hfill \end{array}$

which allows (55) to be transformed into

(60)${\widehat{H}}_{alg}^{U\left(4\right)}=\hslash \omega \left[\left(1-\frac{1}{N}\right){\widehat{n}}_a+\frac{3}{2}-\frac{{\widehat{n}}_a^2}{N}\right].$

Here we have taken into account that the operator $\widehat{N}$ is diagonal in any basis and consequently it can be substituted by its eigenvalue N. This result leads to the identification of the $U\left(3\right)$ chain (42) as the energy representation with eigenkets

(61)${\widehat{n}}_a|\left[N\right]n_{a}LM⟩=n_a|\left[N\right]n_{a}LM⟩,$

and branching rules $n_a=0,1,\dots ,N$ and $L=n_a,n_a-2,\dots ,1\mathrm{or}0$.

Let us now come back to chains (43,44). Since the Casimir operators characterizing the $SO\left(4\right)$ and $\overline{SO}\left(4\right)$ groups are ${\widehat{W}}^2={\widehat{D}}^2+{\widehat{L}}^2$ and ${\widehat{\overline{W}}}^2={\widehat{R}}^2+{\widehat{L}}^2$ respectively, Equation (50) and (51) implies that chain (43) provides the coordinate representation while chain (44) gives the momentum representation. The corresponding eigenvectors satisfy

(62)$\begin{array}{cc}\hfill {\widehat{W}}^2|\left[N\right]\zeta LM⟩=& \zeta (\zeta +2)\left|\right[N]\zeta LM⟩,\hfill \end{array}$

(63)$\begin{array}{cc}\hfill {\widehat{\overline{W}}}^2|\left[N\right]\overline{\zeta }LM⟩=& \overline{\zeta }(\overline{\zeta }+2)|\left[N\right]\overline{\zeta }LM⟩,\hfill \end{array}$

with branching rules $\zeta ,\overline{\zeta }=N,N-2,\dots ,1\mathrm{or}0$, and $L=0,1,\dots ,\zeta \left(\overline{\zeta }\right)$. The importance of these transformations is twofold. On one hand, we have the mapping (49) and consequently it is possible to project any state to the position representation as an expansion of harmonic oscillator functions. On the other hand, an arbitrary potential depending on $\sqrt{{\widehat{\mathcal{Q}}}^2}$ takes a diagonal form in the basis (62), providing a direct way to obtain its matrix representation in the energy representation through a similarity transformation.

Let us now consider a general situation of a 3D Hamiltonian of the form (32). The corresponding algebraic Hamiltonian is obtained by applying the anharmonization procedure (52):

(64)${\widehat{H}}_{alg}^{U\left(4\right)}=\widehat{H}{|}_{q\to \widehat{\mathcal{Q}},p\to \widehat{\mathcal{P}}}=\frac{1}{2m}{\widehat{\mathcal{P}}}^2+\widehat{V}\left(\sqrt{{\widehat{\mathcal{Q}}}^2}\right).$

A practical convenient form for this Hamiltonian is obtained by adding and subtracting a quadratic term in such a way that the harmonic oscillator term (1) is identified:

(65)${\widehat{H}}_{alg}^{U\left(4\right)}=\hslash \omega \left[\left(1-\frac{1}{N}\right){\widehat{n}}_a+\frac{3}{2}-\frac{{\widehat{n}}_a^2}{N}\right]+\kappa {\widehat{V}}^{\prime }\left(\sqrt{{\widehat{\mathcal{Q}}}^2}\right).$

where

(66)$V^{\prime }\left({\widehat{\mathcal{Q}}}^2\right)=-\frac{m{\omega }^2}{2}{\widehat{\mathcal{Q}}}^2+\widehat{V}\left(\sqrt{{\widehat{\mathcal{Q}}}^2}\right).$

Taking into account the scalar character of ${\widehat{\mathcal{Q}}}^2$, we have its following matrix elements in the coordinate representation:

(67)$⟨\left[N\right]{\zeta }^{\prime }LM|{\widehat{\mathcal{Q}}}^2|\left[N\right]\zeta LM⟩=\frac{d^2}{2}\frac{\left[\zeta \right(\zeta +2)-L(L+1\left)\right]}{N}{\delta }_{{\zeta }^{\prime },\zeta }.$

Hence, the matrix elements of the Hamiltonian (for a given L and M) in the energy representation take the general form

(68)${\mathbf{H}}^{\left(E\right)}={\mathsf{\Lambda }}^{\left(E\right)}+\kappa {\mathbf{T}}^{†}{\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}\mathbf{T},$

where ${\mathsf{\Lambda }}^{\left(E\right)}$ is the diagonal contribution of the deformed harmonic oscillator

(69)$\left|\right|{\mathsf{\Lambda }}^{\left(E\right)}\left|\right|=\hslash \omega \left[\left(1-\frac{1}{N}\right)n_a+\frac{3}{2}-\frac{n_a^2}{N}\right]{\delta }_{n_a^{\prime },n_a},$

while ${\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}$ is the diagonal matrix of the potential $V^{\prime }\left(\sqrt{{\widehat{\mathcal{Q}}}^2}\right)$ in the position representation:

(70)$\left|\right|{\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}\left|\right|=\hslash \omega \left[-\frac{1}{2}\frac{\xi {(\zeta ,L)}^2}{2N}+\frac{1}{\hslash \omega }V\left(d\frac{\xi (\zeta ,L)}{\sqrt{2N}}\right)\right]{\delta }_{\zeta ,{\zeta }^{\prime }},$

with $\xi (\zeta ,L)=\sqrt{\zeta (\zeta +2)-L(L+1)}$. The $\mathbf{T}$ matrix stands for the transformation brackets: $\mathbf{T}=\left|\right|⟨\left[N\right]\zeta LM|\left[N\right]n_{a}LM⟩||$, which are obtained by diagonalizing the matrix representation of the Casimir operators ${\widehat{W}}^2$ given by

$\begin{array}{cc}\hfill & ⟨\left[N\right]n_{a}LM\left|{\widehat{W}}^2\right|\left[N\right];n_{a}LM⟩=\hfill \\ \hfill & N(2n_a+3)-2n_a(n_a+1)+L(L+1),\hfill \\ \hfill & ⟨\left[N\right]n_a+2,LM\left|{\widehat{W}}^2\right|\left[N\right];n_{a}LM⟩=\hfill \\ \hfill & -\sqrt{(N-n_a-1)(N-n_a)(n_a-L+2)(n_a+L+3)},\hfill \end{array}$

or in analytic form [18,21]. This matrix representation of ${\widehat{W}}^2$ is tridiagonal and its diagonalization to provide $\mathbf{T}$ is calculated once and for all for a given N, while the diagonal matrix ${\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}$ is modified in accordance with the specific potential. The contribution (69) involves accidental degeneracy due to the quadratic contribution ${\widehat{n}}_a^2$. This degeneracy manifests in the full Hamiltonian (68) by the appearance of degeneracy for the bound states in the N large limit. This degeneracy is removed by substituting the diagonal term (69) by [16]

(71)$\left[\left(1-\frac{1}{N}\right){\widehat{n}}_a+\frac{3}{2}-\delta \frac{{\widehat{n}}_a^2}{N}\right],$

with $\delta =1\mathrm{for}{n}_a<N/2$ and $\delta =0\mathrm{for}{n}_a\ge N/2.$

The expression (68) contains the main features of our approach. The physical information characterizing the system is given by the diagonal matrix ${\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}$ easily generated and given by (70), while the calculation of the transformation brackets allows any system to be analyzed in simple form using the standard tools of matrix diagonalization. It is important to notice that this approach is intended to be considered in the N large limit [16]. This is in contrast to the previous analysis of the algebraic models where $N\to \infty$ represents the classical limit [22].

3. Coulombic Potential: Hydrogen Atom

Because of its importance in dealing with molecular systems the Coulombic potential deserves special attention. For this reason we shall consider the non-relativistic Hydrogen atom together with the effect of an external electric field to show the advantages of the proposed approaches. We start considering the Hamiltonian of the Hydrogen atom.

Following M. Moshinsky, [17], in energy units of the first Bohr orbit $E_B=\frac{1}{2}m{e}^4/{\hslash }^2$, the Hamiltonian takes the form

(72)$\frac{\widehat{H}}{E_B}=\frac{1}{2m{E}_B}{p}^2-\sqrt{\frac{2{\hslash }^2}{m{E}_B}}\frac{1}{r}.$

Here, the mass associated with the harmonic oscillator is identified with the reduced mass of the Hydrogen system. Adding and subtracting the harmonic oscillator potential we obtain

(73)$\frac{\widehat{H}}{E_B}={ϵ}^2{\widehat{\overline{H}}}^{\mathrm{HO}}-\left[\frac{{ϵ}^2}{2}{\overline{r}}^2+ϵ\sqrt{2}\frac{1}{\overline{r}}\right],$

where $\overline{r}$ is assumed to be in units of $d=\sqrt{\hslash /m\omega }$, the harmonic oscillator Hamiltonian ${\widehat{\overline{H}}}^{\mathrm{HO}}$ in units of $\hslash \omega$ and $ϵ$ a dimensionless parameter given by

(74)$ϵ=\sqrt{\frac{\hslash \omega }{E_B}},$

with $\omega$ a variational parameter to be determined.

Let us now establish, on one hand, the application of the HO-DVR approach to this problem. The matrix representation of the Hamiltonian in units of $\hslash \omega$ takes the form

(75)$\mathbf{H}={ϵ}^2{\mathsf{\Lambda }}^{\left(n\right)}-\kappa {\mathbf{T}}^{†}{{\mathsf{\Lambda }}^{\prime }}^{\left(r\right)}\mathbf{T},$

where the diagonal matrix ${\mathsf{\Lambda }}^{\left(n\right)}$ in the energy representation is given by

(76)$\left|\right|{\mathsf{\Lambda }}^{\left(n\right)}\left|\right|=(n+3/2){\delta }_{n^{\prime }n},$

while for the diagonal matrix ${{\mathsf{\Lambda }}^{\prime }}^{\left(r\right)}$ in the position representation we have

(77)$\left|\right|{{\mathsf{\Lambda }}^{\prime }}^{\left(r\right)}\left|\right|=\left[\frac{{ϵ}^2}{2}{\overline{r}}_i^2+ϵ\sqrt{2}\frac{1}{\sqrt{{\overline{r}}_i^2}}\right]{\delta }_{ij}.$

On the other hand, in the framework of the $U\left(4\right)$-DVR, we introduce the realization (50) and (51) into (73) to obtain

(78)$\frac{\widehat{H}}{E_B}={ϵ}^2\left[\left(1-\frac{1}{N}\right){\widehat{n}}_a+\frac{3}{2}-\frac{{\widehat{n}}_a^2}{N}\right]-\left[\frac{1}{2}{ϵ}^2\left(\frac{{\widehat{D}}^2}{2N}\right)+ϵ\frac{2\sqrt{N}}{\sqrt{{\widehat{D}}^2}}\right].$

The corresponding matrix representation takes the form

(79)${\mathbf{H}}^{\left(E\right)}={\mathsf{\Lambda }}^{\left(E\right)}-\kappa {\mathbf{T}}^{†}{\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}\mathbf{T},$

where ${\mathsf{\Lambda }}^{\left(E\right)}$ is the diagonal contribution of the deformed harmonic oscillator

(80)$\left|\right|{\mathsf{\Lambda }}^{\left(E\right)}\left|\right|={ϵ}^2\left[\left(1-\frac{1}{N}\right)n_a+\frac{3}{2}-\frac{n_a^2}{N}\right]{\delta }_{n_a^{\prime },n_a},$

while ${\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}$ is the diagonal matrix containing the potential with the form

(81)$\left|\right|{\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}\left|\right|=\left[\frac{1}{2}{ϵ}^2\left(\frac{1}{2N}\right)\xi {(\zeta ,L)}^2+\frac{ϵ2\sqrt{N}}{\xi (\zeta ,L)}\right]{\delta }_{\zeta ,{\zeta }^{\prime }},$

where the preservation of the angular momentum has been considered implicitly. It is importance to notice that because of the branching rule $\zeta =N,N-2,\dots ,1\mathrm{or}0$, the number of bosons N must be taken to be odd in order to avoid a singularity (obtained when $L=\zeta =0$).

The diagonalization of the Hamiltonian matrices (75) and (79) lead to the correlation diagram displayed in Figure 1 as a function of $\kappa$. The parameter (74) turns out to be $ϵ=1$. Although the correlations diagram is the same for both methods, the corresponding basis dimension is different in order to reach convergence. As noticed, the levels corresponding to the principal quantum number in the Hydrogen atom converge to the well known $n^2$ degeneracy, although to simplify the figure we have only included angular momenta up to 2. The degeneracy $2L+1=2l+1$ corresponding to the projection $(M,m)$ should be assumed but not shown because of the $(M,m)$-independence of the Hamiltonian. The states above the zero energy represent a discretization of the continuum. This kind of the description of the continuum is also present in the Morse and Pöeschl–Teller potentials [23,24,25].

In order to explicitly show the energy convergence, the comparison between the exact energy levels $E_n=-1/n^2$ and the calculated ones are displayed in Table 1 for the HO-DVR method and in Table 2 for the $U\left(4\right)$-UGA. The convergence is similar, but while in the HO-DVR it is reached for $N=2001$, the $U\left(4\right)$-UGA needed a number of bosons $N=4001$. This difference may be explained by the accidental degeneracy involved in the $U\left(4\right)$-UGA but also because the bases are different. In both cases, the space dimension is large enough to obtain a reasonable convergence to the exact values. In general, the convergence can be measured through the fidelity F or overlap of eigenstates involving different parameters N [26]. Along this venue we have calculated the energies for number of bosons $N\pm 50$, registering the energy difference as the error displayed in Table 3 for the HO-DVR method and in Table 4 for the $U\left(4\right)$-UGA. It is interesting that the uniformity of the errors in the latter case. On the other hand, we know the analytical solutions and consequently the wave functions can be directly compared.

In our approach the $\alpha$-th eigenstate in the HO-DVR method takes the form

(82)$|{\psi }_{\alpha ,lm}^N⟩=\sum _{n=0}^N⟨nlm|{\psi }_{\alpha ,lm}^N⟩|nlm⟩,$

while for the $U\left(4\right)$-UGA

(83)$|{\psi }_{\alpha ,LM}^N⟩=\sum _{n_a=0}^N⟨\left[N\right]n_{a}LM|{\psi }_{\alpha ,LM}^N⟩|\left[N\right]n_{a}LM⟩.$

From these expressions we obtain the corresponding position representation for the radial contribution of the wave function. Integrating over the angular coordinates and applying the coordinate projection we obtain $⟨r|{\psi }_{\alpha L}^N⟩$ or $⟨r|{\psi }_{\alpha l}^N⟩$ with

(84)$⟨r|nl⟩=⟨r|\left[N\right]n_{a}L⟩=A_{nl}{r}^l{e}^{r^2/2}{L}_{(n-l)/2}^{l+1/2}\left(r^2\right),$

where $A_{nl}=\sqrt{(2^{l+2}(n-l)!!)/(\sqrt{\pi }(n+l+1)!!)}$. The projections $⟨r|{\psi }_{\alpha L}^N⟩$ or $⟨r|{\psi }_{\alpha l}^N⟩$ can be compared to the exact analytical functions $R_{nL}\left(r\right)$. Here the variable r is given in the distance unit d. Since $ϵ=1$ we have that $\hslash \omega =E_B$, a result that determines the distance units in the MKS units system. In Figure 2 the position projection of the eigenstates (82) and (83) together with the exact solutions for the $1s,2s$ and $3s$ orbitals are displayed, while in Figure 3 the corresponding $2p,3p$ and $3d$ orbitals are shown. Beyond the intrinsic convergence ($N=2001$ for the HO-DVR method and $N=4001$ for the U(4)-UGA) provided by the harmonic oscillator basis, we notice that the wave functions are very well described, a feature that proves the validity of both approaches for the free Hydrogen atom.

4. Stark Effect

In this section, we study the effect of an static external electric field on the spectrum of the Hydrogen atom. The electric field interacts through the dipole producing a splitting of the spectral lines called Stark effect, since this phenomenon was experimentally shown by Stark in 1913 [27]. There is a vast amount of literature on the Stark effect in the Hydrogen atom. This problem has been treated in configuration space using semi-classical methods [28,29,30], as well as variational approaches [31,32], but also algebraically in an elegant way. Since it is impossible to mention all the contributions to this problem we confine our discussion to some relevant works to put in context our algebraic method. We do not pretend that our method represents an improvement to the previous descriptions. Instead we are interested in showing that it is possible to tackle this problem along the venue worked out by M. Moshinsky, where the harmonic oscillator was the basis to solve problems from atoms to quarks [17].

Algebraic methods applied to the Hydrogen atom may be considered to have originated from the identification of the $SO\left(4\right)$ group as the symmetry group of the non-relativistic Hydrogen atom for the bound states and $SO(3,1)$ for the continuum [33,34,35,36,37,38,39]. Based upon the $SO\left(4\right)$ symmetry group the Stark effect for both static and time-dependent electric fields has been described from different points of view, including semiclassical methods taking into account the crucial result that the Runge–Lenz vector is proportional to the coordinate and the energy [40,41,42,43]. In any case, these treatments are confined to the subspace characterized for a given principal quantum number n since the generators of the group commute with the Hamiltonian, a fact that constraints the descriptions to perturbative schemes. In addition the natural basis to deal with this problem, called the Stark basis $|n_1{n}_{2}m⟩$, is based on the parabolic coordinates where the Schrödinger equation is separable. In this context an important role is played by the transformation brackets $⟨n_1{n}_{2}m|nlm⟩$ connecting to the states $|nlm⟩$ with good angular momentum [44].

A full algebraic treatment for the Stark effect needs the introduction of the $SO(4,2)$ dynamical group [38,45,46,47,48] or alternative groups [49]. The $SO(4,2)$ group can be identified by the merging of the $SO(2,1)$ and $SO\left(4\right)$ groups [50]. The $SO(2,1)$ group, or alternatively, the $SU(1,1)$ group deals with the radial contributions that shift the principal quantum numbers [50,51,52]. Even though the $SO(4,2)$ group allows the dipole matrix elements to be calculated, the use of the dynamical group tends to perturbation treatments [53,54]. In contrast, our algebraic approach is not based on the use of the dynamical group, instead, our approach falls in the framework of a discrete variable representation approach, whose salient feature consists in its simplicity as we next present.

In the previous section, we have shown that both algebraic DVR approaches are able to describe the free Hydrogen atom, although with a different convergence degree. The result is that the HO-DVR method needs half of the number of the basis functions, albeit with somewhat less uniformity of errors. In this section, the interaction with a static electric field is introduced whose better description is tempting to be given by the HO-DVR approach. The result, however, is that the $SU\left(4\right)$-UGA is more suitable to describe the Stark Effect given the dimensions taken into account for the free Hydrogen atom, as we show next.

The Hydrogen atom under the action of the electric field takes the form

(85)${\widehat{H}}^S=\widehat{H}+\widehat{V}\left(E\right),$

where $\widehat{H}$ represents the unperturbed Hydrogen atom (72) and $\widehat{V}\left(E\right)$ is the potential energy associated with the interaction of the electron with the external electric field $\left|\mathbf{E}\right|$ along the z direction

(86)$\widehat{V}\left(E\right)=-e\left|\mathbf{E}\right|z,$

and in the convenient units previously used

(87)$\frac{\widehat{V}\left(E\right)}{E_B}=-\frac{\overline{z}}{ϵ}2\sqrt{2}\left|\overline{\mathbf{E}}\right|,$

with $\overline{z}$ in units of d and $\overline{\mathbf{E}}$ in atomic units defined by the electric field of the nucleus given by $e/a_0^2\approx 5\times {10}^9$ V/cm.

Let us now consider the realization of the potential $\widehat{V}\left(E\right)$ in accordance with the two algebraic DVR methods. According to the relations (6) and (7) the z component is given by

(88)$\overline{z}=\frac{1}{2}(a_0^{†}+a_0).$

Hence, in the HO-DVR method we just calculate the matrix elements of the operator $a_0^{†}$ in the basis (9). We obtain

(89)$\begin{array}{cc}\hfill ⟨n+1,l+1,m\left|a_0^{†}\right|nlm⟩& =\sqrt{(n+l+3)}\sqrt{\frac{(l+1)}{(2l+3)}}\sqrt{\frac{(l+m+1)(l-m+1)}{(2l+1)(l+1)}},\hfill \end{array}$

(90)$\begin{array}{cc}\hfill ⟨n+1,l-1,m\left|a_0^{†}\right|nlm⟩& =-\sqrt{(n-l+2)}\sqrt{\frac{l}{(2l-1)}}\sqrt{\frac{(l+m)(l-m)}{l(2l+1)}}.\hfill \end{array}$

On the other hand, in the U(4)-UGA method, taking into account (50) and (51) with the correspondence $\overline{z}\to \widehat{\mathcal{Q}}/d$, we have

(91)$\overline{z}=\frac{1}{\sqrt{2N}}{\widehat{D}}_0.$

Using the Wigner–Eckart theorem, its matrix elements take the form

(92)$\begin{array}{cc}\hfill ⟨\left[N\right]{n^{\prime }}_a{L}^{\prime }M\left|\widehat{D}\right|\left[N\right]n_{a}LM⟩=& ⟨LM;10|L^{\prime }M⟩⟨\left[N\right]{n^{\prime }}_a{L}^{\prime }\left|\right|\widehat{D}\left|\right|\left[N\right]n_{a}L⟩,\hfill \end{array}$

with the reduced matrix elements

(93)$\begin{array}{cc}\hfill ⟨\left[N\right]n_a-1L-1⟩\left|\right|\widehat{D}\left|\right|\left[N\right]n_{a}L⟩& =-\sqrt{\frac{(N-n_a+1)(n_a+L+1)L}{2L-1}},\hfill \end{array}$

(94)$\begin{array}{cc}\hfill ⟨\left[N\right]n_a-1L+1⟩\left|\right|\widehat{D}\left|\right|\left[N\right]n_{a}L⟩& =-\sqrt{\frac{(N-n_a+1)(n_a-L)(L+1)}{2L+3}},\hfill \end{array}$

(95)$\begin{array}{cc}\hfill ⟨\left[N\right]n_a+1L-1⟩\left|\right|\widehat{D}\left|\right|\left[N\right]n_{a}L⟩& =\sqrt{\frac{(N-n_a)(n_a-L+2)L}{2L-1}},\hfill \end{array}$

(96)$\begin{array}{cc}\hfill ⟨\left[N\right]n_a+1L+1⟩\left|\right|\widehat{D}\left|\right|\left[N\right]n_{a}L⟩& =\sqrt{\frac{(N-n_a)(n_a+L+3)(L+1)}{2L+3}}.\hfill \end{array}$

Taking into account these technical ingredients, we carried out the diagonalization of the Hamiltonian (85) in the framework of both methods.

In Figure 4 the energy levels associated with the subspace $n=2$ as a function of the electrical field are depicted. In order to obtain these results, a cutoff of the dimension space determined by the feasibility of the calculations was necessary as well as a constraint in the involved angular momenta. The black solid circles correspond to the spectrum obtained using the $U\left(4\right)$-UGA, where the number of boson was taken to be $N=3501$ up to L = 16, while the blue triangles correspond to HO-DVR method taking N = 3001 up to l = 11. The continuous lines correspond to the energy-splitting obtained by the variational approach described in Ref. [32], where an accurate variational method is presented for the n = 2−5 shells for different $m^{\prime }s$ including a display of the wave functions. Since, for each approach, the convergence of the energy levels for the null field is different the continuous lines were appropriately located to follow the splitting trend. There are two remarkable facts to be highlighted. First, the splitting is correctly described by the $SU\left(4\right)$-UGA, while the HO-DVR method tends to overestimate the splitting as the field intensity increases. This behavior is explained by the different basis dimensions. As a second point, we notice that for small-field intensity the splitting is rather flattened, showing an evident deviation from the linearity predicted by perturbation theory. In order to clearly see this effect a zoom of Figure 4 is depicted in Figure 5 for low field intensities in the interval $[0,{10}^{-4}]$. Both approaches present this anomalous behavior, which we believe corresponds to a lack of convergence. To make clear that this is in fact the case, a similar plot is depicted in Figure 6 for the HO-DVR approach for different boson numbers. The black triangles correspond to $N=3001$, while the red ones to $N=1501$. From these plots, it becomes clear that as the basis dimension increases the splittings approach to the correct results.

It is worth mentioning the complexity of these calculations due to the spherical symmetry breaking. Due to our computational limitations, we are forced to fix a maximum value for the angular momentum. Without this approximation, we would be dealing with a ${\left(\mathrm{3,067,752}\right)}^2$ dimension matrices, which are impossible to calculate given our computational capacity. However, the restriction of the angular momentum components has physical justification in the sense that not all the functions with a given projection are expected to be mixed. Indeed for the $n=2$ and $n=3$ levels the mixing of states with different angular momenta are not expected to be so large (may be up to $L,l=5$), but in this work we have included the maximum angular momenta that allows our calculations to be carried out.

In order to see whether the previous behavior for the p-space is reproduced for other subspaces, in Figure 7 we display the effect of the electric field provided by both algebraic DVR methods over the subspace with $n=3$ for $M=m=0$. The same trend is obtained given the comparison with the exact result given in Ref. [32]. While the $SU\left(4\right)$-UGA reproduces the correct splitting, the HO-DVR method overestimates it.

Finally it is convenient to recall the usual treatment of the Stark effect using perturbation theory. For the first excited state ($n=2$) the fourfold degeneracy must be born in mind. In this context, since the angular momentum projection is preserved and z is odd, only the states $|{\mathsf{\Psi }}_{200}⟩$ and $|{\mathsf{\Psi }}_{210}⟩$ interact through (86). The rest of the matrix element vanish leading to the eigenvalues for $M=m=0$ [55,56]:

(97)$E_{\pm }=E_2\pm 6\left|\overline{\mathbf{E}}\right|,$

where the energy is intended to be in units of $E_B$ and the electric field in atomic units. The magnitude of the splitting of the levels is proportional to the electric field strength, giving rise to the linear Stark effect. The straight lines provided by the equation (97) coincide with the exact values displayed in Figure 4. Deviation from linearity manifests around field intensities of 0.004 au.

5. Conclusions

In this work, we have applied two algebraic discrete variable representation approaches to describe the non-relativistic Hydrogen atom together with the Sark effect. The selection of the Coulombic potential to test our approaches is twofold. On one hand, its study is compulsory to consider more complex situations like molecular systems. On the other hand it is reasonable to evaluate the feasibility of these approaches to take into account additional interactions breaking the symmetry, in particular for the $SU\left(4\right)$-UGA where the condition (71) is imposed to remove the intrinsic degeneracy associated with the deformed harmonic oscillator (60).

The two approaches we have proposed may be, in principle, applied to general potentials to obtain the solutions of the time-independent Schrödinger equation. In both approaches, the matrix representation of the Hamiltonian given by either (36) or (79) is given in terms of two diagonal matrices together with the transformation brackets. The calculation of the transformation brackets is independent of the potential and hence they are calculated only once for a given space dimension. The diagonal matrices ${\mathsf{\Lambda }}^{\left(E\right)}$ corresponding to the simple or deformed harmonic oscillator is trivially generated and fixed for any system, while ${\mathsf{\Lambda }}^{\left(r\right)}$ and ${\mathsf{\Lambda }}^{\left(\mathcal{Q}\right)}$ are diagonal in the coordinate representation which are also easily constructed. Hence, these approaches represent powerful simple alternatives to study systems where the analytical solutions of the Schrödinger equations are either not known or difficult to obtain by integration methods. Even considering moderate values of space dimension, the approach provides solutions reflecting the main salient features of the system, a situation very useful when the solution cannot be obtained by conventional methods. Hence, this approach may be useful in catastrophe theory where a small change in the original stable potential changes the topology of the problem [57], as well as in the study of quantum phase transitions [58], but also provides new ways to obtain Franck–Condon factors between different potentials since in our approach a common basis of harmonic oscillators are used [59].

We have shown the feasibility of both approaches, the HO-DVR method and U(4)-UGA to obtain the solutions of the 3D Coulombic potential. Although both algebraic methods provide the solutions, the HO-DVR approach proves to show a faster convergence, a fact that makes it the best candidate to study this system. A quite different situation is manifested when the dipole interaction is taken into account to describe the Stark effect, which is particularly difficult due to the spherical symmetry breaking. The full dimension of the matrices to be diagonalized is too large to be consider in realistic terms and consequently, we were forced to cutoff the included angular momenta. With this approximation, the calculations become feasible with reasonable convergence. Fixing the total number of bosons obtained to describe the free Hydrogen atom, the $SU\left(4\right)$-UGA turns out to be the suitable method to describe the energy splittings, while the HO-DVR method overestimates the Stark effect. Although both approaches are based on the harmonic oscillator basis, the localized functions associated with the discrete variable representation basis are different, a fact that makes the difference between these approaches. We have thus obtained the unexpected remarkable result that a preference of a model based on the convergence in the zeroth-order Hamiltonian does not imply the same preference when a symmetry breaking interaction is included.

The two presented algebraic approaches provide alternatives to obtain the solutions of 3D systems for any potential where the harmonic oscillator functions represent a good basis. Even though for the Coulombic potential the harmonic oscillator is not by far the best basis, our methods provide reasonable results in a simple way. This opens the possibility of studying systems with different algebraic DVR methods where the exact solutions are not known, providing the possibility of extracting properties by purely algebraic means using the standard tools of matrix diagonalization.

Author Contributions

M.R.-A., M.B.-M.: investigation, software, methodology, reviewing, discussion. R.L.: conceptualization, supervision, investigation, original draft preparation, discussion. J.M.A., J.G.-C.: conceptualization, investigation, writing, discussion. E.O.: investigation, discussion. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

This work is partially supported by DGAPA-UNAM, Mexico, under project IN-212020, by the Consejería de Economía, Conocimiento, Empresas y Universidad de la Junta de Andalucía (Spain) under Group FQM-160, by the Spanish Ministerio de Ciencia e Innovaciíon, ref. FIS2017-88410-P, RTI2018-098117-B-C21 and PID2019-104002GB-C22, and by the European Commission, ref. H2020-INFRAIA-2014-2015 (ENSAR2). First author is grateful to DGAPA-UNAM for the postdoctoral scholarship (Facultad de Química). Second author is also grateful for the scholarship (Posgrado en Ciencia e Ingeniería de Materiales) provided by CONACyT, México.

Conflicts of Interest

The authors declare no conflict of interest.

Figures and Tables

Figure 1. Correlation diagram from the harmonic oscillator (κ=0) to the Coulomb potential (κ=1) corresponding to the Hamiltonians (36) and (68), with (77) and (81), respectively. At the right of the figures the principal quantum number for the electron in the H atom, n, is indicated. To simplify only the levels of the Hydrogen system with L,l=0,1,2 are included. In the bottom panel a zoom is displayed in order to show the convergence of the levels with different angular momenta. The top level in each group corresponds to angular momentum L,l=0. The parameters were taken to be ϵ=1 and N=4001 for the U(4)-unitary group approach (UGA), while N=2001 for the HO-DVR method.

Figure 2. Comparison between exact and calculated radial wave functions Rn0(r);n=1,2,3 for the Hydrogen atom using (75) and (79) with parameters N=2001 and N=4001 respectively, with ϵ=1. The dash lines correspond to the exact wave functions.

Figure 3. Comparison between exact and calculated radial wave functions R21(r),R31(r),R32(r) for the Hydrogen atom using (75) and (79) with parameters N=2001 and N=4001 respectively, with ϵ=1. The dash lines correspond to the exact wave functions.

Figure 4. Effect of the electric field effect over the subspace n=2 with m=M=0 in the Hydrogen atom. The black solid circles correspond to spectrum obtained with the U(4)-UGA taking N=3501 up to L = 16, while the blue triangles correspond to HO-DVR method taking N = 3001 up to l = 11. Since the calculations were carried out independently, the plotted points do not coincide. The continuous line correspond to the variational approach described in Ref. [32]. |E¯| is the electric field in atomic units.

Figure 5. Electric field effect over the n=2 levels of Hydrogen atom for M=m=0. The black solid circles correspond to the results using the U(4)-UGA method taking N = 3501 and L = 16. The blue triangles correspond to HO-DVR method taking N = 3001 and l = 11.

Figure 6. Zoom of the electric field effect over the subspace n=2 with m=M=0 in the Hydrogen atom provided by the HO-DVR method for N=1501 (red triangles) and N=3001 (black triangles).

Figure 7. Electric field effect over the subspace characterized by the principal quantum number n=3 with M=m=0 in the Hydrogen atom. The solid black dots correspond to the U(4)-UGA method taking N = 3501 up to L = 16, while the blue triangles correspond to the HO-DVR method taking N = 3001 up to l = 11. The continuous lines correspond to the exact results provided by Ref. [32].

Table 1

Exact bond energies compared with the energies provided by the HO-DVR method, the basis dimension was taken to be $N=2001$.

	Basis Dimension N
	2001
${\mathit{E}}_{\mathit{n}}$	$\mathit{L}=\mathbf{0}$	$\mathit{L}=\mathbf{1}$	$\mathit{L}=\mathbf{2}$
−1.	−0.9969
−0.25	−0.2496	−0.2500
−0.1111	−0.1111	−0.1111	−0.1111
−0.0625	−0.0625	−0.0625	−0.0625
−0.04	−0.0400	−0.0400	−0.0400
−0.0278	−0.0275	−0.0276	−0.0276
−0.0204	−0.0173	−0.0175	−0.0180

Table 2

Exact bond energies compared with the energies provided by the diagonalization of (68) taking a total number of bosons $N=4001$ with the $U\left(4\right)$-UGA method.

	Total Number of Bosons N
	4001
${\mathit{E}}_{\mathit{n}}$	$\mathit{L}=\mathbf{0}$	$\mathit{L}=\mathbf{1}$	$\mathit{L}=\mathbf{2}$
−1.	−1.00330
−0.25	−0.25088	−0.25039
−0.1111	−0.11237	−0.11222	−0.11220
−0.0625	−0.06463	−0.06457	−0.06455
−0.0400	−0.04344	−0.04340	−0.04336
−0.0278	−0.0330	−0.0330	−0.0330
−0.0204	−0.02541	−0.02555	−0.02586

Table 3

Errors of the energies provided by the HO-DVR method. Errors were calculated considering the difference of energies $\left|E\right(N=2051)-E(N=1951\left)\right|$.

Error for $\mathit{L}=0$	Error for $\mathit{L}=1$	Error for $\mathit{L}=2$
0.0001
0.00002	9 × 10${}^{-9}$
5 × 10${}^{-6}$	3 × 10${}^{-9}$	8 × 10${}^{-13}$
2 × 10${}^{-6}$	1 × 10${}^{-9}$	1 × 10${}^{-11}$
2 × 10${}^{-6}$	5 × 10${}^{-7}$	3× 10${}^{-7}$
0.00008	0.00007	0.00005
0.0006	0.0005	0.0005

Table 4

Errors of the energies provided by the $U\left(4\right)$-UGA method. Errors were calculated considering the difference of energies $\left|E\right(N=4051)-E(N=3951\left)\right|$.

Error for $\mathit{L}=0$	Error for $\mathit{L}=1$	Error for $\mathit{L}=2$
0.00007
0.00002	0.00001
0.00003	0.00003	0.00003
0.00005	0.00005	0.00005
0.00009	0.00009	0.00009
0.0001	0.0001	0.0001
0.00008	0.00006	0.00003