Spectral Analysis of Proton Eigenfunctions in Crystalline Environments

The Schr\"odinger equation and Bloch theorem are applied to examine a system of protons confined within a periodic potential, accounting for deviations from ideal harmonic behavior due to real-world conditions like truncated and non-quadratic potentials, in both one-dimensional and three-dimensional scenarios. Numerical computation of the energy spectrum of bound eigenfunctions in both cases reveals intriguing structures, including bound states with degeneracy matching the site number $N_w$, reminiscent of a finite harmonic oscillator spectrum. In contrast to electronic energy bands, the proton system displays a greater number of possible bound states due to the significant mass of protons. Extending previous research, this study rigorously determines the constraints on energy gap and oscillation amplitude of the previously identified coherent states. The deviations in energy level spacing identified in the computed spectrum, leading to minor splitting of electromagnetic modes, are analyzed and found not to hinder the onset of coherence. Finally, a more precise value of the energy gap is determined for the proton coherent states, ensuring their stability against thermal decoherence up to the melting temperature of the hosting metal.


I. INTRODUCTION
Over the last decades, extensive theoretical research and the development of new experimental techniques have made possible the manipulation of coherent and collective phenomena in quantum optics at the nano-scale level, such as superradiance and coherent population trapping [1].These phenomena often involve quantum phase transitions [2] and an ultrastrong coupling regime of light/matter interaction [3].Their theoretical description must take into account dissipation, as is the case in general for open quantum systems [4][5][6][7][8], which are becoming increasingly relevant for applications to nano-electronics [9].Superradiance and the Dicke effect occur when an ensemble of molecules or quantum oscillators confined within a sub-wavelength region collaboratively emit and absorb coherent radiation.Coherence among these emitters is mediated by the electric radiation field.In some cases, the field-matter coupling can be significantly enhanced by the presence of surface a Email address: luca.gamberale@unimib.itb Email address: giovanni.modanese@unibz.itplasmons, resulting in what is known as the plasmonic Dicke effect [10][11][12][13][14][15][16][17][18].
In this work we continue our investigation of an idealized system comprising N s charges oscillating within a lattice structure, exhibiting a superradiant transition above a specific threshold density.The dynamics of this system involve bulk plasmons and enhanced fieldmatter coupling, topics that we studied analytically in [19], where we demonstrated in the large N s limit the existence of an energy gap for a specific set of wavefunctions displaying coherence in both the matter and field sectors.
The oscillating charges that we have considered in our model are protons, with reference to the concrete physical case of hydrogen loading in metals, where protons are bounded to tetrahedral or octahedral sites.The large mass of protons, compared to electrons, simplifies the treatment because it allows to disregard the fermionic character of their wavefunctions, at least in certain conditions.It also turns out that the large mass of protons plays a crucial role in determining the energy spectrum and the spatial extension of their bounded wavefunctions.Here lies the major new contribution of the present work towards an improvement of the model in Ref. [19].In [19] we have described the bound states of protons as those of ideal harmonic oscillators, introducing an ad-hoc physical cutoff on their oscillation amplitudes.Here we provide a solid theoretical justification to the existence of the limitation in the oscillation amplitude and compute numerically the exact states of the protons in a crystal structure with periodic but finite potential wells.As could be intuitively anticipated, the proton wavefunctions are much more localized than those of electrons.Furthermore, the count of bound excited states is precisely determined by the crystal structure.For instance, in the case of octahedral voids with parameters as detailed in Table I, there are precisely 12 bound states.Only these states actively contribute to the superradiance phenomenon.This elucidates and justifies the previously mentioned cutoff on oscillation amplitude.
In [20], we recall the Hamiltonian of our model and the canonical transformations of the photon field operators used to diagonalize the photonic term.That work succinctly summarizes key findings from [19] in a self-consistent manner, and presents numerical simulations showing the occurrence of a quantum phase transition for small values of N s , also in the presence of dissipation.
The final Hamiltonian in the dipole approximation takes the form: Here, C i and C † i represent linear combinations of photon annihilation and creation operators, obtained by projecting along the three spatial directions êi and summing over all possible momentum directions k.Using these operators, the vector potential simplifies to A = The term H osc in (1) denotes the sum of the Hamiltonians of the proton oscillators: , where the a n are standard destruction operators for the harmonic oscillators.For the definition of the frequencies ω, ω p , ω ′ , please refer below.We utilize units in which ℏ = c = 1.
One of the key features of the model, as extensively elaborated in reference [19], is its consistent definition of the natural frequency denoted as ω for all oscillators within the material.This uniformity is established by employing a periodic electrostatic potential known as the "jellium crystal".During the diagonalization process of the total Hamiltonian, the frequency ω is combined with the plasma frequency ω p = e 2 N s /(mV ) associated with the oscillating charges.This combination yields a modified or "dressed" frequency denoted as ω ′ = ω 2 + ω 2 p .It is noteworthy that the photon momentum, represented as k, remains unchanged throughout this process.Consequently, within this model, there exist no states in which electromagnetic energy generated within the material can propagate into the vacuum, thereby forming a naturally occurring resonating cavity.
Another crucial aspect of the emerging coherent electromagnetic modes, which oscillate with a fixed phase relation to the matter oscillators, is the orientation of their wave vectors in all possible spatial directions.This enhances the coupling between matter and electromagnetic fields, enabling the formation of the energy gap.
While the overall framework seems largely satisfactory in its primary aspects, a notable question remains regarding size of the oscillation amplitude of the coherent state, a parameter directly influencing the energy gap's magnitude.In this investigation, we approach this issue employing a rigorous physical methodology, thereby dispelling any lingering remnants of heuristic speculation.This is the main result of the present work.
In the subsequent sections, we will focus on the bound states of protons bound to crystal sites.In Sect.II we find the spectrum through a numerical solution of the Schrödinger equation and taking into account the Bloch theorem.In Sect.IV we consider a perturbation term which accounts for deviations in energy level spacing, resulting in minor e.m. mode splitting.
Finally, we rigorously calculate the maximum amplitude of the coherent oscillation, which had been manually set using a heuristic argument in [19], and summarize our findings.

