text:
# Dynamical clustering in oscillator ensembles with time-dependent interactions

## Abstract

We consider an ensemble of coupled oscillators whose individual states, in addition to the phase, are characterized by an internal variable with autonomous evolution. The time scale of this evolution is different for each oscillator, so that the ensemble is inhomogeneous with respect to the internal variable. Interactions between oscillators depend on this variable and thus vary with time. We show that as the inhomogeneity of time scales in the internal evolution grows, the system undergoes a critical transition between ordered and incoherent states. This transition is mediated by a regime of dynamical clustering, where the ensemble recurrently splits into groups formed by varying subpopulations.

Synchronization phenomena in interacting dynamical systems have attracted a great deal of attention in the last two decades bocca1 . The spontaneous appearance of coherent evolution as a result of interactions is ubiquitous in a vast class of natural and artificial systems, and can be successfully reproduced by relatively simple models. Most studies of this kind of phenomena have focused on patterns of highly correlated evolution, such as full and phase synchronization. Though these forms of synchronization are essential to certain artificial systems—for instance, secure-communication and chaos-control devices—they are expected to play a more restricted role in such objects as biological populations, where the complexity of their collective function requires a delicate balance between behavioral coherence and diversity Mikh . A clear-cut example is the human brain, where highly coherent functional patterns are only realized under pathological states—notably, during epileptic seizures. Neural tissues, as well as many natural systems ranging from intracellular molecular complexes to social populations, are normally found in segregated configurations, where the system splits into groups of elements with distinct functions. In these clustered states, highly coherent evolution occurs inside each group, while the collective dynamics of different groups are much less correlated. This is typically a dynamical phenomenon, since—as a consequence of the evolution of the internal state of each element—clusters may temporarily disperse to be reconstituted later, and individual elements or small groups can migrate between clusters.

Clustering has been observed in ensembles of coupled dynamical systems with specific individual dynamics kcl ; zcl , or subject to sufficiently complex interaction functions golomb ; okuda . Disintegration of clusters and migration of elements is typically observed under the action of noise kcl ; noise ; kori . In this Letter, we show that such kind of dynamical clustering can occur in an ensemble of coupled periodic oscillators, in the absence of external fluctuations, when the interaction between oscillators depends on individual internal state variables with autonomous evolution. The time scale of this internal evolution differs between oscillators. A critical transition between ordered and incoherent evolution, mediated by dynamical clustering, takes place as the inhomogeneity of such time scales grows.

Consider $`N`$ coupled phase oscillators, each of them characterized by a phase variable $`\varphi _i(t)`$, obeying by the equations

$$\dot{\varphi }_i=N^1\underset{j=1}{\overset{N}{}}J_{ij}(t)\mathrm{sin}(\varphi _j\varphi _i)$$
(1)
($`i=1,\mathrm{},N`$). When the interaction weights $`J_{ij}(t)`$ are constant, $`J_{ij}(t)J`$ for all $`i`$ and $`j`$, these equations reduce to Kuramoto’s model for identical phase oscillators kura , which exhibits full synchronization for any positive value of $`J`$. Time-independent non-identical interaction weights have also been considered, disclosing typical features of disordered systems, such as glassy-like behavior, frustration, and algebraic relaxation towards equilibrium d622 ; d1073 . Here, we analyze the case where the interaction weights depend on the internal state of each oscillator, which is specified by a variable $`\theta _i(t)`$ with autonomous evolution. Specifically, we consider that $`\theta _i`$ is itself a phase variable evolving according to $`\dot{\theta }_i=\mathrm{\Omega }_i`$, where the frequency $`\mathrm{\Omega }_i`$ is chosen at random for each element, from a fixed distribution $`g(\mathrm{\Omega })`$. The interaction weight is given by

$$J_{ij}(t)=J\mathrm{cos}[\theta _j(t)\theta _i(t)],$$
(2)
where $`J>0`$ is a constant factor that can be fixed to $`J=1`$ by rescaling time. With this choice for $`J_{ij}`$, the interaction between the phase $`\varphi _i`$ and $`\varphi _j`$ is attractive when $`\mathrm{cos}(\theta _j\theta _i)`$ is positive, and repulsive otherwise. The sign of $`J_{ij}`$ changes in a time scale of order $`|\mathrm{\Omega }_j\mathrm{\Omega }_i|^1`$, giving rise—at the level of the whole ensemble—to a complex time-dependent interaction pattern. This evolving connection network makes possible the appearance of several dynamical regimes, including clustering.

We start our study of the above model by performing extensive numerical realizations in ensembles ranging in size from $`N=10^2`$ to $`10^6`$. The internal frequencies $`\mathrm{\Omega }_i`$ are drawn at random from a Gaussian distribution, $`g(\mathrm{\Omega })=\mathrm{exp}[(\mathrm{\Omega }\mathrm{\Omega }_0)^2/2\sigma ^2]/\sqrt{2\pi \sigma ^2}`$. Since the dynamics is invariant under a constant shift in these frequencies, we fix $`\mathrm{\Omega }_0=0`$. It turns out that the collective behavior of the ensemble is sensible to the frequency dispersion $`\sigma `$. For small $`\sigma `$ interactions evolve slowly, over long time scales. Consequently, the phase $`\varphi _i`$ can adiabatically follow the evolution of the interaction weights. For long times, the system is found in a state where the phases $`\varphi _i`$ are homogeneously distributed over the interval $`[0,2\pi )`$. A strong correlation is observed in such state between $`\varphi _i`$ and the internal variable $`\theta _i`$, namely, $`\varphi _i\pm \theta _i+\varphi _0`$, where $`\varphi _0`$ is a constant. The sign factor of $`\theta _i`$ and the value of $`\varphi _0`$ are the same for all oscillators, and depend on the initial condition; moreover, $`\varphi _0`$ can slowly change as time elapses. For large $`\sigma `$, on the other hand, many of the interaction weights exhibit large changes over very short times. Each oscillator is thus subject to rapidly fluctuating forces, which again results in a state where the phases $`\varphi _i`$ are homogeneously distributed over $`[0,2\pi )`$. Now, however, there is no correlation between $`\varphi _i`$ and $`\theta _i`$.

To characterize the transition between the regimes of small and large frequency dispersion we introduce, for each oscillator, the two-dimensional complex vector $`𝐦_i=(\mathrm{exp}[\mathrm{i}(\varphi _i\theta _i)],\mathrm{exp}[\mathrm{i}(\varphi _i+\theta _i)])`$, and define the average

$$𝐦=N^1\underset{j=1}{\overset{N}{}}𝐦_j(\mu _+\mathrm{exp}(\mathrm{i}\psi _+),\mu _{}\mathrm{exp}(\mathrm{i}\psi _{})).$$
(3)
The time average $`\mu `$ of $`\sqrt{\mu _+^2+\mu _{}^2}`$ is a suitable order parameter for the transition. If the system spends most of the time in the state where $`\varphi _i\pm \theta _i+\varphi _0`$ for all $`i`$, we have $`\mu 1`$, while $`\mu N^{1/2}`$ in the incoherent state. Figure 1 shows numerical results for $`\mu `$ as a function of $`\sigma `$ for various system sizes. The dependence with $`N`$ suggest the presence of a critical phenomenon for $`\sigma 0.3`$. Within a few assumptions supported by numerical evidence, this critical behavior can be analytically studied in the limit $`N\mathrm{}`$, as follows. First of all, we note that for $`\sigma 0`$ and at sufficiently long times, the internal variables $`\theta _i`$ are homogeneously distributed over $`[0,2\pi )`$. If, in the limit $`\sigma 0`$ (but $`\sigma 0`$), the temporal variation of $`\theta _i`$ is disregarded, $`\varphi _i=\pm \theta _i+\varphi _0`$ is a stationary state of our system, as direct substitution in Eq. (1) shows. This conclusion holds for any distribution of frequencies highly concentrated around $`\mathrm{\Omega }=0`$. For larger values of $`\sigma `$, we write $`\varphi _i=\pm (\theta _i\delta _i)+\varphi _0`$, where $`\delta _i(t)`$ measures the deviation of each oscillator with respect to the small-dispersion state. This deviation satisfies

$$\dot{\delta }_i=\mathrm{\Omega }_i\frac{\mu }{2}\mathrm{sin}\delta _i,$$
(4)
where we have written the order parameter $`\mu `$ in terms of the distribution of the deviations $`\delta _i`$ over the ensemble, $`p(\delta )`$, as

$$\mu =_{\pi /2}^{\pi /2}p(\delta )\mathrm{cos}\delta d\delta .$$
(5)
Here, we have assumed that $`p(\delta )`$ is an even function; as shown below, this amounts to postulate that $`g(\mathrm{\Omega })`$ is itself even.

Equation (4) has been extensively discussed in connection with phase-oscillator ensembles saka . For $`|\mathrm{\Omega }_i|<\mu /2`$ it has a stable fixed point at one of the solutions of $`\mathrm{sin}\delta _i=2\mathrm{\Omega }_i/\mu `$. For $`|\mathrm{\Omega }_i|>\mu /2`$, on the other hand, there are no fixed points and $`|\delta _i|`$ grows indefinitely with time. The ensemble can therefore be thought of as consisting of two subpopulations. Those oscillators for which $`|\mathrm{\Omega }_i|<\mu /2`$ attain, at asymptotically large times, a stationary deviation $`\delta _i`$, which depends on the value of $`\mathrm{\Omega }_i`$. For the oscillators with $`|\mathrm{\Omega }_i|>\mu /2`$, in contrast, the phase $`\varphi _i`$ does not reach a stationary value with respect to $`\theta _i`$. We refer to these two subpopulations as subensembles I and II, respectively. The distribution $`p(\delta )`$ can be split into contributions from each subensemble, $`p=p_\mathrm{I}+p_{\mathrm{II}}`$. The first contribution is found from the distribution of frequencies, taking into account the identity $`p_\mathrm{I}(\delta )d\delta =g(\mathrm{\Omega })d\mathrm{\Omega }`$. As for $`p_{\mathrm{II}}`$, it can be shown that its contribution to the integral of Eq. (5) vanishes. This equation becomes, then,

$$\mu =_{\mu /2}^{\mu /2}g(\mathrm{\Omega })\sqrt{14\mathrm{\Omega }^2/\mu ^2}𝑑\mathrm{\Omega },$$
(6)
making it possible to find $`\mu `$ self-consistently for any (even) distribution of frequencies. Under quite general conditions on $`g(\mathrm{\Omega })`$, this equation describes a second-order transition between disordered ($`\mu =0`$) and ordered ($`\mu 0`$) states. For $`\mu 0`$, two states are possible, either with $`\mu _+0`$ or with $`\mu _{}0`$, which demonstrates the occurrence of symmetry breaking at the transition. The curve in Fig. 1 shows the solution to Eq. (6) for the Gaussian distribution of frequencies considered in our numerical realizations. It vanishes as $`\mu |\sigma \sigma _c|^{1/2}`$ at $`\sigma _c=\sqrt{\pi /32}0.313`$.

Note that Eq. (1) can be written, in terms of the average quantity introduced in Eq. (3), as