II. SINGLE-PARTICLE BOUND STATES IN A FINITE QUADRATIC PERIODIC POTENTIAL
In reference [19], we tackled the issue of one-particle bound states within an ideal periodic infinite potential well.In idealized scenarios, such wells provide a simplified framework for understanding quantum systems.However, when considering real-world conditions, oscillators, or particles, deviate from ideality due to truncated and non-quadratic potentials.This departure from ideal behavior carries significant consequences for the structure of bound energy levels.In the subsequent analysis, we will delve into the profound implications of these deviations, examining how they impact the characteristics of bound energy states and, more importantly, its consequences in the properties of their coherent states.
Consider an ensemble of protons with mass m immersed in a periodic quasi-harmonic truncated potential in 3D of the form where ω is the free vibration frequency, ⃗ x cijk are the 3D positions of the lattice sites and where q max is a parameter that controls the width and depth of the potential well (the argument q of V(q) is a mute variable).The coordinates ⃗ x cijk are selected to represent various periodic crystal structures.We opt to utilize the face-centered cubic (FCC) symmetry, characterized by where i, j, k ), and a denotes the crystal spacing (a = 3.52 Å for Nickel).The Heaviside (step) function θ ensures the potential remains non-positive, as depicted in Fig. 1 along the x-direction with y = z = 0.The Schrödinger equation that determines the eigenstates of the single-particle problem is whose solution can be afforded using the Bloch theorem.By setting and working in the tight binding approximation [21] we can set where ϕ ⃗ n (⃗ x) are the eigenfunctions of the Schrödinger equation for the isolated potential The validity of this approximation stems from the observation that the particles under examination have a mass much larger than that of the electrons, resulting in their wavefunctions displaying significant localization around the centers of the potential wells.This localization, in turn, guarantees that the level of overlap between the wave-functions associated with neighboring sites can be considered practically negligible to a considerable degree at least up to a certain level of excitation of the bound states.
The energy spectrum of the bound eigenfunctions described by Equation ( 5) can, in fact, be derived from the spectrum of the equation independent of the momentum vector ⃗ K, that can take any value.
where the 1D potential is given by ( , where N w is the even number of wells considered in the calculation) Through a rotational symmetry argument, it can be demonstrated that the 1D problem is analogous to the 3D problem and shares the same energy spectrum.Specifically, it can be shown that Eq. ( 8) corresponds to the radial component of the Schrödinger equation in three dimensions (3D) for a single potential well (Eq.( 7)), focusing on the reduced radial solutions characterized by zero angular momentum ("s-wave" solutions).Consequently, the numerical results we aim to present remain valid even within the framework of the complete 3D scenario.
By setting the discretization points x j = εj, −(M − 1)/2 ≤ j ≤ (M − 1)/2 (M arbitrary odd integer) where the discretization step is ε = d Nc where N c is the number of discretization points per cell, Eq. ( 8) is rewritten in discretized form as where T is the second derivative operator in discretized form and is defined by and ψ is a vector whose components are the values of ψ at the points x i : The final parameter requiring specification is the selection of the value for q max .In the context of addressing the challenge of hydrogen loading in metals within the octahedral phase, prior literature indicates that the available radius within the octahedral voids is approximately 0.41 d/2.Thus, our decision is to set q max = 0.41 The eigenvalue problem of Eq. ( 10) is now completely set and can be solved numerically.
The numerical solution with the parameters indicated in Tab.I is shown in Fig. 2 and reveals a very interesting structure.The bound states with negative energy have a degeneracy equal to the site number N w as is expected for the tight binding case and their spectrum is similar to that of the harmonic oscillator.The free states (positive energy) follow a typical Brillouin dispersion relation, since the degeneracy is removed by the fact that the tight binding approximation breaks down, the wave-functions not being confined in their respective lattice sites.The eigenfunctions of the bound states are very similar to the eigenfunctions of a harmonic oscillator centered in a specific site of the lattice (see Fig. 3).FIG.3: Eigenfunctions of the bound states computed by direct diagonalization of Eq. ( 8) with the parameters of Table I.
One might inquire about the difference between the well-established outcomes regarding electronic energy bands, which exhibit a significantly distinct structure compared to the solutions obtained here for protons.The answer lies in the fact that, given that protons are approximately 1800 times more massive than electrons, the lower-energy wave functions are considerably more localized than those of electrons.This localization arises from the spatial scaling factor being 1 √ mω , consequently leading to a significantly greater number of possible bound states.In fact, the same calculation conducted for electrons reveals that the maximum number of bound states is at most 2.
A closer inspection reveals that the negative energy levels are not exactly equally spaced, resulting in several frequencies of the electromagnetic field being involved.Fig. 4 shows the dependence of the photon energies for the various couplings between adjacent energy levels.
Such an issue will be dealt with in Section IV.FIG.4: Linear fit of the photon energies associated to the dipolar transitions between adjacent bound energy levels (the last two points to the right are excluded from the fit).

III. THREE-DIMENSIONAL ANALYSIS
We now extend our numerical analysis to the three-dimensional (3D) case, albeit without delving into the complexities associated with the resulting band structure for free eigenstates.Our primary objective is to demonstrate that even in the 3D scenario, bound states exhibit uniform energy spacing, a crucial factor for coherence establishment.
Our goal is to corroborate the idea that the bound states of an isolated oscillator bound in a single potential well of the form (9) has an energy spectrum practically identical to that of a 3D harmonic oscillator with the same mass and oscillation frequency confined by a full harmonic potential (with no truncation) with minimum potential energy − 1 2 mω 2 R w .Equation (7) describes a single particle interacting with a spherically symmetric potential well centered at the origin of the coordinates.The potential V, always non-positive, delineates a finite set of bound states with negative energy and a continuous set of improper eigenfunctions with positive energy.However, in the spectral sector of our interest, the eigenfunctions of the harmonic oscillator decay exponentially at large distances, becoming negligible even at a distance from the origin of less than d 2 .We thus develop the spectral decomposition of Eq. ( 7) based on the eigenfunctions of the 3D harmonic oscillator with an unbounded potential centered at the origin V h (r) = 1 2 mω 2 (r 2 − R 2 w ) with R w = (mω) −1/2 q max .Consequently, the eigenfunctions and eigenvalues of the harmonic oscillator with negative energies closely resemble those of Equation ( 7) with a high degree of accuracy.
The eigenfunctions of the 3D harmonic oscillator are given by and the corresponding eigenvalues are where n = 2n r + l is the principal quantum number, n r is the radial quantum number, and l and m are the angular quantum numbers.Here, Y lm (x) represents the spherical harmonics, and R h nrl (r) denotes the radial eigenfunctions of the 3D harmonic oscillator.These radial eigenfunctions are given by where ρ = √ mωr, L nr are the generalized Laguerre polynomials, and are normalization factors.
To validate our approach, we numerically solved the discretized eigenvalue equation for a radial potential with truncated harmonic behavior at radius R w .Utilizing a MATLAB algorithm, we discretized the eigenvalue equation using the finite element method on a uniform lattice.The discretization method follows the same line of the 1D case, aside from the different choice of the coordinate range and the adoption of the reduced wave-function ϕ(r) → ϕ(r) r .The radial potential considered comprises a truncated harmonic potential plus the centrifugal term and is given by where Here, N ∞ defines the dimensionality of the eigenvalue problem and is chosen to be much larger than the number of bound eigenvalues of interest.The parameter δ is a small number compared to ε and prevents numerical overflow when evaluating the potential at r = 0. We solved the system for a range of principal quantum numbers n and angular quantum numbers l, with a maximum value of n max = 12 chosen as a cutoff parameter.
We observed that the system's eigenvalues are influenced solely by the principal quantum number n with no dependence on the angular quantum number l, thus reproducing the expected analytical result.Remarkably, this property remains valid even for eigenvalues very close to zero, where the influence of the non-harmonicity of the potential may be significant.
In Table II are shown the computed eigenvalues for the bound states together with the corresponding theoretical values for the 3D harmonic oscillator and the results are found to be in excellent agreement with theoretical predictions.The confirmation of independence from angular momentum, as displayed by the eigenstates of the full harmonic oscillator, persists even for eigenvalues very close to the free spectrum.