$$\dot{\varphi }_i=\frac{\mu _+}{2}\mathrm{sin}(\psi _++\theta _i\varphi _i)+\frac{\mu _{}}{2}\mathrm{sin}(\psi _{}\theta _i\varphi _i).$$
(7)
This form of the equation of motion for $`\varphi _i`$ emphasizes the mean-field nature of interactions. The fact that, as $`N\mathrm{}`$, $`\mu `$ vanishes for $`\sigma >\sigma _c`$, implies that the two terms in the right-hand side of Eq. (7) vanish as well. For finite system sizes, they are of order $`N^{1/2}`$. Hence, beyond the critical dispersion $`\sigma _c`$, the evolution of phases can be thought of as driven by fluctuating forces of order $`N^{1/2}`$. In the thermodynamical limit, for $`\sigma >\sigma _c`$, the dynamics is frozen. Below the transition, evolution becomes progressively slower as the critical point is approached.

Numerical simulations show that for intermediate values of $`\sigma `$, as the predominance of the ordered small-dispersion state decreases and the critical transition is approached, a complex regime of clustering develops. Those oscillators in subensemble II with frequencies just above their lowest value, $`|\mathrm{\Omega }_i|\mu /2`$, become divided into four clusters in $`\varphi `$-space. Though the structure of each cluster is not sharply localized, it is still possible to define its phase by averaging $`\varphi _i`$ over the oscillators in that cluster. Two of the clusters are formed by oscillators with $`\mathrm{\Omega }_i>0`$, and exhibit anti-phase synchronization, so that their phases differ by $`\pi `$ at all times. They move monotonically around $`\varphi `$-space. The other two move in the opposite direction. They are formed by oscillators with $`\mathrm{\Omega }_i<0`$, and also show anti-phase synchronization (see inset of Fig. 1). Due to the opposite motion of the two pairs of clusters, they recurrently cross each other. At these crossings, their populations are temporarily superimposed in phase and, at the same time, their motion decelerates and the clusters become more compact. As a result, two well-localized and relatively long-lived big clusters, their phases differing by $`\pi `$, are recurrently built up out of the ensemble (see Fig. 2c).

The occurrence of two-cluster states is quantitatively disclosed by the complex quantity

$$z=N^1\underset{j=1}{\overset{N}{}}\mathrm{exp}(2\mathrm{i}\varphi _j).$$
(8)
If the whole ensemble splits into two anti-phase well-localized clusters, we have $`|z(t)|1`$, while for incoherent or higher-order clustered states we have $`|z(t)|N^{1/2}`$. Figure 2a shows the evolution of $`|z(t)|`$ in the clustering regime, displaying its irregular oscillations between relatively small and large values. For the same realization, Fig. 2b displays the phase difference $`\mathrm{\Delta }\varphi `$ between two oscillators in subensemble II, both with positive values of $`\mathrm{\Omega }_i`$, as a function of time. When $`|z(t)|`$ is large, $`\mathrm{cos}\mathrm{\Delta }\varphi 1`$ identifies a two-cluster state where the two oscillators belong to the same cluster, while $`\mathrm{cos}\mathrm{\Delta }\varphi 1`$ corresponds to the situation where the oscillators are in different clusters. The intermittent transitions of $`\mathrm{cos}\mathrm{\Delta }\varphi `$ between large positive and negative values reveal the dynamical nature of clustering, where any two oscillators may alternatively belong to the same or to different clusters.

The time average $`\zeta `$ of the quantity $`|z(t)|`$ plays the role of an order parameter detecting the two-cluster state. Figure 3 displays the dependence of $`\zeta `$ with the frequency dispersion $`\sigma `$, for ensembles of several sizes. A peak near $`\sigma _c`$ becomes clearly defined as $`N`$ grows. While for $`\sigma >\sigma _c`$ the value of $`\zeta `$ apparently vanishes for large system sizes, a well-defined profile persists below the critical point. It grows from zero starting at $`\sigma 0.2`$ and attains its maximum, $`\zeta 0.25`$, at $`\sigma \sigma _c`$.

The formation of the four clusters whose recurrent superimposition gives rise to the two-cluster states detected by $`\zeta `$ can be qualitatively understood in terms of previous results for coupled oscillators with disordered interactions. Given an oscillator whose internal variable evolves with frequency $`\mathrm{\Omega }_i`$, the most persistent interactions affecting the dynamics of its phase $`\varphi _i`$ are due to oscillators with similar frequencies, $`\mathrm{\Omega }_j\mathrm{\Omega }_i`$. In this case, in fact, the interaction weights $`J_{ij}`$ are almost constant over very long time scales. These weights, however, are still different from each other, and can take positive and negative values. If the effect of frustration is moderate, the oscillators become entrained into two anti-phase groups of similar sizes d622 . The occurrence of this phenomenon requires that the oscillator phases are not pinned to the internal variables, and that a consistent population with similar frequencies $`\mathrm{\Omega }_i`$ actually exists to trigger condensation. Such conditions are met by those oscillators of subensemble II with either positive or negative frequencies, close to the limit $`|\mathrm{\Omega }_i|=\mu /2`$, where the population density is larger. Note that increasing the frequency dispersion $`\sigma `$ contributes in two cooperating ways to the growth of subensemble II, and thus to the formation of clusters: it implies, at the same time, a larger population for high frequencies and—as $`\mu `$ decreases—a population depletion in subensemble I. As discussed above, however, clustering is delayed by increasingly long transients as $`\sigma `$ approaches the critical point. In the limit $`N\mathrm{}`$, the evolution freezes just beyond the critical dispersion, and cluster formation is thus suppressed.

It is possible to give an approximate mathematical description of the motion of the four clusters, assuming that one of them is centered around a phase $`\varphi ^{}(t)`$ and is formed by oscillators whose internal frequencies are, on the average, $`\mathrm{\Omega }^{}`$. The phases and frequencies of the other three clusters can be derived from symmetry considerations. If the total fraction of the population involved in clusters is $`n^{}`$, the approximate equation of motion for the cluster phase reads

$$2\dot{\varphi }^{}=\mu \mathrm{sin}(\varphi ^{}\mathrm{\Omega }^{}t)n^{}\mathrm{cos}(2\mathrm{\Omega }^{}t)\mathrm{sin}(2\varphi ^{}).$$
(9)
The two terms in the right-hand side correspond, respectively, to the contributions of subensembles I and II. They have competing dynamical effects, and their relative importance changes as $`\sigma `$ and $`n^{}`$ grow while $`\mu `$ decreases. The first term favors the development of correlations between $`\varphi ^{}`$ and the internal variable $`\mathrm{\Omega }^{}t`$. If, on the other hand, the second term prevails, the phase difference between clusters of opposite internal frequencies, $`2\varphi ^{}`$, spends most of the time close to its two fixed points, $`2\varphi ^{}=0`$ or $`\pi `$, with rapid transitions between them when the sign of $`\mathrm{cos}(2\mathrm{\Omega }^{}t)`$ changes. These fixed points correspond, precisely, to the mutual crossing of clusters of opposite frequencies that give origin to the recurrent two-cluster state.

Preliminary numerical results show that dynamical clustering occurs also in ensembles of particles moving in space and coupled through finite-range interactions, with autonomous internal variables as those considered above. In more complex models of interacting dynamical elements with evolving internal states, dynamical clustering may be enhanced if the internal variables themselves undergo segregation into coherent groups. This can be achieved, for example, if interactions between internal states are allowed. Numerical realizations suggest that this is the case when the variables $`\theta _i`$ of the model studied above are coupled through higher-harmonics interaction functions okuda , which give rise to dynamical clustering if noise is added to the internal evolution. The same effect has recently been illustrated for ensembles of interacting particles moving in space, where the internal dynamics of each element is represented by a chaotic map whose evolution is coupled to those of neighboring particles shk . In real systems—specifically, in biological populations—the interaction between internal variables is expected to play an important role when such mechanisms as imitation and social influence are in action.

D. H. Z. acknowledges financial support from the Alexander von Humboldt Stiftung, Germany, and is grateful to the Fritz Haber Institut der Max Planck Gesellschaft, Berlin, for hospitality.


format: text

text:
# Quantum mechanics of the two-dimensional circular billiard plus baffle system and half-integral angular momentum

## Abstract

We examine the quantum mechanical eigensolutions of the two-dimensional infinite well or quantum billiard system consisting of a circular boundary with an infinite barrier or baffle along a radius. Because of the change in boundary conditions, this system includes quantized angular momentum values corresponding to half-integral multiples of $`\mathrm{}/2`$. We discuss the resulting energy eigenvalue spectrum and visualize some of the novel energy eigenstates found in this system. We also discuss the density of energy eigenvalues, $`N(E)`$, comparing this system to the standard circular well. These two billiard geometries have the same area (A=$`\pi R^2`$), but different perimeters ($`P=2\pi R`$ versus $`(2\pi +2)R`$), and we compare both cases to fits of $`N(E)`$ which make use of purely geometric arguments involving only $`A`$ and $`P`$. We also point out connections between the angular solutions of this system and the familiar pedagogical example of the one-dimensional infinite well plus $`\delta `$-function potential.

1. Introduction

Two-dimensional quantum systems can offer a range of interesting new features compared to more familiar 1D systems, including less trivial implementation of boundary conditions leading to energy quantization as well as more sophisticated connections between conserved quantum numbers and symmetries of the system. They can provide model systems to probe the connections between classical and quantum chaos, as in billiard geometries, and they also find proverbial ’real-life’ experimental realizations in a number of surface systems, such as atomic corrals. Two-dimensional systems can also provide the first chance to utilize important semi-classical methods such as periodic orbit theory gutzwiller ; brack , where the detailed structure (specifically the oscillatory component) of the energy level density, $`N(E)`$, can be very directly connected to the closed classical trajectories in the system. Finally, students may also find connections between solutions of the 2D Schrödinger equation for infinite well or billiard systems and more familiar solutions of wave equations for membranes (drumheads) from mathematical physics courses.

Such 2D quantum infinite well or billiard systems provide practice in the practical application of boundary conditions on the quantum wavefunction arising from purely geometric considerations and a number of example of two-dimensional infinite well systems have been discussed in the literature. The 2D infinite square well liboff  can be studied as an example of the product of two separable one-dimensional problems, but also in the context of energy level degeneracy, either due to obvious symmetries or for more subtle reasons accidental\_degeneracy . The problem of a $`45^{}45^{}90^{}`$ isosceles triangle billiard isoceles  can be easily solved using linear combination solutions of the 2D square well, while the problem of the equilateral ($`60^{}60^{}60^{}`$) triangle billiard ’footprint’ has been solved in closed form by a number of authors canadian robinett\_annals , using a variety of methods.

The circular infinite well, defined by the potential energy function