IV. QUASI-HARMONIC CASE
In [19] we have analyzed the conditions for the onset of electromagnetic coherence for a system of N charged oscillators immersed in a periodic harmonic potential and we have heuristically inferred the conditions for the non-divergence of the energy gap and oscillation amplitude.More specifically, the divergence of the cited quantities depends directly on the unlimited number of the bounded energy levels of the harmonic oscillator.
We now demonstrate that when considering the more realistic scenario in which the periodic potential is confined as outlined in Section II, both the energy gap and the oscillation amplitude are naturally constrained without any additional assumption.
Another complicating factor in our analysis arises from the non-uniform spacing of energy levels within the oscillators, leading to the coupling of several electromagnetic field modes with the oscillators, as opposed to just a single mode.
To tackle this intricacy, we pivot towards the simpler 1D scenario, underlining that analogous conclusions apply to the 3D counterpart.We commence by scrutinizing the idealized scenario where all oscillator levels, labeled as E i = (i + 1 2 )ω 0 , exhibit uniform spacing, alongside introducing a minor perturbation term of 1 2 δ ω i(i − 1) to accommodate deviations.This adjustment yields an energy disparity between successive levels of E i+1 − E i = ω 0 + δ ω i (ω 0 denotes the unperturbed frequency).
Commencing with δ ω = 0, we initially observe that the coherent state comprises a solitary frequency ω 0 .Subsequently, we incrementally increase δ ω by an infinitesimal amount, resulting in an exceedingly minor splitting of the various electromagnetic modes.For the sake of continuity, we anticipate that the state characterized by the minimum energy still consists of a solitary electromagnetic mode.
Upon conducting a thorough examination of the phase transition mechanism, we ascertain that, in order to generate a negative contribution of the interaction term throughout a single oscillation period of the oscillators, it is imperative that the phase differential between the modes with frequencies ω i originating from the i-th energy level and the resulting frequency ω ′ associated with the coherent mode remains confined to a value less than ϵπ, where ϵ is the coupling constant of Eq. ( 26) in [19].
The above condition can be implemented by requiring where ω ′ is taken as the averaged oscillation frequencies of the em field and and where δ ω is the coefficient of variation of ω i with the level index (see Fig. 4).
Eq. ( 16) bears resemblance to the critical condition observed in the Kuramoto model [22] for an ensemble of oscillators with marginally varying natural frequencies.The theoretical constructs delineated in [19] could potentially draw parallels with the Kuramoto model, thereby furnishing a theoretical underpinning from the vantage point of quantum field theory.This correlation will be investigated in an upcoming inquiry.
Setting the average frequency ω ′ as the mean position with respect to the involved frequencies so that we can solve Eqs. ( 15), ( 16), (17) for the variable i, obtaining the maximum value as i max = 2ϵω 0 |δω|(1+ϵ) .Utilizing the parameters in Table I and from Appendix A in [19], we determine that the number of bound states is N b = i max = 12 (see Fig. 2) for the 1D case, and N b = i max = 11 for the 3D case, one less than in the 1D case due to the zero point energy contribution (see Table II).Consequently, we can infer that all these bound states contribute to the establishment of coherence, resulting in an averaged frequency of ω ′ = 0.4 eV.
Finally, we can compute the maximum amplitude of the coherent oscillation that is the parameter that ultimately controls the size of the energy gap (see Eqs. (30) and (43b) of [19] for the up-to-1st and up-to-2nd order calculation).
The number of photons within the coherent state corresponds to the available dipolar transitions in the system, which in our context amounts to N b − 1.Additionally, for a Glauber state, the expected number of photons is |α| 2 , thereby setting |α| 2 = N b − 1.
Utilizing Eq. ( 34) from [19], and given ω ′ = 0.41 eV for protons (Eq.(A24) of [19]), we derive for the 3D case an estimate for the oscillation amplitude: which contrasts with the heuristic estimate of 0.4 in Ref. [19], established on intuitive grounds alone.Furthermore, the energy gap value shifts from the estimated 1 eV to the more precise 0.37 eV.