$$V(r,\theta )=V(r)=\{\begin{array}{cc}0\hfill & \text{for }r<R\hfill \\ \mathrm{}\hfill & \text{for }rR\hfill \end{array}$$
(1)
is equivalent to the problem of a classically vibrating circular membrane mathews , but has also been considered in detail in the quantum mechanical case robinett\_circular\_well . This system can provide students with their first introduction to quantized angular momentum values in a simple context and is the starting point for our discussion.

In this note, we will extend the analysis of the ’standard’ circular well (shown in Fig. 1(a)) to a similar infinite well or quantum billiard system consisting of the same circular ’footprint’, but with the addition of an infinitely narrow and infinitely high wall along the $`\theta =0`$ direction, which we will call a baffle, as shown in Fig. 1(c). The presence of this additional constraint changes the application of the boundary conditions, and hence the allowed quantum numbers, both the quantized values of angular momentum, as well as the energy eigenvalue spectrum. In the next section, we will review the solutions of the full circular well, as well as the related special case of the ’half-circular’ well (as shown in Fig. 1(b)), and then use these results to derive the quantum eigensolutions of the circular well plus baffle case. We find that half-integral values (in units of $`\mathrm{}`$) of quantized angular momentum are now allowed and we also visualize some of the novel position-space eigenfunctions

In Sec. 3, we then calculate the distributions of energy levels, $`N(E)`$, for the full-circular well and circular well plus baffle case, comparing both results to predictions made using purely geometric arguments requiring only knowledge of the area and perimeter of the two-dimensional systems, finding excellent agreement with formulae familiar in the mathematical literature. We emphasize the importance and relative straightforwardness of such analyses for 2D billiard systems. We then discuss, in Sec. 4, the similarities between the circle plus baffle system and the familiar pedagogical case of the one-dimensional infinite well plus repulsive $`\delta `$-function potential gettys delta\_revival , which illustrates how the solutions of the full circle problem are continuously related to the new half-integral angular momentum states. We end by presenting our conclusions and suggesting additional exercises of this type, involving other novel 2D billiard systems which can be analyzed in terms of an energy level density.

2. Circular well and variations

We begin by reviewing the derivation of the solutions of the circular infinite well. The two-dimensional, time-independent, free-particle Schrödinger equation, in the relevant polar coordinates, is written in the form

$$\frac{\mathrm{}^2}{2\mu }\left(\frac{^2}{r^2}+\frac{1}{r}\frac{}{r}+\frac{1}{r^2}\frac{^2}{\theta ^2}\right)\psi (r,\theta )=E\psi (r,\theta ).$$
(2)
For notational convenience, we have labeled the particle mass as $`\mu `$, in order to avoid confusion with familiar quantum numbers. We assume a separable solution of the form $`\psi (r,\theta )=R(r)\mathrm{\Theta }(\theta )`$, and obtain the angular and radial equations

$$\frac{d^2\mathrm{\Theta }_{(m)}(\theta )}{d\theta ^2}=m^2\mathrm{\Theta }_{(m)}(\theta )$$
(3)
and

$$\frac{d^2R(r)}{dr^2}+\frac{1}{r}\frac{dR(r)}{dr}\frac{m^2}{r^2}R(r)=k^2R(r)$$
(4)
where $`k=\sqrt{2\mu E/\mathrm{}^2}`$. Initially one can think of $`m^2`$ as simply a separation constant to be determined.

The angular equation has (normalized) solutions of the form

$$\mathrm{\Theta }_{(m)}(\theta )=\frac{1}{\sqrt{2\pi }}e^{im\theta }$$
(5)
and the requirement that the solutions be single-valued under redefinitions of the angular variable, namely $`\mathrm{\Theta }_{(m)}(\theta +2\pi )=\mathrm{\Theta }_{(m)}(\theta )`$, implies that

$$e^{im(\theta +2\pi )}=e^{im\theta }m=0,\pm 1,\pm 2,\pm 3,\mathrm{}$$
(6)
giving the familiar integral values of quantized angular momentum. Because of the central nature of the potential, angular momentum is also conserved, with the $`\mathrm{\Theta }_{(m)}(\theta )`$ being the eigenfunctions of $`\widehat{L}_z=(\mathrm{}/i)(/\theta )`$ with eigenvalues $`m\mathrm{}`$.

The resulting equation for the radial component, when written in terms of the variable $`z=kr`$, becomes,

$$\frac{d^2R(z)}{dz^2}+\frac{1}{z}\frac{dR(z)}{dz}+\left(1\frac{m^2}{z^2}\right)R(z)$$
(7)
which can be recognized as Bessel’s (cylindrical) equation. This has two linearly independent solutions for each value of $`|m|`$, the so-called regular $`J_m(z)`$ (well-behaved as $`z0`$) and the singular $`Y_m(z)`$ (divergent as $`z0`$) solutions. Only the well-behaved $`J_m(z)`$ are used here, with the energy eigenvalues determined by the boundary condition at the circular boundary, namely $`J_m(kR)=0`$. The resulting energy eigenvalues are then given by

$$E_{(m,n_r)}=\frac{\mathrm{}^2k_{(m,n_r)}^2}{2\mu }=\frac{\mathrm{}^2}{2\mu R^2}\left[z_{(m,n_r)}\right]^2$$
(8)
where $`z_{(m,n_r)}`$ is the $`n_r`$-th zero of the regular Bessel function $`J_m(z)`$. The number of radial nodes is then given by $`n_r1`$. The properly normalized radial wavefunctions are given by

$$𝒥_{(m,n_r)}(kr)N_{(m,n_r)}J_{(m,n_r)}(kr)\text{where}\left[N_{(m,n_r)}\right]^2_0^R\left[J_{(m,n_r)}(kr)\right]^2r𝑑r=1$$
(9)
with $`k=z_{(m,n_r)}/R`$.

The complete spectrum for the full circular well corresponds to one set of the $`m=0`$ solutions, but is doubly degenerate for each value of $`|m|>0`$, since both positive and negative values of $`m`$ give the same contributions; this corresponds physically to the equivalence of clockwise versus counter-clockwise motions. Some of the low-lying energy eigenvalues (solid lines) are illustrated in Fig. 2, ordered by angular momentum values.

The related case of the half-circular infinite well, shown in Fig. 1(b), can also be solved directly using these results. The angular eigenfunctions in Eqn. (5) can also be written in the (normalized) form

$$\stackrel{~}{\mathrm{\Theta }}_{(m)}=\{\begin{array}{cc}1/\sqrt{2\pi }\hfill & \text{for }m=0\hfill \\ \mathrm{cos}(m\theta )/\sqrt{\pi }\hfill & \text{for }m=1,2,3,\mathrm{}\hfill \\ \mathrm{sin}(m\theta )/\sqrt{\pi }\hfill & \text{for }m=1,2,3,\mathrm{}\hfill \end{array}$$
(10)
by taking appropriate linear combinations of the $`\mathrm{exp}(im\theta )`$ solutions, with the same pattern of degeneracies obvious (namely one $`m=0`$ state and doubly degenerate $`|m|0`$ states.)

The $`\mathrm{sin}(m\theta )`$ solutions will still satisfy the Schrödinger equation inside the half-circular well, as well as the boundary conditions that the wavefunction vanish along both the $`\theta =0`$ and $`\theta =\pi `$ directions. Thus, with a trivial change in normalization (an extra factor of $`\sqrt{2}`$ to account for the smaller ’footprint’), the $`\mathrm{sin}(m\theta )`$ solutions are also appropriate for the half-circular well, with the result that the energy eigenvalue spectrum for that problem consists of one copy of the integral $`m>0`$ energy eigenvalues in Fig. 2. The fact that the ’half-circle’ state has roughly half the states of the full circle problem will be elaborated upon in a detailed way in the discussion of energy eigenvalue density in Sec. 3.

Turning now to the case of the circular well plus baffle, as shown in Fig. 1(c), we can also make use of the forms in Eqn. (10), but we must also reinterpret the required boundary condition as imposed by the baffle. The vanishing of the angular wavefunctions along the $`\theta =0`$ line implies that the $`\mathrm{cos}(m\theta )`$ combinations (and the $`m=0`$ state) are not allowed, while the remaining $`\mathrm{sin}(m\theta )`$ solutions must now also vanish along the $`\theta =2\pi `$ line (the same infinite wall baffle) which implies that we must have

$$\mathrm{sin}(2\pi m)=0m=\frac{1}{2},1,\frac{3}{2},2,\frac{5}{2},\mathrm{}$$
(11)
and half-integral values of the angular momentum quantum number, $`m`$, are now allowed. The integral $`m=1,2,3,\mathrm{}`$ solutions from Eqn. (10) are still allowed, this time with the same normalization as for the full circular well, but there is a new class of angular (and hence radial solutions) characterized by $`m=(2j+1)/2=1/2,3/2,5/2,\mathrm{}`$.

As we will see in detail in Sec. 4, the angular eigenfunctions corresponding to half-integral values of $`m`$ correspond to the even functions of $`\theta `$ (about $`\theta =0`$), so we can also write these solutions in the form

$$\stackrel{~}{\mathrm{\Theta }}_{(m)}(\theta )=\frac{1}{\sqrt{\pi }}\mathrm{sin}(m|\theta |)\text{for }m=(2j+1)/2=1/2,3/2,\mathrm{}$$
(12)
if we choose to define them over the angular interval $`(\pi ,+\pi )`$. The first four lowest lying angular solutions corresponding to $`m=1/2,1,3/2`$, and $`2`$ are shown in Figs. 3(a), (b), (c), and (d) respectively, both over the interval $`(0,2\pi )`$ and $`(\pi ,+\pi )`$, illustrating both the symmetry as well as the ’cusp’ in $`\stackrel{~}{\mathrm{\Theta }}_{(m)}(\theta )`$ at $`\theta =0`$, induced by the infinite wall, for the $`m=(2j+1)/2`$ values. Because of the additional ’angular barrier’, the potential is no longer purely central, and angular momentum is not conserved. The $`\stackrel{~}{\mathrm{\Theta }}_{(m)}(\theta )`$ are, however, still eigenfunctions of $`\widehat{L}_z^2`$, with eigenvalues given by $`(m\mathrm{})^2`$; this corresponds classically to the fact that particles would ’rebound’ from the baffle, reversing direction, and hence the sign of angular momentum, but keeping the same magnitude of $`L_z`$.

The resulting energy eigenstates, $`\psi (r,\theta )`$, for the integral $`m`$ values are identical to the odd-parity states of the full circular well (since they already vanish completely along the entire $`\pm x`$ axis), but the (normalized) wavefunctions for the half-integral $`m=(2j+1)/2`$ angular momentum states states are now given by

$$\psi (r,\theta )=𝒥_{(m=j+1/2)}(kr)\stackrel{~}{\mathrm{\Theta }}_{(2j+1)/2}(\theta )\text{with}k=z_{(j+1/2,n_r)}/R.$$
(13)
We can use the well-known connection between the (cylindrical) Bessel functions, $`J_m(z)`$, and the spherical Bessel functions, $`j_m(z)`$ (obtained from the 3D version of the free particle Schrödinger equation gasiorowicz ) to write the half-integral solutions in the form

$$j_m(z)=\sqrt{\frac{\pi }{2z}}J_{m+1/2}(z)$$
(14)
Thus, half of the solutions of the 2D circular well plus baffle, those corresponding to half-integral angular momentum, can actually be described by integral angular momentum solutions, but of the corresponding 3D problem. For example, the two lowest $`m`$-value spherical Bessel functions are given by

$`j_0(z)`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}(z)}{z}}`$ (15)
$`j_1(z)`$ $`=`$ $`{\displaystyle \frac{\mathrm{sin}(z)}{z^2}}{\displaystyle \frac{\mathrm{cos}(z)}{z}}.`$ (16)
The energy eigenvalues for the $`m=1/2`$ case then correspond to the zeros of the 3D $`j_0(z)`$ function, namely $`z=n_r\pi `$, with $`n_r=1,2,\mathrm{}`$ with the result that the $`m=1/2`$ energies are given by the very simple formula

$$E_{(m=1/2,n_r)}=\frac{\mathrm{}^2\pi ^2n_r^2}{2\mu R^2}.$$
(17)
We illustrate the energy eigenvalues corresponding to half-integral values of $`m`$ in Fig. 2, where they nicely interpolate between the standard integral $`m`$ results. The number of states for the circular well plus baffle case is then, in some sense, exactly the same as that for the standard circular well, and we will in fact argue in Sec. 4, that the half-integral angular momentum states correspond very directly to the even-parity integral-$`m`$ states of the full circle, transforming into them as the effect of the baffle is slowly added.

We then plot in Figs. 4 and 5 the position-space probability density ($`|\psi (r,\theta )|^2`$) for the two lowest-lying energy eigenstates for the $`m=1/2`$ and $`m=3/2`$ states respectively, which can be compared to the more familiar looking (drumhead-like) case of $`m=1`$ shown in Fig. 5.

For all cases with $`m>1/2`$, because of the fairly obvious symmetries of the solutions, (with $`s=2m`$ ’spokes’ along nodal lines at $`\theta =2\pi /s`$), the expectation values of the position variables can be easily shown to vanish, namely, $`x=y=0`$ for all $`n_r`$ values when $`m>1/2`$. For the special case of $`m=1/2`$, however, the expectation value of $`\stackrel{~}{\mathrm{\Theta }}_{(1/2)}|\mathrm{cos}(\theta )|\stackrel{~}{\mathrm{\Theta }}_{(1/2)}_{m=1/2}=1/2`$ combines with the simple form of the radial solutions from Eqn. (15) to give $`𝒥_{(1/2)}(kr)|r|𝒥_{(1/2)}(kr)=R/2`$ to yield $`x_{(m=1/2,n_r)}=R/4`$ and this asymmetry is also obvious from Fig. 4.

3. Energy level density and completeness of the energy eigenstates

The calculation of the number of normal modes of solutions of the wave equation (scalar, electromagnetic, etc.) in 2-D or 3-D cavities is an important problem in a number of branches of physics, with many discussions building on early work by Weyl. Students are perhaps most familiar with this problem in the context of the derivation of the Planck formula for blackbody radiation for a simple 3D cubical cavity which is often included in the modern physics curriculum, but many interesting mathematical results exist for rather general two-dimensional shapes.

A standard reference text morse  on theoretical physics shows that the number of normal mode wavenumbers in the range $`(k,k+dk)`$ for a 2D shape of area $`A`$ and perimeter $`P`$ is given by

$$dN(k)=\left[\frac{A}{2\pi }k\frac{P}{4\pi }\right]dk$$
(18)
which upon integration gives

$$N(k)=\frac{A}{4\pi }k^2\frac{P}{4\pi }k$$
(19)
or, in the context of solutions of the Schrödinger equation of relevance here,

$$N(E)=\frac{A}{4\pi }\left(\frac{2\mu }{\mathrm{}^2}E\right)\frac{P}{4\pi }\sqrt{\frac{2\mu }{\mathrm{}^2}E}.$$
(20)
The ’experimental’ energy spectrum is, of course, a discontinuous ’staircase’ function of the form

$$N(E)=\underset{i}{}\theta (EE_i)$$
(21)
and so the Weyl-like result of Eqn. (20) will be an approximation to a smoothed out version of the ’data’.

This type of analysis can be extended monograph  to include an additional constant (and hence subleading) term which arises from the consideration of such geometrical effects as corners, curvature, and the connectivity of the 2D domain. Perhaps more importantly, the oscillatory behavior of $`N(E)`$ is the central theme of periodic orbit theory gutzwiller , brack  and simple pedagogical examples of this have been given robinett\_periodic\_orbit  for the square, circular well, and half circle ’footprints’.

The result of Eqn. (20) can be rather easily tested on familiar 2D billiard/infinite well systems such as the square or rectangular wells, as well as both the $`45^{}`$ isosceles and $`60^{}`$ equilateral triangles, all because of their extremely simple energy eigenvalue formulae (consisting of simple quadratic powers of two integral quantum numbers.) In fact, two studies of the equilateral triangle billiard canadian , berry  made use of this relation as a cross-check on the completeness of the energy eigenstates they derived, confirming that their calculated energy spectrum saturated the Weyl-like prediction for $`N(E)`$.

For the cases we consider here, we wish to compare the theoretically predicted spectrum for $`N(E)`$ against the geometric ’fit’ of Eqn. (20). From the differing ’footprints’ for the geometries we have considered, the relevant area ($`A`$) and perimeter ($`P`$) values are given by

| shape | $`A`$ | $`P`$ |
| --- | --- | --- |
| full circular well | $`\pi R^2`$ | $`2\pi R`$ |
| half circular well | $`\pi R^2/2`$ | $`(2+\pi )R`$ |
| well plus baffle | $`\pi R^2`$ | $`(2\pi +2)R`$ |

so that the full circular well and the well plus baffle have the same area, but differing perimeters (due to the ’intrusive’ baffle.) Note that the effective perimeter for the well plus baffle case is $`2R`$ larger than $`2\pi R`$ due to the fact that both the ’top’ and the ’bottom’ of the baffle are inside the well, each adding a factor of $`R`$.

We can easily collect (using numerical calculations of the required zeros of various order Bessel functions) large numbers of the low-lying energy eigenvalues for the cases we consider and we plot the resulting $`N(E)`$ versus $`E`$ using the first hundred or so lowest $`E`$ values for the three cases above in Fig. 7(a), as the three ’staircase’ functions shown there. The smooth curves are the predictions of Eqn. (20), while the two dotted curves are the predictions using only the area ($`A`$) terms which are clearly not a very good fit by themselves. The correspondence to the observed energy spectra is extremely good, and the clear difference between the full circle and half circle cases is obvious. A smaller region (shown as a dashed box) is enlarged on the right in Fig. 7(b) and there one can here easily see the nice distinction between the full circle case (upper data and dashed curve) and the circle plus baffle (lower data and solid curve) indicating that the effective values of $`A`$ and $`P`$ used above are correct. This case is especially interesting since it’s one of the few for which the ’footprint’ area is naturally the same, while the perimeter is different. The case of the 2D infinite square well ($`L\times L`$) or an equal area rectangular well ($`fL\times L/f`$) with different perimeters is another such simple example.

We wish to especially emphasize how this type of analysis can be rather easily applied to a wide variety of 2D billiard systems, perhaps with an eye towards introducing students to numerical methods, and eventually to the interesting topic of periodic orbit theory. We also note that the discussion can be made more quantitative by applying least square fits to the $`N(E)`$ versus $`E`$ data used in Fig. 7, using a functional form $`N(E)=aE+b\sqrt{E}`$, to obtain numerical values to compare directly to the results of Eqn. (20), with good agreement being obtained using only a few hundred low-lying states.

4. Relationship to 1D infinite well plus $`\delta `$ function problem

In order to better understand the structure of the half-integral angular wavefunctions derived above, and to make connection with a familiar one-dimensional quantum system, we wish to model the effect of introducing the baffle wall and continuously increasing its height. This approach is, then, very similar to the often discussed gettys delta\_revival  problem of a 1D infinite square well with a repulsive $`\delta `$-function potential placed at the center. We will briefly review the methodology and results of that problem, and then apply it to case of the circular baffle.

As an example of how the additional of a singular potential of arbitrary ’strength’ can change the energy level structure of a simple system, we consider a symmetric infinite well potential with walls at $`(L,+L)`$, which has even- and odd-parity energy eigenfunctions and eigenvalues given by

$`\psi _{(n)}^{(+)}(x)={\displaystyle \frac{1}{\sqrt{L}}}\mathrm{cos}\left({\displaystyle \frac{(2n1)\pi x}{2L}}\right)`$ with $`E_{(n)}^{(+)}={\displaystyle \frac{\mathrm{}^2\pi ^2(2n1)^2}{8\mu L^2}}`$ (22)
$`\psi _{(n)}^{()}(x)={\displaystyle \frac{1}{\sqrt{L}}}\mathrm{sin}\left({\displaystyle \frac{(2n)\pi x}{2L}}\right)`$ with $`E_{(n)}^{()}={\displaystyle \frac{\mathrm{}^2\pi ^2(2n)^2}{8\mu L^2}}`$ (23)
where once again we denote the particle mass by $`\mu `$.

We then introduce a singular, repulsive potential of the form $`V(x)=\lambda \delta (x)`$ to the center of the well and ask how the energy eigenvalues and eigenfunctions are changed. The odd parity $`\psi _{(n)}^{()}(x)`$, which have a node at $`x=0`$, are unaffected by the addition of $`V(x)`$, whatever the strength of the singularity, $`\lambda `$. The even parity solutions are changed and we can write a very general even solution in the form

$$\psi _{(\lambda )}^{(+)}(x)=\{\begin{array}{cc}A\mathrm{cos}(kx)+B\mathrm{sin}(kx)\hfill & \text{for }0x+L\hfill \\ A\mathrm{cos}(kx)B\mathrm{sin}(kx)\hfill & \text{for }Lx0\hfill \end{array}$$
(24)
The boundary conditions at either wall (where $`\psi `$ must vanish) and at the singularity (where $`\psi ^{}(0)`$ is discontinuous) are

$$\psi _{(\lambda )}(L)=\psi _{(\lambda )}(+L)=0\text{and}\psi _{(\lambda )}^{}(0_+)\psi _{(\lambda )}^{}(0_{})=\frac{2\mu \lambda }{\mathrm{}^2}\psi _{(\lambda )}(0)$$
(25)
and give the relations

$$A\mathrm{cos}(kl)+B\mathrm{sin}(kL)=0\text{and}2Bk=\frac{2\mu \lambda }{\mathrm{}^2}A$$
(26)
which combine to yield the energy eigenvalue condition

$$\lambda \left(\frac{\mu L}{\mathrm{}^2}\right)=\frac{kL\mathrm{cos}(kL)}{\mathrm{sin}(kL)}.$$
(27)
For the case of $`\lambda =0`$ (no additional $`\delta `$ perturbation), the solutions are given by $`kL=(2n1)\pi /2`$ and the standard $`\psi _{(n)}^{(+)}(x)`$ states are reproduced. In the limit of $`\lambda +\mathrm{}`$, however, a constant horizontal line of $`\lambda (\mu L/\mathrm{}^2)`$ cuts the right hand side of Eqn. (27) at $`kL=n\pi `$, so that the even energy eigenvalue solutions approach those of the (unchanged) odd parity solutions from below, namely $`E_{(\lambda )}^{(+)}(\lambda \mathrm{})E_{(n)}^{()}`$ and the energy spectrum is now doubly degenerate for each $`n`$ value.

For the case of the baffle added to the standard circular well, we will model the effect of continuously ’turning on’ the baffle by re-writing the angular Schrödinger equation in Eqn. (3) in the form

$$\frac{d^2\mathrm{\Theta }_{(m)}(\theta )}{d\theta ^2}+g\delta (\theta )=m^2\mathrm{\Theta }_{(m)}(\theta )$$
(28)
with the angles defined over the symmetric interval $`(\pi ,+\pi )`$. The odd parity (in $`\theta `$) angular $`\mathrm{sin}(m\theta )`$ solutions of Eqn. (10) are unaffected by the additional repulsive $`\delta `$ interaction, while we can write the even solutions, for arbitrary values of $`g`$, in a very similar form to that in Eqn. (24), namely

$$\mathrm{\Theta }_{(m)}^{(+)}(\theta )=\{\begin{array}{cc}A\mathrm{cos}(m\theta )+B\mathrm{sin}(m\theta )\hfill & \text{for }0\theta +\pi \hfill \\ A\mathrm{cos}(m\theta )B\mathrm{sin}(m\theta )\hfill & \text{for }\pi \theta 0\hfill \end{array}$$
(29)
The relevant boundary conditions are now on the continuity of $`\mathrm{\Theta }^{}(\theta )`$ at $`\theta =\pm \pi `$ and on the correct discontinuity at $`\theta =0`$ which give

$$2Bm\mathrm{cos}(m\pi )=2Am\mathrm{sin}(m\pi )\text{and}2mB=gA$$
(30)
respectively. These combine to give the condition for the quantized angular momentum quantum numbers, $`m`$, as

$$g=\frac{2m\mathrm{sin}(m\pi )}{\mathrm{cos}(m\pi )}f(m)$$
(31)
We plot the right-hand-side of Eqn. (31) in Fig. 8 and note that solutions of this eigenvalue problem once again correspond to the intersections of horizontal lines of constant $`g`$ with the various branches of $`f(m)`$. For $`g=0`$, namely the case of the standard circular well without the baffle, the intersections are at integral values of $`m`$, including $`m=0`$, and the even and odd angular wavefunctions are doubly degenerate (except for $`m=0`$) as noted above. For the case of $`g+\mathrm{}`$, the intersections arise at half-integral values of $`m`$, so that $`m=1/2,3/2,\mathrm{}`$ and the angular solutions are no longer degenerate. One can, in fact, write closed-form expressions for the (appropriately normalized) even parity angular eigenfunctions in the form

$$\stackrel{~}{\mathrm{\Theta }}_{(m)}^{(+)}(\theta )=\frac{\mathrm{cos}(m|\theta |\varphi _m)}{\sqrt{\pi (1+\mathrm{sin}(2m\pi )cos(2\varphi _m)/2m\pi )}}$$
(32)
where $`\varphi _m`$ is given by $`\mathrm{tan}(\varphi _m)=g/2m`$ and $`m`$ is given by the solutions of Eqn. (31). This form actually gives the correct normalization for both the $`g0`$ limit of integral angular momentum values (both $`m=0`$ and $`m>0`$) and the $`g+\mathrm{}`$ limit of half-integral values.

5. Conclusions and discussion

We have analyzed the circular well plus baffle system in detail, making use of the Weyl-like energy eigenvalue distribution of Eqn. (20) for two-dimensional quantum billiard systems. The system studied here can also be approached as the limiting case of a circular ’slice’ or ’wedge’ potential robinett\_slice . A circular infinite well, where the angle subtended by the billiard is given by $`\mathrm{\Phi }=(1+f)\pi `$, can interpolate between the half-circle case (when $`f=0`$) and the well circular well plus baffle case (when $`f=1`$). The boundary conditions ($`\mathrm{sin}(m\mathrm{\Phi })=0`$) imply that the angular momentum values must satisfy $`m=n/(1+f)`$ for integral $`n`$ and the one can still easily find the zeroes of $`J_{(m=n/(1+f))}(z)`$ using standard mathematical packages. Using this type of data, one can confirm that the energy eigenstate distribution, $`N(E)`$, still tracks the result of Eqn. (20), with continuously varying areas and perimeters now given by $`A_f=(1+f)\pi R^2`$ and $`P_f=(2+(1+f)\pi )R`$. One can also consider $`1<f<0`$ to handle such cases as the ’quarter-circle’ billiard and even smaller slices.

Another related case which can be analyzed using these methods is the annular infinite well or billiard robinett\_annular , robinett\_annular\_other , where a second concentric infinite barrier at $`R_{in}=fR`$ (with $`0<f<1`$) is added, and the particle is confined to the radial region $`R_{in}=fR<r<R`$. The angular solutions (with or without a baffle) are easily obtained as above, while the radial eigenfunctions are linear combinations of both cylindrical Bessel functions,

$$R(r)=\alpha J_m(kr)+\beta Y_m(kr),$$
(33)
where the divergent $`Y_m(kr)`$ is included since the singular $`r=0`$ point is explicitly excluded. Imposing the boundary condition that $`R(R_{in}=rR)=0`$ as well as vanishing at the outside boundary, gives the eigenvalue condition

$$J_m(kr)Y_m(fkR)J_m(fkR)Y_m(kR)=0$$
(34)
and one can once again evaluate $`N(E)`$ and compare it to Eqn. (20) using $`A_f=\pi R^2(1f^2)`$ and $`P_f=2\pi R(1+f)`$ with good agreement. Interested students can find a variety of relatively simple 2D quantum billiard geometries for which such an analysis is fairly straightforward.

Acknowledgments This work was supported in part by the National Science Foundation under Grant DUE-9950702.

Figure Captions

* Geometric ’footprints’ for the (a) circular infinite well (or billiard), (b) the half-circular well, and (c) the circular well with baffle, studied here.
* Energy spectrum versus quantized angular momentum values relevant for the circular, half-circular, and circular well with baffle cases. Energy values (in units of $`\mathrm{}^2/2\mu R^2`$) for integral values of angular momentum ($`m=0,1,2,\mathrm{}`$) are shown as solid, while results for half-integral ($`m=1/2,3/2,\mathrm{}`$) are shown as dashed. The spectrum for the full circular well consists of one copy of the $`m=0`$ values and two sets of the integral $`m>0`$ values, while for the half-circular well, one set of the integral $`m>0`$ values gives the entire spectrum. For the circular well plus baffle case, one set of both the integral and half-integral values with $`m>0`$ constitute the entire energy spectrum.
* Angular wavefunctions, $`\stackrel{~}{\mathrm{\Theta }}_{(m)}(\theta )`$ versus $`\theta `$, for half-integral and integral values of quantized angular momentum for (a) $`m=1/2`$, (b) $`m=1`$, (c) $`m=3/2`$, and (d) $`m=2`$, over the ranges $`(0,2\pi )`$ and $`(\pi ,+\pi )`$. The ’cusps’ at $`\theta =0`$ for the half-integral cases, induced by the infinite baffle along $`\theta =0`$, are evident.
* Normalized position-space probability densities, $`|\psi (r,\theta )|^2`$, for the two lowest-lying $`m=1/2`$ energy eigenstates ($`n_r=1`$ on top and $`n_r=2`$ at the bottom).
* Same as Fig. 4, but for the two lowest-lying $`m=3/2`$ states.
* Same as Fig. 4, but for the two lowest-lying $`m=1`$ states. These states are also solutions for the standard full circular well.
* Energy eigenvalue distribution, $`N(E)`$ versus $`E`$, for the full circle and circle plus baffle (top two curves in (a)) and half-circle case (bottom curve in (a)). The dashed and solid lines are the Weyl-like predictions of Eqn. (20), while the dotted lines are the result of only using the first ’area’ term. The dashed square part of (a) is expanded at the right in (b) to show the fine details and the clear differences between the full circle (top, dashed line) and circle plus baffle (bottom, solid line) cases.
* Plot of the ’eigenvalue’ condition of Eqn. (31) versus $`m`$ for even-parity angular wavefunctions. Horizontal lines of constant $`g`$ cut the curve at solutions of the eigenvalue condition for allowed values of $`m`$. For $`g=0`$ (corresponding to no baffle), the allowed eigenvalues are $`m=0,1,2,\mathrm{}`$ as expected, while for $`g+\mathrm{}`$, the intersections are at half-integral values given by $`m=1/2,3/2,\mathrm{}`$ corresponding to the circle plus baffle case.


text:
# Accuracy estimate for a relativistic Hamiltonian approach to bound-state problems in theories with asymptotic freedom

## I Introduction

This article describes a foretaste study of accuracy for a recently proposed Hamiltonian weak-coupling expansion procedure that in principle can start from a local asymptotically-free quantum field theory and produce sufficiently small relativistic effective Hamiltonian eigenvalue problems that may be soluble on a computer and yield the wave functions of bound states in that theory. The study is needed to determine the prospects of reaching a reasonable accuracy in the expansion since the strong coupling constant rises when the renormalization group scale is lowered toward the scale of binding mechanism af1 ; af2 . So far, successful approaches required different formulations of the theory at the bound-state and high-energy scales, such as Wilson’s lattice and Feynman’s diagrammatic techniques Wilson:1976kh ; match1 ; match2 ; match3 . They also used approximations such as the non-relativistic limit that applies in the case of heavy quarkonia charm1 ; charm2 ; charm3 ; charm4  and helps in constructing effective theories in analogy with QED CasLep ; krp ; nrqcd1 ; nrqcd2 .

The approach discussed here is the renormalization group procedure for effective particles (RGEP) RGEP  that stems from the application of the similarity renormalization group procedure similarity  to the light-front Hamiltonian of QCD long . In one and the same scheme, RGEP produces the asymptotically free coupling constant in the Hamiltonians for quarks and gluons g\_QCD , provides the conceptual framework for constructing the whole renormalized Poincare algebra in terms of the creation and annihilation operators for effective particles algebra , and leads to a simple first approximation in the case of heavy quarkonia ho . The procedure is boost invariant and raises a hope for connecting the constituent model for hadrons at rest Hagiwara  and the parton picture in the infinite momentum frame partons . The key ingredient of the procedure is the vertex form factor $`f`$ that multiplies all interaction terms. It falls off quickly to zero when the change of an invariant mass in an interaction vertex exceeds the width parameter $`\lambda `$. This width variable is also the renormalization-group evolution parameter. The effective particles that correspond to a small value of $`\lambda `$ cannot be copiously created because the form factor makes the interactions effectively weak even if the coupling constant becomes large. The internal structure of the effective particles at such small $`\lambda `$ may still be given by the parton-like picture in terms of constituents that correspond to a much larger value of $`\lambda Q`$ in the same flow, where $`Q`$ is the momentum scale of the external probe largep .

Qualitative accuracy studies of a similarity scheme were performed before modelafbs  using the elegant Wegner flow equation Wegner1 ; Wegner2 , which was invented independently for solving Hamiltonian eigenvalue problems in solid state physics. Unfortunately, it was found that the Wegner equation was not suitable for a straightforward weak-coupling expansion beyond second order. The coefficients of the expansion grew too fast and alternated in sign, which lead to erratic results for bound state energies with no signs of convergence. But the Wegner equation can be modified within the similarity scheme optimization  and the improved equation provides the benchmark here for estimating the accuracy prospects for the relativistic RGEP procedure. The comparison with the altered Wegner equation and the fact that improvements are also needed in condensed matter physics Mielke1 ; Mielke2 ; Mielke3 , imply that the effective particle approach may also find application outside QCD, i.e., wherever the dynamical couplings increase in the flow of Hamiltonians toward the region of physical interest.

The RGEP strategy that is tested here is to start with a regulated $`H`$ of the theory to be solved (this $`H`$ provides an initial condition for the differential equations of RGEP at $`\lambda =\mathrm{}`$), and to evaluate $`H_\lambda `$ with $`\lambda `$ on the order of a bound-state energy as a series in powers of a coupling constant $`g_\lambda `$. After evaluation of $`H_\lambda `$, one calculates its matrix elements in the basis defined by eigenstates of the Hamiltonian $`H_{0\lambda }`$, which is obtaind from $`H_\lambda `$ by setting $`g_\lambda =0`$; $`H_{0\lambda }|n=E_n|n`$ for all values of the label $`n`$. Suppose that the labels are ordered so that $`E_m<E_n`$ implies $`m<n`$ and the initial Hamiltonian is regulated by forcing matrix elements $`H_{mn}=m|H|n`$ to vanish unless $`Mm,nN`$ with certain ultraviolet cutoff number $`N`$ and infrared cutoff $`M`$. The RGEP procedure is designed in such a way that the matrix elements $`m|H_\lambda |n`$ quickly tend to 0 when $`|E_mE_n|`$ grows above $`\lambda `$, see Fig. 1. The next step is to focus attention on the window $`W_\lambda `$ of matrix elements of $`H_\lambda `$ among the basis states that have energies similar to the energy scale $`E`$ of the physical problem at hand, i.e.,

$`W_{\lambda mn}=m|H_\lambda |n,`$ (1)
with $`\stackrel{~}{M}m,n\stackrel{~}{N}`$ and $`E_{\stackrel{~}{M}}EE_{\stackrel{~}{N}}`$. Since only states within the width $`\lambda `$ on the energy scale can directly interact with each other, states that differ in energy from $`E`$ by much more than $`\lambda `$ are usually not important modelafbs ; optimization . They can be important as long as the coupling strength can overcome the difference in energy, but it does not matter here because the coupling constant is assumed to not grow to large values. The next step is to solve the non-perturbative eigenvalue equation for the matrix $`W_{\lambda mn}`$ by diagonalizing it on a computer. The middle-size eigenvalues of $`W`$ are expected to be close to the exact solutions with accuracy that depends on many factors in the procedure. These dependencies need to be estimated. The main question addressed here concerns the accuracy of the weak-coupling expansion for $`W_\lambda `$.

An exact RGEP procedure provides $`W_{\lambda mn}`$ whose structure depends on $`\lambda `$ but the eigenvalues of the middle size in the window spectrum do not. Becasue the couplings to the states outside the window are ignored, the eigenvalues of sizes at the limits of the window spectrum cannot be accurate even if the window is calculated exactly. Now, when one calculates $`W_\lambda `$ in the expansion in powers of $`g_\lambda `$ to some order, and extrapolates the result to $`g_\lambda 1`$, considerable errors can ensue because of the missing terms in the series. Moreover, one never knows the right value of $`g_\lambda `$ at given $`\lambda `$ from the theory. Therefore, one must fit $`g_\lambda `$ to the bound-state observables and perform consistency checks for a whole set of them. The smaller is $`\lambda `$ the larger is $`g_\lambda `$ and the more significant are the errors of perturbation theory in evaluation of $`W_\lambda `$. But at the same time the larger is $`\lambda `$, the larger must be the range $`[\stackrel{~}{M},\stackrel{~}{N}]`$ of basis states needed in the computer diagonalization of $`W_\lambda `$. A compromise must be made and the critical question is how close one can come to a true solution using RGEP equations. The question is essential for the prospects of applying RGEP to QCD. A well-known asymptotically free matrix model is used here to find out what level of accuracy can in principle be expected. One should stress at this point that even if the test gave a promising result in the model, the utility of the procedure would remain only tentative until the actual calculations in realistic theories are performed and display signs of stability as functions of the order of the expansion, size of $`\lambda `$, variation of the window, and other more specific features of the eigenvalue problems at hand.

Section II defines the test model and quotes equations used for calculating $`H_\lambda `$ in the benchmark of the altered Wegner’s flow and in the RGEP case. Section III describes results obtained using six different ways of fitting the coupling constant to the exact spectrum. These ways are labeled throughout the paper by letters A, B, C, D, E, and F. Each of these versions is studied in six successive orders of the weak-coupling expansion. Section IV provides a brief conclusion. Appendix A contains the generic analytic RGEP formula for 4th order calculations that applies to arbitrary $`H`$, and Appendix B provides the numbers that illustrate in detail what happens in the expansions including from one to six orders.

## II Model

The matrix elements of the model Hamiltonian used here for estimating the accuracy of RGEP are modelafbs ; optimization

$$H_{mn}=E_n\delta _{mn}g\sqrt{E_mE_n},$$
(2)
where $`E_n=b^n`$, $`b>1`$, and $`n`$ is an integer, with a convention that the energy equal 1 corresponds to 1 GeV. The model diverges and needs to be regulated by the ultraviolet cutoff $`\mathrm{\Lambda }=b^N`$ that limits the allowed energies $`E_n`$ from above. Also a low energy cutoff, $`b^M`$, is introduced for making the Hamiltonian matrix finite and enable exact computations, but the lower bound is of no physical consequence for the results reported here. Thus, the subscripts in Eq. (2) are limited to the range $`[M,N]`$, $`M`$ being large negative and $`N`$ large positive. The ultraviolet renormalizability of the model, its asymptotic freedom, and its lack of sensitivity to the infrared cutoff, were described in modelafbs ; optimization .

Two values of the cutoff $`N`$ and the corresponding coupling constant $`g=g_N`$ are used in this study of RGEP: $`g_{16}=0.060600631`$ and $`g_{20}=0.04878048667`$. The two cutoffs are introduced to verify the accuracy of renormalizability in the RGEP scheme. It is known that for such large values of $`N`$ the effective dynamics with $`\lambda 1`$ is practically independent of $`N`$ in the case of altered Wegner’s equation, and one can verify how well the RGEP approach satisfies this condition. At the same time, this condition puts constraints on the range of changes that one can consider in varying the RGEP generator without interference with renormalizability (see below). The coupling constants $`g_{16}`$ amd $`g_{20}`$ are fitted to obtain the exact bound-state eigenvalue $`E=1`$ with 8 digits of accuracy for $`b=2`$ and $`M=21`$, as in optimization . With these choices, the Hamiltonian $`H`$ is a $`(NM+1)`$ $`\times `$ $`(NM+1)`$ matrix with $`NM`$ positive eigenvalues, and one negative, spanning the range of energies between $`E_M0.5`$ KeV and $`E_N65`$ or 1000 TeV.

The similarity renormalization group procedure for Hamiltonians that leads to the altered Wegner equation Wegner1 ; optimization  and provides the benchmark here, can be written in the differentail form as

$$\frac{d}{d\lambda }=[F\{\},],$$
(3)
with the initial condition $`=H`$ when $`\lambda =\mathrm{}`$. The initial condition contains counterterms but the similarity analysis showed similarity  that the structure of the counterterms is simple and the presence of them is equivalent for large $`N`$ to making $`g`$ in $`H`$ depend on $`N`$. This is precisely what is done by the fitting mentioned earlier that guarantees that the eigenvalue $`E=1`$ stays unchanged for different $`N`$s. All small eigenvalues are then independent of $`N`$. The generator of the similarity transformation can be written as ($`𝒟_m=_{mm}`$)

$$m|F\{\}|n=h_{mn}(𝒟_m𝒟_n)_{mn}.$$
(4)
Different choices of $`h_{mn}`$ lead to different matrix elements of the renormalized Hamiltonians. Assuming

$$h_{mn}=\varphi _{mn}\frac{ds}{d\lambda },$$
(5)
one obtains Wegner’s equation when $`\varphi _{mn}1`$, and $`s=1/\lambda ^2`$ plays the role of the original Wegner parameter l Wegner1 . The altered Wegner equation has optimization

$`\varphi _{mn}={\displaystyle \frac{1}{1+|mn|}},`$ (6)
and this new equation is referred to as the benchmark.

In the plain perturbative RGEP procedure, the matrix elements of $`H_\lambda `$ in the effective basis states associated with the scale $`\lambda `$, are obtained from the matrix elements of an auxiliary Hamiltonian $``$ in the initial basis RGEP . The structure of $``$ is given by

$`=_0+f𝒢_I.`$ (7)
The Hamiltonian $`_0`$ is equal to the initial $`H_0`$, which is the free part of $`H`$, i.e., $`H`$ with $`g=0`$. $`f`$ denotes the form factor that can be written in all matrix elements using eigenvalues of $`_0`$, and reads

$`f_{mn}=\mathrm{exp}\left[\varphi _{mn}{\displaystyle \frac{(E_mE_n)^2}{\lambda ^2}}\right].`$ (8)
The RGEP equation for $`𝒢_I`$ is

$$\frac{d𝒢_I}{d\lambda }=[f𝒢_I,\left\{\frac{d}{d\lambda }(1f)𝒢_I\right\}],$$
(9)
where the curly bracket around an operator has the following meaning in terms of the matrix elements,

$$\left\{A\right\}_{mn}=\frac{A_{mn}}{E_nE_m}.$$
(10)
A priori, the optimal choice for $`\varphi _{mn}`$ in the RGEP could be different from the one that optimized the benchmark optimization . But we have studied various factors $`\varphi _{mn}`$ a la optimization  and found that $`c1`$ is also the best choice to make in RGEP, for similar reasons. In addition, it is useful for the test that these factors are made equal since then the first order calculations give identical results in both approaches. The factor $`\varphi _{mn}`$ is included in the analytic 4th order RGEP formulae provided in Appendix A.

The RGEP Eq. (9) cannot be integrated exactly on a computer as easily as Wegner’s equation can, because it contains the derivative of $`𝒢_I`$ on its right-hand side. One has to solve a complex linear problem to extract $`𝒢_I^{}`$. But the altered Wegner equation provides a perturbative benchmark pattern that is known to approximate an exact solution well and one can estimate the accuracy of RGEP by comparison.

In the perturbative evaluation of $`H_\lambda `$, the RGEP calculus is free from the problem of extracting $`𝒢_I^{}`$ because the derivative is computed order by order and all terms needed on the right-hand side are known at each successive order from the lower order results. In fact, the new procedure is designed for a perturbative approach. It is simpler than the Wegner case in the sense that the form factors guarantee the band structure in perturbation theory so that this structure does not have to be recovered from and controlled in the evolution of specific matrix elements. Small energy denominators that might otherwise lead to infrared singularities are excluded by design of RGEP. It also does not generate any terms with inverse powers of $`\lambda `$ in the coefficients of products of the interaction Hamiltonian (see Appendix A and Eq. (3.2) in optimization ). Such terms can in principle lead to a variety of mixing effects in the evolution of $`H_\lambda `$, when one uses the Wegner equation. On top of these purely quantum mechanical features, the perturbative RGEP calculus is capable of respecting seven Poincare symmetries in an economic way and has a potential to obtain the remaining three algebra . At the same time, it preserves the cluster decomposition property Weinberg  in the effective interactions. All these features are desired for the description of relativistic particles using theories with asymptotic freedom. But the accuracy of the weak-coupling expansion for windows $`W_\lambda `$ in RGEP must be measured against the benchmark to estimate the cost of the apparent advantages in terms of precision.

## III Accuracy of RGEP

The results of weak-coupling expansions up to the first-order terms in the RGEP and altered Wegner equation are identical, and read

$`_{mn}`$ $`=`$ $`E_m\delta _{mn}g_\lambda \sqrt{E_mE_n}f_{mn}.`$ (11)
The form factor $`f_{mn}`$ causes that the interaction Hamiltonian matrix is narrow on the energy scale and has width $`\lambda `$, see Figure 1. It is clear that in the case of Eq. (11) the effective coupling constant can be extracted from the matrix $``$ using the formula

$$g_\lambda =1\frac{_{M,M}(\lambda )}{E_M},$$
(12)
which is analogous to the Thomson limit in QED when $`M`$ is large and negative. The same Eq. (12) is used for defining $`g_\lambda `$ also in higher order calculations.

In the weak-coupling expansion in powers of $`g_\lambda `$ to order $`k`$ (see e.g. optimization ),

$`=_0+g_\lambda _1+g_\lambda ^2_2+\mathrm{}+g_\lambda ^k_k.`$ (13)
Appendix A contains analytic expressions for $`_k`$ with $`k=1`$, 2, 3, and 4. The right value of $`g_\lambda `$ can be found by solving the flow equations with the initial condition $`g_{\mathrm{}}=g_N`$ exactly. Then, one can check the accuracy of a perturbatively computed $`H_\lambda `$ by diagonalizing it for the exact value of $`g_\lambda `$, and by comparing the resulting bound-state eigenvalue with $`E=1`$ that was secured to exist by the initial choice of $`g_N`$ modelafbs ; optimization . The accuracy test for RGEP is carried out differently because the exact solution of Eq. (9) is not known. This situation is analogous to QCD where one can use perturbation theory to calculate $`H_\lambda `$ but an exact value of $`g_\lambda `$ is not available as a function of $`g_N`$. For the purpose of the accuracy test, the exact spectrum of the model is treated as experimental data. An approximate value of $`g_\lambda `$ is found by fitting some eigenvalue of the perturbatively calculated $`W_\lambda `$ to the data, or by performing a fit for a whole group of eigenvalues. Then one checks how well the bound-state eigenvalue $`E=1`$ is reproduced using the best-fit value for $`g_\lambda `$. In principle, one could fit $`g_\lambda `$ at one value of $`\lambda `$ that is most convenient for that purpose, evolve this value using RGEP to the new $`\lambda `$ that is most suitable for the bound-state calculation, and then compare the calculated spectrum with data g1 ; match3 ; g2 ; g3 . In fact, the model used here can be used for testing accuracy of such procedures in a comprehensive way. This type of tests may help in narrowing the current spread of estimates for $`\alpha _s`$ that come from various sources Hagiwara , by distinguishing theoretical procedures of least ambiguity. But the accuracy of RGEP is checked here using one and the same scale for fitting $`g_\lambda `$ and calculating the bound-state energy, for simplicity. The scale we choose is $`\lambda =2`$.

There exist infinitely many options for how one can fit $`g_\lambda `$ so that the spectrum of the window $`W_\lambda `$ at $`\lambda =2`$ approximates the exact one as closely as possible. We display results for six options that illustrate the dependence of results on such choices, labeled by A, B, C, D, E, and F. All methods are based on the minimalization of certain function $`K(g_\lambda )`$. We use

$$K_r(g_\lambda )=Z\underset{n}{}\left(\frac{v_n}{v_{en}}1\right)^2,$$
(14)
and

$$K_s(g_\lambda )=Z\underset{n}{}\left(\frac{v_nv_{n+1}}{v_{en}v_{en+1}}1\right)^2,$$
(15)
where $`v_{en}`$ is the exact eigenvalue of number $`n`$, and $`v_n`$ is the corresponding eigenvalue of $`W_\lambda `$ when $`W_\lambda `$ is derived in a given order $`k`$. $`Z`$ is the normalization constant equal to the inverse of the number of terms in the sum. $`Z`$ is not important in the minimalization of $`K(g_\lambda )`$ as function of $`g_\lambda `$. The subscript $`r`$ refers to the ratios of eigenvalues used in $`K_r`$, and the subscript $`s`$ refers to the splittings between the eigenvalues used in $`K_s`$. The six choices differ by which function $`K(g_\lambda )`$ is used and what is the range of summation over $`n`$ in Eqs. (14) or (15).

The eigenvalues of $`H`$ are numbered from $`M`$ to $`N`$ in the order in which they appear on the diagonal of $`H_\lambda `$ when $`\lambda 0`$ in the benchmark calculation. This numbering is also applied to the corresponding rows and columns of the matrix $`H_{\lambda mn}`$. The window $`W_\lambda `$ is always chosen to extend from $`\stackrel{~}{M}=8`$ to $`\stackrel{~}{N}=2`$ optimization , and it has 11 eigenvalues. These are numbered in the same order as for $`H_\lambda `$, with the bound-state eigenvalue having number $`n_0`$. We distinguish two eigenvalues that are closest in modulus to the bound-state energy $`E=1`$, one just smaller than 1, with number $`n_s`$, and one just larger than 1 with number $`n_l`$. The six fitting procedures are designated as follows:

$`\mathrm{A}`$ $``$ $`K_r,n=n_s,`$ (16)
$`\mathrm{B}`$ $``$ $`K_r,n=n_l,`$ (17)
$`\mathrm{C}`$ $``$ $`K_r,n\{n_s,n_l\},`$ (18)
$`\mathrm{D}`$ $``$ $`K_s,n=n_l,`$ (19)
$`\mathrm{E}`$ $``$ $`K_r,n[\stackrel{~}{M}+2,\stackrel{~}{N}2],nn_0,`$ (20)
$`\mathrm{F}`$ $``$ $`K_s,n[\stackrel{~}{M}+2,\stackrel{~}{N}2],nn_0.`$ (21)
In cases E and F, the two smallest and two largest eigenvalues are dropped because they are too much distorted by the edge of $`W_\lambda `$, as explained in Section I. The description of results obtained in the benchmark and tested approach in six successive orders of perturbation theory in each of these fitting schemes involves 72 results. The benchmark results are labeled ”Wegner” and the tested case is labeled ”RGEP.” Appendix B contains all pertinent numbers. An example is given below to help in reading the figures and tables in the appendix. Otherwise, only main features of the results are explained.

Fig. 2 corresponds to the case C of Eq. (18). The shape of $`K(g_\lambda )`$ clearly selects the value of $`g_\lambda `$ preferred by a given fitting procedure. Similar plots in other cases are given in Fig. 6 in Appendix B. The numbers that result from the fits for the coupling constants, along with the corresponding bound-state eigenvalues of $`W_\lambda `$, are also tabulated in Appendix B.

The summary of results for the coupling constants obtained in the fits is given in Fig. 3. The numbers 1 to 6 on the horizontal axis correspond to the order of perturbation theory in the evaluation of $`W_\lambda `$. The columns labeled A to F correspond to the algorithms given in Eqs. (16) to (21). The RGEP calculation equals the benchmark in the first order results. It is also visible that the fits of $`g_\lambda `$ in the benchmark case consistently reproduce the exact value of $`g_\lambda =0.2852`$ at $`\lambda =2`$ GeV optimization . The RGEP

displays similar stability and coalesces around 0.3. However, when the nearst neighbor level with energy larger than 1 is included in a fit as the only one (case B), or together with a nearest lower level using ratios of eigenvalues (case C), or together with the nearest lower level but using ratios of splittings between the two eigenvalues (case D), a visible variation in the fits occurs. This variation can be attributed to the lack of accuracy in the calculation of the nearst higher level, since in the cases E and F that include additional five lower levels, the higher level becoming much less significant, the fits resemble case A where the higher level is absent. This result suggests a rule for fits in future calculations that they should be focused on states with eigenvalues as far as possible from the bounds implied by the window choice and $`\lambda `$, in order to avoid the lack of convergence. Note that even in the benchmark case the 4th order calculation has to be corrected in the orders 5th and 6th to bring the accuracy into the few percent range around the exact value of $`g_\lambda `$. The same effect is observed in the bound-state eigenvalues themselves.

The summary of results for the bound-state eigenvalues is provided in Fig. 4 in a one-to-one correspondence to Fig. 3. The eigenvalues are obtained by diagonalization of windows $`W_\lambda `$ that are calculated in the six consecutive orders of perturbation theory (indicated on the horizontal axis in the individual columns), using the corresponding values of $`g_\lambda `$ from Fig. 3 (the numbers are given in Appendix B). It is visible that the RGEP procedure matches accuracy of the altered Wegner equation. One also sees the dramatic effect of the attempts to include the next higher level. The results clearly point out that both procedures should not be used for calculations of energy levels close to the size of $`\lambda `$. The fourth order calculations achieve accuracy on the order of 5%. This is encouraging, because orders 5 and 6 lead to even better results and one can expect that fits for $`g_\lambda `$ can be improved by focusing on the properties of low energy levels, including properties other than just eigenvalues. It is hard to imagine that such focus could lead to worsening of the accuracy.

The last question concerning accuracy of the RGEP method concerns renormalizability, which can be studied in analogy to optimization . However, one has to measure the sensitivity of the effective theory at the scale $`\lambda `$ to the cutoff $`\mathrm{\Lambda }=b^N`$ in a whole range of coupling constants, because the exact value of $`g_\lambda `$ is not known in the RGEP procedure. Using $`g_{16}`$ and $`g_{20}`$ from Section II, we calculate $`W_\lambda `$ starting from $`H`$ with the two different values of $`\mathrm{\Lambda }`$ for $`N=16`$ and $`N=20`$, and we compare the resulting matrix elements of $`W_\lambda `$ at $`\lambda =2`$ GeV. For $`b=2`$ the cutoff is changed by the factor 16, i.e. from about 65 TeV to about 1000 TeV. The divergence in the bare theory is logarithmic. The results could change at the rates implied by the change in the logarithm of $`\mathrm{\Lambda }`$, i.e., about 25%, but in the presence of a proper set of counterterms similarity  one expects no change to occur. Fig. 5 shows the measure of changes in $`W_\lambda `$ that are obtained here with only one counterterm, which amounts to the change of the coupling constant from $`g_{16}`$ to $`g_{20}`$. The function plotted in 5 is defined as

$`R(g_\lambda )={\displaystyle \underset{n,m=\stackrel{~}{M}}{\overset{\stackrel{~}{N}}{}}}\left[{\displaystyle \frac{W_{\lambda mn}(N=16)}{W_{\lambda mn}(N=20)}}\mathrm{\hspace{0.17em}1}\right]^2.`$ (22)
The results change with the order of perturbation theory used for evaluation of $`W_\lambda `$. The sensitivity of results to the change of $`\mathrm{\Lambda }`$ in 1st order is the same in RGEP as in the benchmark case optimization . Inclusion of higher orders, from 2 to 6, exhibits slight variations to the advantage of one or the other method, as shown in Fig. 5 (the solid line for the benchmark and the dashed one for RGEP). The last diagram in Fig. 5 shows the corresponding ratio $`R_W/R_R`$, and it demonstrates that the renormalizability condition is satisfied with the same accuracy in both approaches. A closer comparison requires that the benchmark case (or $`R_W`$) is taken for $`g_\lambda =0.2852`$, and RGEP (or $`R_R`$) for some $`g_\lambda 0.3\pm 0.005`$ (see Appendix B), but these corrections are negligible at the current stage of the analysis.

## IV Conclusion

The RGEP weak-coupling expansion achieves precision comparable with that of the best available benchmark method of the altered Wegner equation in the case of a simple matrix model for the Hamiltonian formulation of theories with asymptotic freedom and bound states. The few percent accuracy is reached by introducing a factor of $`\varphi `$ in Eq. (6), which is similar in both methods. Since the RGEP procedure is designed for application to relativistic quantum field theory, the model test calculation provides a comprehensive outline of steps that can be repeated in realistic cases, especially in the theory of quarks and gluons. On the other hand, knowing that Wegner’s equation is useful in condensed matter physics, one can expect that the altered Wegner equation and the RGEP procedure may also find applications in many-body theory. The model study indicates a possibility that the weak-coupling expansion may lead to a systematic approximation scheme despite the growth of the coupling when the characteristic scale $`\lambda `$ of the effective theories is lowered. This is indicated in columns E and F in Figs. 3 and 4, where one can see the convergence of the coupling constant to a stable value and the corresponding convergence of the bound-state energy to the exact result. One should also note that already second order calculations render effective window Hamiltonians $`W_\lambda `$ that can produce the bound-state eigenvalue with accuracy better than 10%.

The fact that both methods have similar accuracy suggests that they already show the range of calculational power that is available through a plain expansion in a single coupling constant such as the one defined in Eq. (12). In order to move beyond the 1% accuracy level one has to achieve better understanding of the structure of $`H_\lambda `$. For example, there exist terms in $`H_\lambda `$ with specific dependence on the eigenvalues of $`H_0`$ and $`\lambda `$, like $`(E_mE_n)^2/\lambda ^2`$ in the model. One may hope to understand how the coefficients of these terms depend on a more suitably defined coupling constant than the one given Eq. (12). It is conceivable that such understanding may further improve the weak-coupling expansion. Note also that none of the renormalization group universality features were so far explicitly employed in the plain expansion tested here. Such options for improvement may depend on a theory. One has to study specific theories using concrete versions of the RGEP procedure in order to identify the dominant terms and their universal behavior.

The accuracy study for the RGEP procedure shows also that one should be able to carry out similar tests for any other approach to the bound-state problem in asymptotically free theories. For such a test to become possible, the approach in question would have to be understood sufficiently well to determine the steps that that approach implies for handling of the initial $`H`$ in the model. But such understanding is demanded of most formulations of relativistic quantum theories for fundamental reasons and the accuracy tests of similar kind can be consider a challenge for any scheme intended to solve the bound-state problem in theories with asymptotic freedom.

## Appendix A RGEP of 4th order

The RGEP weak-coupling expansion for terms of order $`g^n`$, is written using terms of orders $`k<n`$,

$`𝒢_n^{^{}}`$ $`=`$ $`{\displaystyle \underset{k=1}{\overset{n1}{}}}[f𝒢_k,\left\{\left(1f\right)𝒢_{nk}\right\}^{^{}}].`$ (23)
If the interaction Hamiltonian $`H_I=HH_0`$ is proportional to $`g`$, the expansion in powers of $`g`$ is the same as the expansion in powers of $`H_I`$. If $`H_I`$ contains a polynomial in $`g`$ with operator coefficients, then the expansion in powers of $`H_I`$ is an intermediate step for obtaining an expansion for $`H_\lambda `$ in powers of $`g`$. Therefore, this appendix provides a generic set of coefficients for expansion of $`𝒢_k(\lambda )`$ in powers of $`H_I`$. The argument $`\lambda `$ is omitted in what follows. The procedure of rewriting the expansion in powers of the bare $`g`$ into the expansion in powers of $`g_\lambda `$ requires a definition of $`g_\lambda `$ in the structure of $`H_\lambda `$. This definition, in analogy to Eq. (12), provides then a series expression for $`g_\lambda (g)`$. This expression is inverted and substituted into the formal expansion of $`H_\lambda `$ in powers of the bare $`g`$, producing the desired expansion in powers of $`g_\lambda `$.

The first four terms in the perturbative expansion are written as

$`\left[𝒢_1\right]_{mn}`$ $`=`$ $`H_{mn},`$ (24)
$`\left[𝒢_2\right]_{mn}`$ $`=`$ $`a_{m1n}H_{m1}H_{1n},`$ (25)
$`\left[𝒢_3\right]_{mn}`$ $`=`$ $`b_{m12n}H_{m1}H_{12}H_{2n},`$ (26)
$`\left[𝒢_4\right]_{mn}`$ $`=`$ $`c_{m123n}H_{m1}H_{12}H_{23}H_{3n}.`$ (27)
It is uderstood that the indices 1, 2, and 3, are summed over the entire range available for them in the theory.

In order to write down expressions for the coefficients $`a`$, $`b`$, and $`c`$, we introduce a set of auxiliary symbols. Their meaning becomes sucessively self-evident when one decifers them in the order they are given here. The variables $`s`$ and $`t`$ have the meaning of the renormalization group parameter $`1/\lambda ^2`$.

$`E_{mn}`$ $`=`$ $`E_mE_n,`$ (28)
$`G_{mn}`$ $`=`$ $`\phi _{mn}\left(E_mE_n\right),`$ (29)
$`F_{mn}`$ $`=`$ $`\phi _{mn}\left(E_mE_n\right)^2,`$ (30)
$`G_{m1n}`$ $`=`$ $`G_{1m}+G_{1n},`$ (31)
$`F_{m\mathrm{}lkn}`$ $`=`$ $`F_{m\mathrm{}lk}+F_{kn},`$ (32)
$`A_{ml\mathrm{}kn}(t)`$ $`=`$ $`f_{ml}\mathrm{}f_{kn},`$ (33)
$`A_{m\mathrm{}l+k\mathrm{}n}(t)`$ $`=`$ $`A_{m\mathrm{}l}(t)A_{k\mathrm{}n}(t),`$ (34)
$`B_{m\mathrm{}k}(s)`$ $`=`$ $`{\displaystyle _0^s}𝑑tA_{m\mathrm{}k}(t).`$ (35)
In this notation one obtains

$`a_{m1n}`$ $`=G_{m1n}`$ $`B_{m1n}.`$ (36)
$`b_{m12n}`$ $`=`$ $`{\displaystyle \frac{G_{m12}}{E_{m2}}}\left[B_{m12n+m2}B_{m12n}\right]`$ (37)
$`+`$ $`{\displaystyle \frac{G_{n12}}{E_{n2}}}\left[B_{m21n+n2}B_{m21n}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{m12}}{F_{m12}}}\left[B_{m2n+m12}B_{m2n}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{n12}}{F_{n12}}}\left[B_{m2n+n12}B_{m2n}\right]`$
$`c_{m123n}`$ $`=`$ $`{\displaystyle \frac{G_{23m}}{E_{m1}E_{m2}}}\left[B_{m321n}B_{m321n+1m}B_{m321n+2m}+B_{m321n+1m2}\right]`$ (38)
$`+`$ $`{\displaystyle \frac{G_{23n}}{E_{n1}E_{n2}}}\left[B_{m123n}B_{m123n+1n}B_{m123n+2n}+B_{m123n+1n2}\right]`$
$`+`$ $`{\displaystyle \frac{G_{132}}{E_{12}E_{m1}}}\left[B_{m231n}B_{m231n+21}B_{m231n+m1}+B_{m231n+m12}\right]`$
$`+`$ $`{\displaystyle \frac{G_{132}}{E_{12}E_{n1}}}\left[B_{m132n}B_{m132n+21}B_{m132n+n1}+B_{m132n+n12}\right]`$
$`+`$ $`{\displaystyle \frac{G_{12m}G_{23m}}{F_{23m}E_{m1}}}\left[B_{m21n}B_{m21n+m1}B_{m21n+m32}+B_{m21n+1m32}\right]`$
$`+`$ $`{\displaystyle \frac{G_{12n}G_{23n}}{F_{23n}E_{n1}}}\left[B_{m12n}B_{m12n+n1}B_{m12n+n32}+B_{m12n+1n32}\right]`$
$`+`$ $`{\displaystyle \frac{G_{12m}G_{13n}}{F_{13n}E_{m1}}}\left[B_{m21n}B_{m21n+m1}B_{m21n+n31}+B_{m21n+m13n}\right]`$
$`+`$ $`{\displaystyle \frac{G_{12n}G_{13m}}{F_{13m}E_{n1}}}\left[B_{m12n}B_{m12n+n1}B_{m12n+m31}+B_{m12n+m31n}\right]`$
$`+`$ $`{\displaystyle \frac{G_{12m}G_{132}}{F_{132}E_{m1}}}\left[B_{m21n}B_{m21n+m1}B_{m21n+132}+B_{m21n+m132}\right]`$
$`+`$ $`{\displaystyle \frac{G_{12n}G_{132}}{F_{132}E_{n1}}}\left[B_{m12n}B_{m12n+n1}B_{m12n+132}+B_{m12n+n132}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{23m}}{E_{m2}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+m321}}{F_{m321}}}{\displaystyle \frac{B_{m1n}B_{m1n+m321+m2}}{F_{m321}+F_{2m}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{23n}}{E_{n2}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+n321}}{F_{n321}}}{\displaystyle \frac{B_{m1n}B_{m1n+n321+n2}}{F_{n321}+F_{2n}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{132}}{E_{12}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+m231}}{F_{m231}}}{\displaystyle \frac{B_{m1n}B_{m1n+m2132}}{F_{m231}+F_{21}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{132}}{E_{12}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+n231}}{F_{n231}}}{\displaystyle \frac{B_{m1n}B_{m1n+n2132}}{F_{n231}+F_{21}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{12m}G_{23m}}{F_{23m}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+m21}}{F_{m21}}}{\displaystyle \frac{B_{m1n}B_{m1n+m321+m2}}{F_{m21}+F_{m32}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{12n}G_{23n}}{F_{23n}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+n21}}{F_{n21}}}{\displaystyle \frac{B_{m1n}B_{m1n+n321+n2}}{F_{n21}+F_{n32}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{12m}G_{132}}{F_{132}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+m21}}{F_{m21}}}{\displaystyle \frac{B_{m1n}B_{m1n+m2312}}{F_{m21}+F_{231}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{12n}G_{132}}{F_{132}}}\left[{\displaystyle \frac{B_{m1n}B_{m1n+n21}}{F_{n21}}}{\displaystyle \frac{B_{m1n}B_{m1n+n2312}}{F_{n21}+F_{231}}}\right]`$
$`+`$ $`{\displaystyle \frac{G_{m1n}G_{12n}G_{13m}}{F_{12n}F_{13m}}}\left[B_{m1n}B_{m1n+n21}B_{m1n+m31}+B_{m1n+m312n}\right]`$
In relativistic applications, these formulae need a replacement of the differences of energies like $`E_mE_n`$ by the changes of invariant masses in the interaction vertices, and multiplication of $`G_{mn}`$ by the parent momenta $`P_{mn}^+`$ in the vertices, see e.g. g\_QCD ; ho .

## Appendix B Numerical details

Numerical integration of Eqs. (3) and (9) was performed using fourth-order Runge-Kutta procedure and an automated algorithm for generating expressions of order $`n`$ from expressions of orders $`k<n`$. The numerically calculated matrix elements $`H_{\lambda mn}`$ were checked using analytic formula from Appendix A.

The fits of the coupling constants $`g_\lambda `$ that were used in the perturbative evaluation of the windows $`W_\lambda `$, were obtained using the functions $`K(g_\lambda )`$ from Eqs. (14) and (15) that are plotted in Fig. 6. The corresponding numbers for the coupling constants are given in Table I.