V. CONCLUSION
In this investigation, we computed the maximum amplitude of coherent oscillation of protons in a metal hydride, constrained by the periodic electrostatic potential generated by the electron distribution within the metallic crystal.This calculation, derived from fundamental principles, plays a fundamental role in the theoretical framework describing the quantum coherence of a plasma comprising protons confined to crystal lattice sites by a quasi-harmonic potential.Moreover, this parameter dictates the ultimate energy gap attained by the coherent state.
Previously, the determination of this maximum amplitude relied on heuristic considerations, guided by the understanding that protons, tethered to tetrahedral or octahedral sites within each crystal unit cell, cannot transition to adjacent cells during their oscillatory motion, particularly within the realm of rapid dynamics.However, it is acknowledged that over extended time scales, proton diffusion within the crystal lattice occurs (indeed, this process underpins the initial hydrogen loading of the system).
It is noteworthy that the energy gap obtained is significantly greater than the average thermal energy characteristic of the crystal's operational temperature range (an average thermal energy of 0.37 eV corresponds to a temperature above 4000 K).This ensures the stability of the coherent state against thermal decoherence.In essence, these coherent structures, once established within the metallic crystal, constitute a crucial element of proton dynamics and necessitate consideration in situations where collective phenomena hold paramount significance.
In order to find the maximum number of quasi-harmonic energy levels involved in the coherence process, that are directly related to the coherent oscillation amplitude, we have solved the Schrödinger equation in the crystal (a periodic solution, according to the Bloch theorem), after reducing the full 3D problem to 1D and defining a potential which is still harmonic but with a finite depth, and therefore has a spectrum with both bound and free states.
The calculated one-particle energy spectrum, employing realistic physical parameters, exhibits distinctive bands that markedly differ from electronic bands.This disparity arises from the considerable mass discrepancy between protons and electrons, coupled with the significantly more localized nature of proton wavefunctions in comparison to electrons.
The nonuniform spacing of energy levels, stemming from the non-harmonicity of the potential, introduces the complexity of multiple electromagnetic modes rather than a singular one.To address this challenge, we introduced a perturbation term to accommodate deviations in energy level spacing, leading to slight level splitting.Importantly, our analysis demonstrates that this splitting does not impede the onset of coherence, and all bound modes actively contribute to the process.
To summarize, the limited number of energy levels involved in the coherent transition immediately resolves the problem related to the unlimited energy gap foreseen in Ref. [19] and the related unboundedness of the oscillation amplitude.In Ref. [19] we imposed "by hand" an upper bound the coherent oscillation through a heuristic physical argument.In the present paper we theoretically substantiate such an ansatz through an explicit calculation and we find a result essentially consistent with the previous result.
The energy gap identified in our analysis guarantees that the coherent proton states spontaneously formed within metal hydrides possess sufficient stability against thermal decoherence, even up to the melting temperatures of the metal itself.

FIG. 2 :
FIG. 2: One particle energy spectrum calculated using the parameters of Table I.Note the harmonic-oscillator-like energy spacing of the bound states and their degeneracy, equal to the number of lattice sites N w .The free eigenvalues have a Brillouin-like structure.Different bans correspond to different colors.The finite number of wells causes intermediate eigenvalues to appear between adjacent bands (see sub-figure).
FIG.5: Eigenvalues of the radial equation for all quantum numbers of the bound states.The confirmation of independence from angular momentum, as displayed by the eigenstates of the full harmonic oscillator, persists even for eigenvalues very close to the free spectrum.

TABLE I :
Chosen parameters for the numerical calculation