text:
# Monte Carlo study of linear chain submonolayer structures. Application to Li/W(112)

## I Introduction

Chemisorption on metal surfaces has been studied for over 30 years naumov99 ; persson92 ; zangwill ; braun89 ; naumov77 ; tsong88 ; pfnuer459 ; pfnuer62 ; pfnuer460 . One of interesting problems concerns with structures and phase transitions in chemisorbed submonolayers. In particular alkaline, alkaline–earth, and rare-earth elements adsorbed on W(112) and Mo(112) form linear chain structures naumov99 ; pfnuer459 . In case of low coverages ($`\mathrm{\Theta }<0.5`$), the commensurate phases with chains normal to the troughs of the substrate are observed. The period of these structures changes from 2 to 9naumov99 . At higher coverages ($`\mathrm{\Theta }>0.5`$), incommensurate surface phases occur pfnuer62 ; pfnuer460 . It has been suggested naumov99 ; braun89  that formation of chain structures might be due to indirect interaction between adatoms via conduction electrons. Recently we proposed BLO  the lattice gas model to account for linear chain structures on W(112) and Mo(112). The model includes dipole–dipole interaction as well as indirect oscillatory interaction. As the interactions are of long range, the model is difficult to study theoretically. Preliminary studies performed in molecular field approximation confirmed that competition of dipole-dipole interaction and indirect interaction leads to formation of linear chain structures. Further extensive analysis of ground states in related effective one-dimensional model with infinite range of interactions LO  pointed out that competition between these two interactions is crucial in formation of linear chain structures. In very recent investigation of phase transitions in Li/Mo(112) and Sr/Mo(112), H. Pfnür et al. pfnuer459 ; pfnuer457 ; pfnuer67  have shown that the lattice gas model could be useful in description of long periodic phases observed at low coverages. Results of their investigation support an idea that surface electronic state are responsible for the oscillatory indirect interaction between adatoms.

It is intention of this paper to study phase diagrams by Monte Carlo simulations. This method was successfully applied to lattice gas models with short range interactions selke83 ; binder80 ; binder82 ; binder76 . On the other hand, one can expect more difficult problems in case of long range interactions, e.g., occurrence of complicated structures roelofs87 ; naumov92 ; bak82 ; sasaki94 . The paper is organized as follows. In section 2 the lattice gas model is described. Two difficulties connected with application of Monte Carlo simulations to this model are discussed in section 3. In section 4 phase diagrams are presented and section 5 contains discussion and conclusions.

## II The model

In the previous paper BLO  we introduced the lattice gas model of linear chain submonolayer structures on the (112) surface of W and Mo. Adsorption sites form the rectangular lattice (see Fig. 1) with $`a_x=a\sqrt{2}`$ and $`a_y=a\sqrt{3}/2`$, where $`a`$ denotes the lattice constant of the bcc elementary cell.

The model Hamiltonian within the grand canonical ensemble, is defined as

$$=\frac{1}{2}\underset{\genfrac{}{}{0pt}{}{i,j}{ji}}{}V_{ij}n_in_j\mu \underset{i}{}n_i\text{,}$$
(1)
where $`n_i`$ is the occupation variable at each adsorption site $`i`$ with $`n_i=1`$ if the $`i`$th site is occupied by an adatom and $`n_i=0`$, if not. The chemical potential is denoted by $`\mu `$ and $`V_{ij}=V(\stackrel{}{r_i}\stackrel{}{r_j})`$ is the pairwise interaction consisting of electrostatic and indirect interactions

$$V(\stackrel{}{r})=\frac{2d^2}{|\stackrel{}{r}|^3}+\{\begin{array}{cc}A\mathrm{cos}(2k_F|y|+\phi )\frac{1}{|y|}\delta (x,0)\hfill & \hfill \mathrm{for}y0\\ & \\ E_b\delta (|x|,a_x)\hfill & \hfill \mathrm{for}y=0\end{array}\text{.}$$
(2)
The first part of Eq. (2), $`\mathrm{\hspace{0.33em}2}d^2/|\stackrel{}{r}|^3`$, describes a repulsive dipole-dipole interaction between two adatoms with identical dipole moments $`d`$. The second term in Eq. (2) represents an asymptotic part of the indirect interaction between adatoms via conduction electrons of the substrate. This particular interaction is highly anisotropic and is closely related to the existence of nearly flattened segments of the Fermi surface being perpendicular to the $`[\overline{1}\overline{1}1]`$ axis directed along the atomic troughs of the substrate braun89 ; BLO . It is also possible that contribution to this interaction comes from quasi-one-dimensional surfaces states brenig95 . Very recently, G. Godzik et al.pfnuer67  suggested that oscillatory indirect interaction of this form is, in case of Sr/Mo(112), induced by charge density waves involving surfaces states. An amplitude $`A`$ and a phase $`\phi `$ can be treated as phenomenological parameters and $`k_F`$ denotes a wavevector of electrons at the Fermi surface. Finally, for $`y=0`$ there is an attractive indirect interaction $`E_b`$ between the nearest neighbours along the x-direction which is responsible for the chain formation.

In this work we choose the set of model parameters the same as proposed in the previous paper BLO  for adsorption of lithium on W(112) surface, i.e., $`d=1.5D`$, $`A=0.137eV`$Å, $`k_F=0.41\AA ^1`$, $`\phi =1.26\pi `$, $`E_b=0.05eV`$, and $`a=3.16`$Å. We will use a reduced chemical potential $`\mu /\mu _0`$ with $`\mu _0=k_BT`$, $`T=100K`$. Although the interaction described by Eq. (2) is long-range we introduce finite range of this interaction in Monte Carlo study. The dipole-dipole interaction is accounted for distances between adatoms smaller than $`R_m`$ and the indirect interaction along the y-direction for distances smaller than $`Y_m`$ (see Fig. 1).

The adatom coverage $`\theta `$ is defined as

$$\theta =\frac{1}{N}\underset{i}{}n_i,$$
(3)
where $`N`$ is the number of adsorption sites and angle brackets mean a thermodynamic average.

It is convenient to transcript the lattice gas model as Ising model by introducing spin variables $`s_i=12n_i`$

$$=\frac{1}{8}\underset{\genfrac{}{}{0pt}{}{i,j}{ji}}{}V_{ij}s_is_jH\underset{i}{}s_i\text{,}$$
(4)
where the magnetic field $`H`$ is expressed by the chemical and interaction potentials:

$$H=\frac{\mu }{2}+\frac{1}{4}\underset{j}{}V_{ij}.$$
(5)
Using Ising Hamiltonian, Eq. (4), we can perform simulation in the canonical ensemble which is equivalent to simulation of the lattice gas, Eq. (1), in the grand canonical ensemble. Of course, the former method is much easier.

## III Monte Carlo method

In Monte Carlo simulation we employ standard Metropolis algorithm to Ising Hamiltonian, Eq. (4), on finite rectangular lattice of $`L_x\times L_y`$ sites with periodic boundary conditions (PBC). Quantities such as magnetisation $`M`$, specific heat $`C_H`$, and energy $`E`$ are obtained as

$$M=\frac{1}{N}\underset{i}{}<s_i>,$$
(6)
$$C_H=\frac{1}{Nk_BT^2}\left(<^2><>^2\right),$$
(7)
$$E=\frac{1}{N}<>,$$
(8)
where $`N=L_xL_y`$ is the number of spins. We will use the relation between coverage $`\theta `$ and the magnetisation $`M`$ : $`\theta =(1M)/2`$. Since we are interested in linear chain structures we measure also magnetisation at each row of the lattice

$$M_j=\frac{1}{L_x}\underset{i}{}<s_{𝐫_i^{(j)}}>,$$
(9)
where

$$𝐫_i^{(j)}=(ia_x,ja_y),$$

Having $`M_j`$ one can calculate coverage at jth lattice row $`\theta _j=(1M_j)/2`$.

To study phase diagram we calculate the correlation function $`<s_js_l>`$ and the structure factor defined as selke83

$$S(q)=\frac{1}{N^2}\underset{j,l}{}<s_js_l>\mathrm{exp}(i𝐪(𝐫_j𝐫_l)),$$
(10)
where $`q`$ is a wavevector from reciprocal space. This function is used to recognize an ordered phase.

We encountered two difficulties in examination of phase diagrams by Monte Carlo method:

1. The dependence of the results on range of the indirect interaction.
2. The dependence of the resulting spin structure, at low temperatures, on an initial spin configuration.

In the previous paper BLO  we studied this model within the molecular field approximation (MFA) finding several ordered phases of $`(p\times 1)`$ type in phase diagrams. The $`(p\times 1)`$ structure at $`T=0`$, has every $`jp`$ row (for $`j=1,2,\mathrm{}\mathrm{},L_y/p`$) occupied by adatoms and remaining rows are empty. Therefore, coverage in this structure $`\theta =1/p`$ and the structure factor $`S(q)`$ reaches maximum at $`q=(0,2\pi /(pa_y))`$. The phase diagram constructed for the Li/W(112) system contains the $`(2\times 1)`$, $`(3\times 1)`$, $`(4\times 1)`$, $`(6\times 1)`$, and $`(9\times 1)`$ linear chain structures. Performing Monte Carlo simulations we observed that additional ordered structures can occur in a phase diagram. It may happen when the range of indirect interaction $`Y_m`$ is shorter than the linear lattice size along the y-direction $`L_y`$.

We start investigation of this problem by checking the influence of indirect interaction range on the ground states. As an example we present in Fig. 2 a part of ground-state diagram for constant range of interaction $`Y_m=72`$ and different values of the lattice size $`L_y`$. We see that only two phases, $`(2\times 1)`$ and $`(3\times 1)`$, are present in this part of phase diagram when the interaction range $`Y_m`$ is close to the lattice size ($`L_y<108`$). An additional phase denoted as ($`2a`$) appears for $`L_y>108`$. This is linear chain structure consisting of two domains of the $`(2\times 1)`$ phase and separated by two empty rows. Thus it has coverage $`\theta =(L_y2)/(2L_y)`$. It is worth noting that the $`(2\times 1)`$ phase has two domains: one with adsorbate chains placed in even lattice rows and the second with chains in odd rows. The configuration of the $`(2a)`$ phase at $`L_y=144`$ may be represented by sequence of distances between subsequent chains in the following way $`\{2_{35},3_1,2_{34},3_1\}`$, where $`l_j`$ denotes sequence of $`j`$ identical distances $`l`$. When range of interaction is nearly $`1/4`$ of $`L_y`$ then another phase ($`2b`$), containing four domain walls, occurs in the phase diagram (see Fig. 2). It looks like double ($`2a`$) phase, e.g., the ($`2b`$) phase can be represented by following sequences of distances $`\{2_{34},3_1,2_{34},3_1,2_{31},3_1,2_{27},3_1\}`$ for $`L_y=264`$. In a part of diagram corresponding to lower coverage ( not shown in Fig. 2) there are also phases containing walls and domains of $`(4\times 1)`$, and $`(6\times 1)`$ long periodic structures but detailed discussion will not be presented here. It worth to mention at this place, that the($`2a`$) phase can occur in the temperature phase diagram, even if the ($`2a`$) phase is not a ground state. In order to check the role of the boundary conditions, we performed additional calculation of the ground states using free edge boundary condition along the y-direction. The same defected phases where obtained but the number of domain walls was reduced by one. Above discussion shows that truncation of long-range indirect interaction in our model causes some technical difficulties in applying Monte Carlo method, e.g. it can be impossible to study finite-size scaling. On the other hand, the truncation of the interaction interpreted as screening effect might be important from physical point of view because it generates defected phases containing domain walls.

The second problem observed in our Monte Carlo calculation, dependence of resulting configuration on the initial spin configuration, is due to presence of short range attractive interaction between the nearest neighbours along the x-direction. When linear chain structure is formed at low temperature, it is difficult to change the period of such structure by changing temperature or chemical potential using single spin flip Metropolis algorithm. We study a phase diagram performing MC runs at constant magnetic field (or chemical potential) by cooling down or heating up the sample. The simulation start from high temperature in the cooling down process and after thermalization and measurement temperature is lowering the by $`\mathrm{\Delta }T`$ then such procedure is repeated till $`T=0`$. At each temperature, all quantities are measured after $`2\times 10^4`$ MC steps per spin during subsequent 5000 MC steps. The final spin configuration at $`T`$ is used as initial one at $`T\mathrm{\Delta }T`$. For some values of H the cooling down process finishes at $`T=0`$ with configuration having energy greater than energy of ground-state. Hence this way can produce metastable phases. To overcome this problem we use heating up process which starts at temperature $`\mathrm{\Delta }T`$ with ground-state configuration as the initial spin configuration. In subsequent temperature steps, $`T_i=T_{i1}+\mathrm{\Delta }T`$, the specific heat $`C_H`$ and the mean energy $`<>`$ are measured. This allow us to calculate the entropy

$$S(T)=_0^T\frac{C_H(T^{})}{T^{}}𝑑T^{},$$

by numerical integration. Then the free energy can be obtained as

$$F(T)=<>TS(T).$$

Repeating simulations at given $`H`$ with different ground-state configurations as initial configurations, we can find a phase with the lowest free energy at temperature T.

## IV Phase diagrams

To study phase diagrams by Monte Carlo simulations we used method described in the previous section, i.e., simulations at constant chemical potential (or magnetic field) with variable temperature. The phase transition between ordered phases is determined by examination of free energies of appropriate phases. Temperature of transition from disordered phase to an ordered one is determined from behaviour of the structure factor $`S(q)`$, specific heat $`C_H`$ or sublattice magnetisation $`M_j`$. Lattice sizes $`L_x`$, $`L_y`$ used in our simulations are rather small due to long-range interaction and large number of MC runs needed to construct phase diagrams. Most simulations were performed for $`20\times 72`$ lattice sites, but some calculations were performed on $`40\times 72`$ and $`20\times 144`$. Periodic boundary conditions were assumed, therefore $`L_y`$ was chosen as multiplicity of 36 in order to allow formation of ideal chain structures with period 2, 3, 4, 6, and 9, which exist in the ground-state phase diagram.

It has been shown in the previous section that defected ordered structure can occur when $`Y_m<L_y`$. We now discuss the thermal stability of such structure for two sets of parameters $`L_y`$ and $`Y_m`$. We see in Fig. 3a that the ($`2a`$) phase dominates the part of phase diagram corresponding to intermediate coverage and it pushes down the phase $`(2\times 1)`$. Hence, it is impossible to pass directly from the ordered $`(2\times 1)`$ phase to the disordered phase. More interesting case is shown in Fig. 3b where the ($`2a`$) phase occurs at high temperatures, although it is not stable at T=0. In this case, the ($`2a`$) phase is also placed between the $`(2\times 1)`$ phase and the disordered phase. We think, that this is important result because it shows that analysis of the ground states is not sufficient to find all ordered phases. The phase ($`2a`$) has two interesting features

1. The wavevector $`q=(0,2\pi (L_y2)/(2L_ya_y)`$, at which the structure factor has maximum, does not depend on temperature.
2. The sublattice order parameters vanish at some temperatures whereas the behaviour of $`S(q)`$ indicates the existence of ordered structure.

The first feature indicates that the phase ($`2a`$) is not an incommensurate phase. In order to understand the second feature we performed simulations with different Monte Carlo measurement times. If the time of measurement is not too long (e.g., $`5\times 10^4`$ MC time steps at $`T=300K`$) we observe (see Fig. 4a) that order in this structure changes spatially along y-direction.

The disorder appears near domain walls but as one moves form the wall the chain structure emerges with increasing ordering. On the other hand we found that domain walls and the whole structure are moving along the y direction during simulation. This is a reason why the sublattice order parameters vanish if averaging time is long enough (see Fig. 4 where results of simulations with different averaging time are presented) but the structure factor remains finite.

### IV.1 Application to Li/W(112)

Linear chain structures $`(2\times 1)`$, $`(3\times 1)`$, and $`(4\times 1)`$ were observed in the experiment of adsorption of lithium atoms on W(112) surface naumov73 . Using MFA method, we have shownBLO  that lattice gas model, Eq. (1), describes such sequence of structures. Now we present results of Monte Carlo simulations. To avoid occurrence of phases with domain walls of ($`2a`$) type, we chose the range of indirect interaction equal to the linear size of the lattice along y direction, i.e., $`Y_m=L_y`$. Most simulations were performed on $`20\times 72`$ lattice and sometimes the lattice of $`40\times 72`$ sites was used. The range of dipole-dipole interaction was $`R_m=12a_y`$. The ranges of these interactions used in previous MFA studies BLO  were both equal to $`36a_y`$. Therefore we also repeated the MFA calculations for present values of interaction ranges in order to compare the MFA results with Monte Carlo simulation results.

As previously, MFA calculation were performed for all structures of ($`p\times 1`$) type with $`p`$ up to 12. In Fig. 5 we plot results of MFA and MC methods in the (T,$`\mu `$) phase diagram. A low temperature part of the diagram does not depend on the method. However in the upper part of the diagram situation is quite different. Thermal fluctuations, not present in MFA, destroy ordered linear chains structures at temperatures much smaller than MFA critical temperature (the MFA critical temperature for the ($`2\times 1`$) phase can be even two times greater!). Moreover, MFA results do not allow to reach long periodic phases $`(4\times 1)`$ and $`(6\times 1)`$ directly from the disordered phase. On the other hand results of MC simulations show that all ordered phases can be reached from disordered phase in wide range of chemical potential. It is interesting to look at results of MC study presented in the (temperature - coverage) phase diagram (see Fig. 6).

Ordered phases of ($`p\times 1)`$ are separated by mixed phases $`(p,p^{})`$ which denotes mixture of ($`p\times 1)`$ phase and ($`p^{}\times 1)`$ phase. It is worth noting that ordered chain structures with period 6 and 9, present in the phase diagram, were not observed in experimentsnaumov73  which were performed at temperatures above 77K. We expect that these phases could be observed experimentally at temperatures below 77K.

Now, we are going to discuss the order of the phase transitions. The phase transitions between the ordered phases, e.g., $`(2\times 1)`$ $``$ $`(3\times 1)`$, are of the first order because sublattice order parameters and structure factors change discontinuously at the transition points. The analysis of the phase transition from an ordered phase to the disordered phase is much more complicated. Checking temperature dependence of the structure factor and heat capacity (some results are shown in Fig. 7)

we found that the transition to $`(9\times 1)`$ and $`(6\times 1)`$ phases is the first-order phase transition because of discontinuity of the structure factor. The transition to the $`(3\times 1)`$ and $`(4\times 1)`$ phases seems to be continuous transition. It is difficult to state the order of the phase transition form the disordered phase to the $`(2\times 1)`$ phase. Although, the structure factor seems to be continuous function of temperature, the heat capacity is smeared near the transition temperature indicating the first-order phase transition. On the oder hand, in the MFA calculation this phase transition is of the second order. One of possible explanation of this discrepancy is based on observation of additional, very weak maxima of the structure factor with the wave vector close to $`q=(0,\pi /(a_y))`$ and temperature in the narrow range of the transition temperature. This can suggest that defected phases like the 2(a) phase are still present close to the transition temperature. In order to answer whether this problem is connected with the finite range of indirect interaction one should performed finite-size calculation with infinite range of the interaction. It can be done by employing a numerical method to accelerate the convergence of a Fourier series oly96 .

## V Discussion and conclusions

Our Monte Carlo analysis revealed that truncation of the indirect interaction changes a phase diagram even if the interaction range is very large. New phases containing domain walls occur at high temperatures between the disordered phase and the $`(p\times 1)`$ phase. This problem complicates study of critical behaviour via finite-size scaling. One possible solution of such problem is accounting for infinite range of indirect interaction. It can be carried out by employing a numerical method to accelerate the convergence of a Fourier seriesoly96  in the same way as in earlier paper LO . However, it would be time consuming way to study the critical behaviour by Monte Carlo simulations. On the other hand it would be interesting to extend an analysis of defected phases (e.g., ($`2a`$) phase) generated by truncation of indirect interaction. The structure factor of the ($`2a`$) phase has maximum at the wavevector slightly shifted from the position which corresponds to maximum in the ordered $`(2\times 1)`$ phase. Hence it might be important in analysis of LEED patterns.

Present MC calculations demonstrate that results of the mean field approximation are very bad especially at high temperatures. We have shown that MFA transition temperature from the disordered phase to the $`(2\times 1)`$ phase is almost two times greater. Moreover, the high temperature part of the phase diagram has different topology than that obtained by MC simulations. However, it is interesting to note that both method give very similar results at low temperatures. Present MC calculations confirm earlier obtained results BLO  that at low coverages of lithium on the W(112) the long periodic structures ($`9\times 1)`$ and ($`6\times 1)`$ are stable at temperatures below 77K. It would be very interesting to determine experimentally the whole (T,$`\mathrm{\Theta }`$) phase diagram to check theoretical prediction. We are aware that present lattice gas model contains some phenomenological parameters, e.g., short range interaction $`E_b`$ which has influence on the transition temperature. We think that the better way to estimate the model parameters is to perform the first-principles calculations of the interaction potentials in the Li/W(112) system.

Recently Yakovkin studied yakov93  chain structures adsorbed on Mo(112) and Re(1010) using lattice gas model in which the indirect interaction takes nonzero value only for one from all lattice distances. Performing MC simulations in canonical ensemble he was able to reconstruct the sequence of phases ($`4\times 1`$) and ($`2\times 1`$) as observed experimentally in Li/Mo(112). It is interesting whether such simple model could describe correctly the phases observed in adsorption of lithium on W(112).

We observed that Monte Carlo study of phase diagram in canonical ensemble is very difficult due to dependence of final MC results on initial lattice gas configuration. Although similar difficulty appears in grand canonical ensemble it is rather easy to find ground-state configurations and calculate free energies for different phases. However, it would be interesting to look for different algorithm than single spin flip in order to overcome the dependence MC results on initial configurations.


format: text

text:
# 1 Introduction

## 1 Introduction

A great challenge in understanding the hard interactions has been posed by the HERA small-$`x`$ data  in order to describe the strong growth of the inclusive structure function $`F_2`$ as the Bjorken scale $`x`$ decreases. This feature is also supplemented by the scaling violations on the hard scale given by the photon virtuality $`Q^2`$. The currently high accuracy reached turns out the measurements having quite small statistic uncertainties, probing a large interval of virtualities, i.e. from a real photon to thousands of GeVs. As far the Regge approach is concerned, at high energy the $`ep`$ scattering process is dominated by the exchange of the Pomeron trajectory in the $`t`$-channel. From the hadronic phenomenology, this implies that the structure function would present a mild increasing on energy ($`sQ^2/x`$) since the soft Pomeron intercept ranges around $`\alpha _{IP}(0)1.08`$. Such behavior is in contrast with high energy $`ep`$ data, where the effective intercept takes values $`\lambda _{\mathrm{eff}}1.31.4`$. In the Regge language this situation can be solved by introducing the idea of new poles in the complex angular momentum plane, for instance rendered in the multipoles models , producing quite successful data description. Other proposition is given by the two Pomeron model , introducing an additional hard intercept and corresponding residue. However, a shortcoming from these approaches is the poor knowledge about the behavior on virtuality, in general modeled in an empirical way through the vertex functions. An interesting mixed procedure is to fix the initial conditions for the QCD parton distributions (valence and sea/gluon) in a sufficiently large initial scale $`Q_0^215`$ GeV2 from the Regge phenomenology and then perform the QCD evolution up to higher virtualities .

On the other hand, the high photon virtuality allows the applicability of the QCD perturbative methods. Two main approaches emerge from the QCD formalism: the DGLAP and BFKL evolution schemes. The DGLAP formalism  is quite successful in describing most of the measurements on structure functions at HERA and hard processes in the hadronic colliders. This feature is even intriguing, since its theoretical limitations at high energy are well known . Other perturbative approach is the BFKL formalism , well established at LO level but not yet completely understood at NLO accuracy. The main issue in the NLO BFKL effects is the correct account of the sub-leading corrections in the all orders resummation  (for a pedagogical review, see for instance ). The LO BFKL approach can describe HERA structure function in a limited kinematical range, i.e. at not so large $`Q^2`$ and small-$`x`$. The main assumption at LO is that the processes described through the so-called hard Pomeron are driven by gluon ladder diagrams with infinite rungs of $`s`$-channel gluons at asymptotic energies. A consistent treatment considering higher order resummations is currently in progress and should be available soon. An open issue is whether the available energies in the nowadays colliders actually are asymptotic allowing the use of the complete BFKL series.

In this work we perform a global data fitting to the HERA structure function $`F_2`$ using as a model for the hard Pomeron the finite sum of gluon ladders . The model is based on the truncation of the BFKL series keeping only the first few orders in the strong coupling perturbative expansion, where subleading contributions can be absorbed in the adjustable parameters. From the phenomenology on hadronic collisions , just three orders, $`[\alpha _s\mathrm{ln}(1/x)]^2`$, are enough to describe current accelerator data. The hard Pomeron model should be supplemented by a soft piece accounting for the non-perturbative contributions to the process. The description therefore turns out similar to the two Pomeron model , with the advantage of a complete knowledge of the behavior on $`x`$ and $`Q^2`$. Concerning the energy dependence, the hard piece has the logarithmic growth in contrast with the effective powerlike behavior on the model . The original model contains a reduced number of adjustable parameters: the normalization $`𝒩_p`$ and the non-perturbative scale $`\mu ^2`$ from the proton impact factor and the parameter $`x_0`$ scaling the logarithms on energy. More details on those parameter are addressed in the next section. In a previous study at Ref. , two distinct choices for the soft Pomeron were analised. The resulting two Pomeron model was successful in describing data on structure function $`F_2`$ and its derivatives (slopes on $`Q^2`$ and $`x`$) for $`x2.510^2`$ and $`0.045Q^21500`$ GeV2. Here, we extend the kinematical range of fitting and the models considered for the soft Pomeron.

This paper is organized as follows. In the next section, we shortly review the main expressions for the hard piece given by the summation up to the two-rung ladder contribution. In Sec. (3), an overall fit to the recent deep inelastic data is performed based on the hard contribution referred above supplemented by the remaining soft Pomeron and non-singlet contributions. The role played by the hard and soft contributions are investigated in details. There, we also draw up our conclusions.

## 2 The hard contribution: summing gluon ladders

In this section, we review the elements needed to compute the structure function using the finite sum of gluon ladders in the $`ep`$ collision, with the photon-proton center of mass energy labeled by $`W`$. Before that, let us shortly motivate the picture of the finite sum of gluon ladders. At finite energies, the LLA and NLLA summation implies that the amplitude is represented by a finite sum on terms, where the number of terms increases like $`\mathrm{ln}s`$, rather than by the solution of the BFKL integral equation. The interest in to take the firts terms in the completete series in the truncation is related to the fact that the energies reached by the present accelators are not high enough to accommodate a big number of gluons in the ladder rungs that eventually hadronize. In the energy range of HERA there can be only a few real gluons produced in any scattering event. On the other hand, in the LO BFKL resummation real gluons can be emitted without any cost in energy, while in reality the production of a real gluon requires an energy equivalent to one or one and half units of rapidity. This violation of energy conservation is probably cured in the NLO BFKL resummation or by relying on consistency constraint implementing energy conservation in numerical simulations of LO BFKL evolution. The truncation of the whole series could be similar to this consistency condition. Corrobating the truncation hypothesis, for example the coefficient weighting the term $`\mathrm{ln}^3s`$ turn out to be compatible with zero considering even the Tevatron data , in contrast with the expected from the complete resummation.

The proton inclusive structure function, written in terms of the cross sections for the scattering of transverse or longitudinal polarized photons, reads as

$`F_2(x,Q^2)`$ $`=`$ $`{\displaystyle \frac{Q^2}{4\pi ^2\alpha _{\mathrm{em}}}}\left[\sigma _T(x,Q^2)+\sigma _L(x,Q^2)\right],`$ (1)
$`\sigma _{T,L}(x,Q^2)`$ $`=`$ $`{\displaystyle \frac{𝒢}{(2\pi )^4}}{\displaystyle \frac{d^2𝒌}{𝒌^2}\frac{d^2𝒌^{\mathbf{}}}{𝒌_{}^{\mathbf{}}{}_{}{}^{2}}\mathrm{\Phi }_{T,L}^\gamma ^{}(𝒌)𝒦(x,𝒌,𝒌^{\mathbf{}})\mathrm{\Phi }_p(𝒌^{\mathbf{}})},`$ (2)
where $`𝒢`$ is the color factor for the color singlet exchange and $`𝒌`$ and $`𝒌^{\mathbf{}}`$ are the transverse momenta of the exchanged reggeized gluons in the $`t`$-channel. The $`\mathrm{\Phi }_{T,L}^\gamma ^{}(𝒌)`$ is the virtual photon impact factor and $`\mathrm{\Phi }_p(𝒌^{\mathbf{}})`$ is the proton impact factor. The first one is well known in perturbation theory at leading order, while the latter is modeled since in the proton vertex there is no hard scale to allow pQCD calculations. The kernel $`𝒦(x,𝒌,𝒌^{\mathbf{}})`$ contains the dynamics of the process, for instance, the BFKL kernel.

The amplitudes can be calculated order by order: for instance the Born contribution coming from the two gluon exchange and the one-rung ladder contribution read as,

$`𝒜^{(0)}`$ $`=`$ $`{\displaystyle \frac{2\alpha _sW^2}{\pi ^2}}{\displaystyle \underset{f}{}}e_f^2{\displaystyle \frac{d^2𝒌}{𝒌^4}\mathrm{\Phi }_{T,L}^\gamma ^{}(𝒌)\mathrm{\Phi }_p(𝒌)},`$
$`𝒜^{(1)}`$ $`=`$ $`{\displaystyle \frac{6\alpha _s^2W^2}{8\pi ^4}}{\displaystyle \underset{f}{}}e_f^2\mathrm{ln}\left({\displaystyle \frac{W^2}{W_0^2}}\right){\displaystyle \frac{d^2𝒌}{𝒌^4}\frac{d^2𝒌^{\mathbf{}}}{𝒌_{}^{\mathbf{}}{}_{}{}^{4}}\mathrm{\Phi }_{T,L}^\gamma ^{}(𝒌)𝒦(𝒌,𝒌^{\mathbf{}})\mathrm{\Phi }_p(𝒌^{\mathbf{}})},`$
where $`\alpha _s`$ is considered fixed since we are in the framework of the LO BFKL approach. The perturbative kernel $`𝒦(𝒌,𝒌^{\mathbf{}})`$ can be calculated order by order in the perturbative expansion . The Pomeron is attached to the off-shell incoming photon through the quark loop diagrams, where the Reggeized gluons are attached to the same and to different quarks in the loop. The virtual photon impact factor averaged over the transverse polarizations reads as ,

$`\mathrm{\Phi }_{T,L}^\gamma ^{}(𝒌)={\displaystyle \frac{1}{2}}{\displaystyle _0^1}{\displaystyle \frac{d\tau }{2\pi }}{\displaystyle _0^1}{\displaystyle \frac{d\rho }{2\pi }}{\displaystyle \frac{𝒌^2(12\tau \tau ^{})(12\rho \rho ^{})}{𝒌^2\rho \rho ^{}+Q^2\rho \tau \tau ^{}}},`$ (3)
where $`\rho `$, $`\tau `$ are the Sudakov variables associated with the momenta in the photon vertex and $`\tau ^{}(1\tau )`$ and $`\rho ^{}(1\rho )`$.

Gauge invariance requires the proton impact factor vanishing at $`𝒌^{\mathbf{}}`$ going to zero and is modeled in a simple way,

$`\mathrm{\Phi }_p(𝒌^{\mathbf{}})=𝒩_p{\displaystyle \frac{𝒌_{}^{\mathbf{}}{}_{}{}^{2}}{𝒌_{}^{\mathbf{}}{}_{}{}^{2}+\mu ^2}},`$ (4)
where $`𝒩_p`$ is the unknown normalization of the proton impact factor and $`\mu ^2`$ is a scale which is typical of the non-perturbative dynamics. These scales will be considered adjustable parameters in the analysis. Considering the electroproduction process, summing the first orders in perturbation theory we can write the expression for the inclusive structure function,

$`F_2^{\mathrm{hard}}(x,Q^2)`$ $`=`$ $`{\displaystyle \frac{8}{3}}{\displaystyle \frac{\alpha _s^2}{\pi ^2}}{\displaystyle \underset{f}{}}e_f^2𝒩_p[I^{(0)}(Q^2,\mu ^2)+{\displaystyle \frac{3\alpha _s}{\pi }}\mathrm{ln}{\displaystyle \frac{x_0}{x}}I^{(1)}(Q^2,\mu ^2)`$ (5)
$`+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{3\alpha _s}{\pi }}\mathrm{ln}{\displaystyle \frac{x_0}{x}}\right)^2I^{(2)}(Q^2,\mu ^2)],`$
where the functions $`I^{(n)}(Q^2,\mu ^2)`$ correspond to the $`n`$-rung gluon ladder contribution. The quantity $`x_0`$ gives the scale normalizing the logarithms on energy for the LLA BFKL approach, which is arbitrary and enters as an additional parameter. The contributions are written explicitly as,

$`I^{(0)}`$ $`=`$ $`{\displaystyle \frac{1}{2}}\mathrm{ln}^2\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{7}{6}}\mathrm{ln}\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{77}{18}},`$ (6)
$`I^{(1)}`$ $`=`$ $`{\displaystyle \frac{1}{6}}\mathrm{ln}^3\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{7}{12}}\mathrm{ln}^2\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{77}{18}}\mathrm{ln}\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{131}{27}}+2\zeta (3),`$ (7)
$`I^{(2)}`$ $`=`$ $`{\displaystyle \frac{1}{24}}\mathrm{ln}^4\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{7}{36}}\mathrm{ln}^3\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{77}{36}}\mathrm{ln}^2\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)`$ (8)
$`+\left({\displaystyle \frac{131}{27}}+4\zeta (3)\right)\mathrm{ln}\left({\displaystyle \frac{Q^2}{\mu ^2}}\right)+{\displaystyle \frac{1396}{81}}{\displaystyle \frac{\pi ^4}{15}}+{\displaystyle \frac{14}{3}}\zeta (3),`$
where $`\zeta (3)1.202`$. It is worth mentioning the clear behavior on $`x`$ and virtuality. The main result in Ref.  is in a good agreement, in terms of a $`\chi ^2/\mathrm{dof}`$ test, for the inclusive structure function in the range $`0.045Q^21500`$ GeV2, once considering the restricted kinematical constraint $`x0.025`$. The non-perturbative contribution (from the soft dynamics), initially considered as a background, was found to be not negligible. In the next section, we perform a global analysis, extending the range on $`x`$ fitted by adding the non-singlet contribution modeled through the usual Regge parameterizations and analyzing also the multipole models for the soft Pomeron.

## 3 Fitting results and conclusions

In order to perform the fitting procedure, for the hard piece one uses Eq. (5) and for the soft piece we have selected some typical models, as addressed below. First, one takes into account the latest version  of the CKMT model , having a economical number of parameters,

$`F_2^{\mathrm{soft}}(x,Q^2)=A\left({\displaystyle \frac{1}{x}}\right)^{\mathrm{\Delta }(Q^2)}\left({\displaystyle \frac{Q^2}{Q^2+a}}\right)^{\mathrm{\Delta }(Q^2)+1}\left(1x\right)^{n_s(Q^2)},`$ (9)
$`\mathrm{\Delta }(Q^2)=\mathrm{\Delta }_0\left(1+{\displaystyle \frac{Q^2\mathrm{\Delta }_1}{\mathrm{\Delta }_2+Q^2}}\right),n_s(Q^2)={\displaystyle \frac{7}{2}}\left(1+{\displaystyle \frac{Q^2}{Q^2+b}}\right),`$ (10)
where $`\mathrm{\Delta }(Q^2)`$ is the Pomeron intercept. The non-singlet term takes the following form,

$`F^{\mathrm{ns}}(x,Q^2)`$ $`=`$ $`A_Rx^{1\alpha _R}\left({\displaystyle \frac{Q^2}{Q^2+a_R}}\right)^{\alpha _R}(1x)^{n_{ns}(Q^2)}`$ (11)
$`n_{ns}(Q^2)`$ $`=`$ $`{\displaystyle \frac{3}{2}}\left(1+{\displaystyle \frac{Q^2}{Q^2+d}}\right).`$ (12)
The CKMT model was originally constructed to interpolate between the soft hadronic Pomeron phenomenology, where $`\alpha _{pom}1.08`$ and the growth for the low $`Q^2`$ proton structure function in deep inelastic scattering, where $`\alpha _{pom}1.2`$. This is obtained through absorbtive corrections to the bare soft Pomeron, leading to a $`Q^2`$-dependent intercept. Therefore, in phenomenological application considering a two-Pomeron approach, we should be careful in verifying the soft Pomeron intercept coming out.

Another possibility is to select a soft piece which corresponds to the Pomeron with intercept equal to one and has the form of non-perturbative truncated $`\mathrm{ln}\left(\frac{Q^2}{x}\right)`$ series (soft multipole Pomeron) , given by

$$F_2^{\mathrm{soft}}(x,Q^2)=Q^2\left[A\left(\frac{a}{Q^2+a}\right)^\alpha +B\mathrm{ln}\left(\frac{Q^2}{x}\right)\left(\frac{a}{Q^2+a}\right)^\beta \right]$$
(13)
for the dipole Pomeron and

$$F_2^{\mathrm{soft}}(x,Q^2)=Q^2\left[A\left(\frac{a}{Q^2+a}\right)^\alpha +B\mathrm{ln}^2\left(\frac{Q^2}{x}\right)\left(\frac{a}{Q^2+a}\right)^\beta \right]$$
(14)
for the tripole Pomeron.

The multipole means the soft Pomeron is considered as multiple poles instead of just a single pole. This approach has been successful in hadron-hadron and lepton-hadron (mostly vector meson photoproduction) scattering, where the increasing of the cross section on energy can be produced with a unit Pomeron intercept. The multipole model has shown to be quite stable in fitting simultaneously hadronic and lepton-hadron data.

Concerning the hard piece, Eq. (5), we selected the overall normalization factor as a free parameter, defined as $`𝒩=\frac{8}{3}\frac{\alpha _s^2}{\pi ^2}e_f^2𝒩_p`$, considering four active flavors. That contribution was also supplemented by a factor $`(1x)^{n_s}`$ accounting for the large $`x`$ effects. The remaining parameters are the scales $`\mu ^2`$ and $`x_0`$ as well as the coupling constant $`\alpha _s`$ (considered at a fixed scale in the BFKL formalism). As we will see in the following, its value presents small variation in the fitting, suggesting that it can be considered fixed as $`\alpha _s=0.2`$, which is typical in the HERA kinematical region.

For the fitting procedure we consider the data set containing all available HERA data for the proton structure function $`F_2`$ -, updated with new HERA experiment  and adding only the latest (E665 and NMC) data set of fixed target experiments . This choice for the selection of data sets follows the procedures used in Ref. . The whole data set contains a total number of $`1059`$ experimental points, covering the complete available $`x`$ and virtualities $`0.045Q^230000`$ GeV2. It should be stressed that the inclusion of the older fixed target experimental data (SLAC  and BCDMS ) requires a more sophisticated study for the background, for example, following the ALLM parameterization or as in Ref. . However, for that purposes the number of free parameters would be significantly larger.

Here, we have considered the following distinct fitting procedures:

* (I) Overall fit considering the hard piece and the CKMT soft Pomeron plus reggeon piece, Eqs. (9) and (11);
* (II) Overall fit restricting the soft Pomeron intercept, modifying the expression for $`\mathrm{\Delta }(Q^2)`$, see Eq. (10), in the CKMT soft Pomeron. It reads now as,

$`\mathrm{\Delta }(Q^2)=\mathrm{\Delta }_0\left({\displaystyle \frac{Q^2}{\mathrm{\Delta }_1+Q^2}}\right)+\mathrm{\Delta }_2,`$ (15)
* (III) The fit is performed considering only the hard piece plus the non-singlet contribution;
* (IV and V) Overall fit using the hard Pomeron piece and the soft multipoles Pomeron, Eqs. (13) and (14). The main feature is the soft pomeron having an intercept equal one.

The best fit parameters for these procedures are presented in Table (1). The kinematical range covered is $`0x1`$ and $`0.045Q^230000`$. In the following we discuss each procedure, pointing out the relative contribution from the hard and soft pieces and its quality.

Lets start by the procedure (I). The parameters for the hard piece remains conssistent with the previous analysis in , i.e. $`\alpha _s0.2`$, $`x_010^2`$ and $`\mu ^2`$ of order 1 GeV1. Concerning the soft Pomeron, the bare intercept $`\mathrm{\Delta }_0`$ coming out quite small in the whole range of $`Q^2`$. This fact seems corroborate a soft Pomeron having an intercept close to one in this two-Pomeron analysis. The same happens for the procedure (II), with a modified form for the CKMT intercept. For low $`Q^2`$ the intercept comes out small, but reaches at higher virtualities to the upper limit $`\mathrm{\Delta }=0.08`$.

The procedures (IV and V), considering the soft multipole Pomeron, provide the best quality of fit $`\chi ^2/\mathrm{dof}1`$. The tripole Pomeron is preferred as a soft background, giving a slightly better result than the dipole pomeron. Concerning the hard piece, the value for the fixed coupling constant is quite stable, $`\alpha _s0.2`$, consistent with the overall HERA value.

The situation differs only in procedure (III), where a small value for $`\alpha _s`$ is found, suggesting in this case that the best choice would be $`\alpha _s=0.119.`$ which coincides with $`\alpha _s(M_Z)`$. The high $`\chi ^2/\mathrm{dof}`$ for this case is not quite sizeable, in view of the smaller number of parameters considered (8 against 13 from the remaining analysis).

Let start by the procedure (I). The parameters for the hard piece remains consistent with the previous analysis in , i.e. $`\alpha _s0.2`$, $`x_010^2`$ and $`\mu ^2`$ of order 1 GeV2. Concerning the soft Pomeron, the bare intercept $`\mathrm{\Delta }_0`$ comes out quite small in the whole range of $`Q^2`$. This fact seems to corroborate a soft Pomeron having an intercept close to one in this two-Pomeron analysis. The same happens for the procedure (II), with a modified form for the CKMT intercept. For low $`Q^2`$ the intercept comes out small, but reaches at higher virtualities to the upper limit $`\mathrm{\Delta }=0.08`$.

The procedures (IV and V), considering the soft multipole Pomeron, provide the best quality of fit $`\chi ^2/\mathrm{dof}1`$. The tripole Pomeron is preferred as a soft background, giving a slightly better result than the dipole pomeron. Concerning the hard piece, the value for the fixed coupling constant is quite stable, $`\alpha _s0.2`$, consistent with the overall HERA value.

The situation differs only in procedure (III), where a small value for $`\alpha _s`$ is found, suggesting in this case that the best choice would be $`\alpha _s=0.119,`$ which coincides with $`\alpha _s(M_Z)`$. The high $`\chi ^2/\mathrm{dof}`$ for this case is not quite sizeable, in view of the smaller number of parameters considered (8 against 13 from the remaining analysis).

In Figs. 1 and 2, we present the fit result for the inclusive structure function for small and large virtualities. The plots for the different procedures lie on top of each other. In Fig. 3, we present the relative role of the distinct contributions for the fit. The hard and soft pieces are presented separately. We present them explicitly for the virtualities 15, 150 and 1500 GeV2, where the contributions from fitting (I), (II) and (V) are shown. It is verified that the region where the hard Pomeron starts to dominate depends on the vituality. For instance, at $`Q^2=15`$ GeV2 it stays on $`x10^3`$ for the different procedures, whereas is shifted to $`x10^2`$ at $`Q^2=1500`$ GeV2. A two-Pomeron picture is supported, comparable with those using the two-Pomeron analysis in Refs. .

Additionally, we performed the fit restricting the domain of the Bjorken variable to $`x0.07`$. This procedure is similar to that considered in analysis of Ref.. For this purpose we restricted the number of free parameters (equal to 10) for Soft+Hard Pomeron model, as in the Donnachie-Landshoff model reanalysis considered in . To this aim, we fixed the constant $`\alpha _s=0.2`$ as well as the powers $`n_s=7`$ and $`n_{ns}=3`$. The following procedures were considered: (A): hard plus soft terms; (B): hard + soft dipole pomeron terms; (C): hard + soft tripole pomeron terms. Results of the fit are presented in Table 2. The quality of our approach is similar those ones found in Ref. , at least in this specific domain, $`x0.07`$. The parameters for the multipole Pomeron seem stable in this kinematical region. There is a change in the case (A), where the bare Pomeron intercept is higher than in the procedure (I) producing an usual soft result $`\mathrm{\Delta }0.08`$. The parameters for the hard piece remain stable in all fitting procedures.

In conclusion, we verify that the fitting procedure is equivalent to the model using the two-Pomeron approach , with the advantage of clear understanding of the behaviors on $`x`$ and $`Q^2`$ of the corresponding hard content. In particular, the behavior on $`Q^2`$ of the hard piece sheds light about on further implementations of the hard Pomeron residue in two-Pomeron fitting. The fitting procedure shows that the soft Pomeron has an intercept small than the usual $`\alpha _{IP}=1.08`$, suggesting that a suitable choice is the multipole pomeron having intercept equal to one. In particular, the tripole pomeron presented the best fit result in all kinematical ranges considered here, followed by the dipole one. The fitting using only the hard piece and the non-singlet contributions is not completely ruled out, producing a not so high $`\chi ^2/\mathrm{dof}`$ result.

## Acknowledgments

M.V.T.M. thanks the support of the High Energy Physics Phenomenology Group (GFPAE, IF-UFRGS) at Institute of Physics, Porto Alegre.


text:
# Mapping the Wigner distribution function of the Morse oscillator into a semi-classical distribution function

## I Introduction

The Wigner description in phase-space  provides a tool for comparing quantum and classical dynamics . In a previous paper  we introduced a mapping relating the Wigner distribution function (WDF) corresponding to a given wave function, solution of the Schrödinger equation to a semiclassical distribution function (SDF) satisfying the classical Liouville equation with the same potential. So far this mapping was applied to the ground state of the square well , the Pöschl-Teller potential  and to a modified harmonic oscillator potential .

Here this mapping is extended to the ground state of the unidimensional Morse oscillator. We vary the depth of the potential in order to study the effect of the level density on the mapping.

The Morse oscillator  plays an important role in many areas of Physics, because it is a good approximation to diatomic molecular potentials . Today the Morse potential is used in molecular spectroscopy, even for polyatomic molecules . Some authors have used the Morse potential for studying the semiclassical limit of Quantum Mechanics .

The unidimensional Morse oscillator in phase-space has been already studied by many authors. The classical motion  as well as the quantum representation  are well known for this potential.

In Sec. II of this work we summarize the method developed in Ref.  for mapping WDF into SDF. In Sec. III we study the Morse oscillator and in Sec. IV we present our results and conclusions.

## II Description of the method

In this section we present a short derivation of the method used for the mapping . We start with the quantum Liouville equation for the WDF $`\rho `$:

$$\frac{\rho }{t}+\frac{p}{m}\frac{\rho }{q}+K(q,pp^{})\rho (q,p^{},t)𝑑p^{}=0,$$
(1)
where the kernel $`K`$ is given by

$$K(q,pp^{})=\frac{i}{\mathrm{}}\frac{dv}{2\pi \mathrm{}}\mathrm{exp}\left[\frac{i}{\mathrm{}}(pp^{})v\right]\left[V\left(q\frac{v}{2}\right)V\left(q+\frac{v}{2}\right)\right],$$
(2)
$`V(q)`$ being the potential. The corresponding classical Liouville equation may be written similarly

$$\frac{\rho _c}{t}+\frac{p}{m}\frac{\rho _c}{q}+K_c(q,pp^{})\rho _c(q,p^{},t)𝑑p^{}=0,$$
(3)
where

$$K_c(q,pp^{})=\frac{i}{\mathrm{}}\frac{dv}{2\pi \mathrm{}}\mathrm{exp}\left[\frac{i}{\mathrm{}}(pp^{})v\right]v\frac{V}{q}=\frac{V}{q}\frac{}{p}\delta (pp^{}).$$
(4)
We relate solutions $`\rho `$ and $`\rho _c`$ of Eqs. (1) and (3) through the integral equation

$`\rho (q,p,t)`$ $`=`$ $`\rho _c(q,p,t){\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑t^{}{\displaystyle 𝑑q^{}𝑑p^{}G_c(q,p;q^{},p^{},tt^{})}`$ (5)
$`\times `$ $`{\displaystyle 𝑑p^{\prime \prime }[K(q^{},p^{}p^{\prime \prime })K_c(q^{},p^{}p^{\prime \prime })]\rho (q^{},p^{\prime \prime },t^{})},`$ (6)
where $`G_c`$ is the retarded Green’s function corresponding to the classical Liouville equation (3). In fact Eq. (6) is the defining equation of $`\rho _c`$ since the WDF $`\rho `$ is already determined by fixing the wave function solution of the Schrödinger equation. Using the differential equation satisfied by $`G_c`$  one may easily verify that $`\rho _c`$ satisfies Eq. (3) provided $`\rho `$ obeys Eq. (1).

As $`K_c`$ is the term of lowest order in the expansion of $`K`$ in a power series of $`\mathrm{}`$, $`\rho \rho _c`$ contains only corrections of the first and higher order terms in $`\mathrm{}`$.

It can be shown from results given in Ref. , for a WDF generated from solutions of the Schrödinger equation for which the initial wave function at the time $`t=t_0`$ is given; assuming $`\rho (q,p,t)`$ to vanish for $`t<t_0`$, that $`\rho _c`$ from Eq. (6) gives simply the classical evolution of $`\rho (q,p,t_0)`$, that is

$$\rho _c(q,p,t)=\rho (Q(q,p,t_0t),P(q,p,t_0t),t_0),tt_0.$$
(7)
Here $`(Q(q,p,\tau ),P(q,p,\tau ))`$ describes the classical trajectory in phase-space of a particle subject to the potential $`V(q)`$ which at the time $`\tau =0`$ occupies the point $`(q,p)`$, that is

$$q=Q(q,p,0),p=P(q,p,0).$$
(8)
In what follows we shall consider mostly the mapping of the time independent WDF $`\rho (q,p)`$ corresponding to a bound state. In this case, from Eq. (6) one obtains  the following prescription for $`\rho _c`$,

$$\rho _c(q,p)=\underset{ϵ0_+}{lim}ϵ_{\mathrm{}}^o𝑑\tau e^{ϵ\tau }\rho (Q(q,p,\tau ),P(q,p,\tau )),$$
(9)
which is also time independent. Here $`e^{ϵ\tau }`$ is a convergence factor introduced explicitly in the Green’s function.

It can be shown  that for points $`(q,p)`$ on closed classical orbits the expression given by Eq. (9) is equivalent to

$$\rho _c(q,p)=\frac{1}{T(q,p)}_{T(q,p)}^0𝑑t\rho (Q(q,p,t),P(q,p,t)),$$
(10)
where $`T=T(q,p)`$ is the period of the orbit. For open trajectories we have $`T\mathrm{}`$ and as for bound states the integral in Eq. (10) converges, we get $`\rho _c=0`$. As it will be discussed further, in the calculation of averages of physical quantities the quantity of interest is rather $`T(q,p)\rho _c(q,p)`$ which does not vanish for the open trajectories.

We observe here that one may verify directly that the time-independent Liouville equation is obeyed by applying the operator $`p/q(V/q)/p`$ on the right hand side of Eq. (9) or (10) and using the fact that $`\rho (Q(q,p,t),P(q,p,t))`$ satisfies the classical Liouville equation (3).

An important property of the stationary SDF is that it is constant along the classical trajectories  which can be seen by noticing that for any two points $`(q,p)`$ and $`(q^{},p^{})`$ on the same path in phase space one has a time interval $`\mathrm{\Delta }`$ such that

$$Q(q^{},p^{},t)=Q(q,p,t+\mathrm{\Delta }),P(q^{},p^{},t)=P(q,p,t+\mathrm{\Delta }).$$
(11)
($`\mathrm{\Delta }`$ is the time which the particle takes to go from the point $`(q,p)`$ to the point $`(q^{},p^{})`$ in phase space or vice versa). Thus we get

$`\rho _c(q^{},p^{})`$ $`=`$ $`{\displaystyle \frac{1}{T}}{\displaystyle _T^0}𝑑t\rho (Q(q^{},p^{},t),P(q^{},p^{},t))`$ (12)
$`=`$ $`{\displaystyle \frac{1}{T}}{\displaystyle _T^0}𝑑t\rho (Q(q,p,t+\mathrm{\Delta }),P(q,p,t+\mathrm{\Delta }))`$ (13)
$`=`$ $`\rho _c(q,p),`$ (14)
where we used the periodicity of the integrand in the last step.

Let $`(Q(E,\nu ,t),P(E,\nu ,t))`$ represent, for $`t`$ in the interval $`T(E,\nu )/2tT(E,\nu )/2`$ ($`T`$ may be infinity), the points of a trajectory in phase-space corresponding to energy $`E`$ and fixed value of $`\nu `$. Here the discrete parameter $`\nu `$ ($`\nu =1,\mathrm{},\nu _{\mathrm{max}}`$) distinguishes between the different trajectories corresponding to the same value of $`E`$. By varying $`E`$ and $`\nu `$ one covers the entire allowed phase space. Thus we may consider the transformation of points ($`E,\nu ,t`$) in the appropriate domains $`D_\nu `$ of the energy-time space onto the points ($`p,q`$) of the phase space

$$q=Q(E,\nu ,t),p=P(E,\nu ,t).$$
(15)
Let us define the function

$$R_c(E,\nu )=_{T(E,\nu )/2}^{T(E,\nu )/2}𝑑t\rho (Q(E,\nu ,t),P(E,\nu ,t)).$$
(16)
According to Eq. (10) one has

$$\rho _c(q,p)=\frac{R_c(E(q,p),\nu )}{T(E(q,p),\nu )},$$
(17)
by choosing the value of $`\nu `$ appropriate to the trajectory containing the point ($`q,p`$) as follows from Eq. (15). We observe here that the relationship (16) between $`\rho (q,p)`$ and $`R_c(E,\nu )`$ is analogous to that between $`R(q)`$, the probability density in coordinate space and the WDF $`\rho (q,p)`$ :

$$R(q)=\rho (q,p)𝑑p.$$
(18)
The average of any Weyl function  $`𝒪(q,p)`$ corresponding to a certain operator $`O`$ may be written

$$O=𝒪(q,p)\rho (q,p)𝑑p𝑑q=\underset{\nu }{}_{𝒟_\nu }o(E,\nu ,\tau )r(E,\nu ,\tau )𝑑E𝑑\tau ,$$
(19)
where we introduced the functions

$`o(E,\nu ,\tau )=𝒪(Q(E,\nu ,\tau ),P(E,\nu ,\tau )),`$ (20)
$`r(E,\nu ,\tau )=\rho (Q(E,\nu ,\tau ),P(E,\nu ,\tau )),`$ (21)
and we used the fact that the Jacobian $`J`$ of the transformation (15) is unity

$$J=\left|\begin{array}{cc}\frac{P}{E}& \frac{Q}{E}\\ \frac{P}{t}& \frac{Q}{t}\end{array}\right|=\frac{P}{E}\frac{E}{p}+\frac{Q}{E}\frac{E}{q}=1.$$
(22)
In the special case in which the Weyl function $`𝒪`$ is a constant of motion depending on $`(q,p)`$ only through the energy $`E(q,p)`$, we obtain from (19), (18) and (16)

$$O=\underset{\nu }{}_{D_\nu }o(E,\nu )r(E,\nu ,\tau )𝑑E𝑑\tau =\underset{\nu }{}o(E,\nu )R_c(E,\nu )𝑑E.$$
(23)
In particular the normalization condition for $`R_c(E,\nu )`$

$$\underset{\nu }{}R_c(E,\nu )𝑑E=1,$$
(24)
follows from Eq. (23) by taking $`𝒪(q,p)=1`$ since the Wigner function $`\rho (q,p)`$ is normalized.

As $`\rho _c`$ is an approximation correct in zeroth order of the expansion in powers of $`\mathrm{}`$ of the Wigner function $`\rho `$, the average

$$O_c=𝒪(q,p)\rho _c(q,p)𝑑q𝑑p,$$
(25)
is also correct in the same order in $`\mathrm{}`$. The expression (25) may be also written, by making use of the transformation (15), as

$`O_c={\displaystyle \underset{\nu }{}}{\displaystyle _{D_\nu }}o(E,\nu ,t){\displaystyle \frac{R_c(E,\nu )}{T(E,\nu )}}𝑑E𝑑t={\displaystyle \underset{\nu }{}}{\displaystyle \overline{o}(E,\nu )R_c(E,\nu )𝑑E},`$ (26)
where we used Eqs. (17) and (17) and $`\overline{o}(E,\nu )`$ denotes the average

$$\overline{o}(E,\nu )=_{T(E,\nu )/2}^{T(E,\nu )/2}o(E,\nu ,t)\frac{dt}{T(E,\nu )}.$$
(27)
Thus the average $`O_c`$ is obtained by replacing $`o(E,\nu ,t)`$ in Eq. (19) by the average (27). If, for fixed $`E`$, $`o(E,\nu ,t)`$ depends weakly on $`t`$, the average $`O_c`$ is expected to be a good approximation.

## III  Wigner distribution function for the Morse potential

In 1929 Morse  suggested the potential $`U_0(1e^{\alpha r})^2`$ for studying diatomic molecules. The Schrödinger equation for this potential does not have an exact solution, but for the one dimensional case the problem can be solved analytically .

In order to obtain the Wigner Distribution Function (WDF)

$$\rho (q,p,t)=(2\pi \mathrm{})^1_{\mathrm{}}^{\mathrm{}}𝑑ve^{\frac{i}{\mathrm{}}pv}\psi (q\frac{1}{2}v,t)\psi ^{}(q+\frac{1}{2}v,t),$$
(28)
where $`\psi `$ is the solution of the Schrödinger equation, we need the eigenfunctions for the one-dimensional Morse potential

$$𝒱(x)=D(1e^{ax})^2,$$
(29)
$`a`$ and $`D`$ being constant parameters. Starting with the Schrödinger equation

$$\left[\frac{\mathrm{}^2}{2m}\frac{d^2}{dx^2}+D(1e^{ax})^2\right]\psi _n=_n\psi _n,$$
(30)
we introduce the dimensionless parameter $`\lambda `$

$$\lambda =\frac{\sqrt{2mD}}{a\mathrm{}},$$
(31)
and the dimensionless coordinate

$$q=ax,$$
(32)
obtaining an eigenvalue equation depending only on one parameter

$$\left[\frac{1}{\lambda ^2}\frac{d^2}{dq^2}+(1e^q)^2\right]\psi _n=ϵ_n\psi _n,$$
(33)
where

$$ϵ_n=\frac{_n}{D}.$$
(34)
This equation is solved most conveniently using the variable

$$\xi =2\lambda e^q,\mathrm{}<q<\mathrm{}.$$
(35)
The eigenfunctions and eigenvalues are

$`\psi _n(\xi )=N(n,\lambda )e^{\frac{\xi }{2}}\xi ^{\lambda n\frac{1}{2}}L_n^{\lambda n\frac{1}{2}}(\xi ),`$ (36)
$`ϵ_n={\displaystyle \frac{2}{\lambda }}\left(n+{\displaystyle \frac{1}{2}}\right)\left[1{\displaystyle \frac{1}{2\lambda }}\left(n+{\displaystyle \frac{1}{2}}\right)\right],`$ (37)
where the quantum number $`n`$ takes the values

$$n=0,1,2,\mathrm{}[\lambda \frac{1}{2}].$$
(38)
Here $`[\lambda 1/2]`$ denotes the largest integer smaller than $`\lambda 1/2`$, $`L_n^s`$ is the polynomial

$$L_n^s=\underset{j=0}{\overset{n}{}}\left(\begin{array}{c}n+2s\\ nj\end{array}\right)\frac{(\xi )^j}{j!},$$
(39)
and the normalization factor $`N`$ is given by

$$N(\lambda ,n)=\left[\frac{(2\lambda 2n1)\mathrm{\Gamma }(n+1)}{\mathrm{\Gamma }(2\lambda n)}\right]^{\frac{1}{2}},$$
(40)
where the normalization condition $`\psi _n^{}\psi _n𝑑q=1`$ is assumed. In Appendix A, Eqs. (33) and (37) are derived.

From Eq. (37) one verifies that, for $`\lambda 1`$ and $`n\lambda `$, the energy spectrum of the Morse oscillator is written approximately $`_n\mathrm{}\omega _0(n+1/2)`$, which is the spectrum of a harmonic oscillator with frequency

$$\omega _0=\frac{2D}{\mathrm{}\lambda }=a\sqrt{\frac{2D}{m}}=\lambda \mathrm{}\frac{a^2}{m}.$$
(41)
As the semiclassical distribution functions $`\rho _c`$ for the harmonic oscillator coincide with the Wigner functions $`\rho `$ we expect that if $`\lambda 1`$ $`\rho _c`$ does get close to $`\rho `$ for the low lying levels. Thus it may be appropriate to use, instead of the variable $`q`$ and its canonical conjugate momentum $`p`$, the variables which treat harmonic oscillators on the same footing, namely

$$Q=\left[\frac{m\omega _0}{\mathrm{}}\right]^{\frac{1}{2}}x=\left[\frac{2mD}{a^2\mathrm{}^2\lambda }\right]^{\frac{1}{2}}q=\sqrt{\lambda }q,$$
(43)
$$P=\frac{1}{\sqrt{\lambda }}p.$$
(44)
The coordinate $`Q`$ and the momentum $`P`$ have been used in the figures of Sec. IV. We define the dimensionless potential $`V(Q)`$

$$V(Q)=(\mathrm{}\omega _0)^1𝒱(x)=\frac{\lambda }{2}\left(1e^{\frac{Q}{\sqrt{\lambda }}}\right)^2,$$
(45)
which is also used in the figures. For the energy levels there we use similarly

$$E_n=(\mathrm{}\omega _0)^1_n=\left(n+\frac{1}{2}\right)\frac{1}{2\lambda }\left(n+\frac{1}{2}\right)^2.$$
(46)
In the special case of our interest, namely the ground state, $`n=0`$, the wave function is given by

$$\psi _0(\xi )=\left[\frac{2\lambda 1}{\mathrm{\Gamma }(2\lambda )}\right]^{\frac{1}{2}}\xi ^{\lambda \frac{1}{2}}e^{\frac{\xi }{2}}.$$
(47)
In order to calculate the Wigner function replace $`\psi `$ by $`\psi _n(2\lambda e^q)`$ into Eq. (28), obtaining

$$\rho _n^{(\lambda )}(q,p)=(2\pi \mathrm{})^1_{\mathrm{}}^{\mathrm{}}𝑑v\psi _n(2\lambda e^{[q+\frac{v}{2}]})\psi _n^{}(2\lambda e^{[q\frac{v}{2}]})e^{\frac{i}{\mathrm{}}pv}.$$
(48)
Introducing the new integration variable

$$\tau =e^{\frac{v}{2}},$$
(49)
Eq. (48) becomes

$$\rho _n^{(\lambda )}(q,p)=(\pi \mathrm{})^1_0^{\mathrm{}}\psi _n^{}\left(\frac{\xi }{\tau }\right)\psi _n(\xi \tau )\tau ^{2ip\mathrm{}^1}\frac{d\tau }{\tau },$$
(50)
where $`\xi `$ is given by Eq. (35). Substituting $`\psi _n`$ in Eq. (50) by the expression (33) $`\rho _n^{(\lambda )}`$ becomes

$`\rho _n^{(\lambda )}(q,p)=\left[{\displaystyle \frac{\pi \mathrm{}}{2}}\right]^1[N(\lambda ,n)]^2\xi ^{2\lambda 2n1}{\displaystyle \underset{r=0}{\overset{n}{}}}{\displaystyle \underset{s=0}{\overset{n}{}}}b(\lambda ,n,r)b(\lambda ,n,s)\xi ^{r+s}K_{sr+2ip\mathrm{}^1}(\xi ),`$ (51)
where

$$b(\lambda ,n,j)=\left(\begin{array}{c}2\lambda n1\\ nj\end{array}\right)\frac{(1)^j}{j!}$$
(52)
and $`K_\nu (\xi )`$ is defined by

$$K_\nu (\xi )=\frac{1}{2}_0^{\mathrm{}}e^{\frac{\xi }{2}\left(\tau +\frac{1}{\tau }\right)}\tau ^\nu \frac{d\tau }{\tau },$$
(53)
$`\nu `$ being a complex variable. In the particular case in which $`n=0`$ we get

$$\rho _0^{(\lambda )}(q,p)=\left(\frac{\pi \mathrm{}}{2}\right)^1\frac{2\lambda 1}{\mathrm{\Gamma }(2\lambda )}\xi ^{2\lambda 1}K_{2ip\mathrm{}^1}(\xi ).$$
(54)
The numerical method we used to calculate the function $`K_\nu (\xi )`$ will be given in Appendix B.

The symmetries obeyed by the function $`K_\nu `$ are

$`K_\nu ^{}(\xi )={\displaystyle \frac{1}{2}}{\displaystyle _0^{\mathrm{}}}e^{\frac{\xi }{2}(\tau +\frac{1}{\tau })}(\tau ^\nu )^{}{\displaystyle \frac{d\tau }{\tau }}=(K_\nu (\xi ))^{},`$ (55)
$`K_\nu (\xi )={\displaystyle \frac{1}{2}}{\displaystyle _0^{\mathrm{}}}e^{\frac{\xi }{2}(\tau +\frac{1}{\tau })}\tau ^\nu {\displaystyle \frac{d\tau }{\tau }}=K_\nu (\xi ),`$ (56)
where the last step is obtained by making $`\tau \tau ^1`$. From Eqs.(49) and (48) we get

$$K_{a+ib}=K_{a+ib}^{}.$$
(57)
Thus Eq. (51) may be written

$`\rho _n^{(\lambda )}(q,p)=\left[{\displaystyle \frac{\pi \mathrm{}}{2}}\right]^1[N(\lambda ,n)]^2\xi ^{2\lambda 2n1}{\displaystyle \underset{r=0}{\overset{n}{}}}{\displaystyle \underset{s=0}{\overset{n}{}}}b(\lambda ,n,r)b(\lambda ,n,s)\xi ^{r+s}Re\left(K_{sr+2ip\mathrm{}^1}(\xi )\right),`$ (58)
showing explicitly that $`\rho _n^{(\lambda )}`$ is real.

In order to calculate the semiclassical distribution function $`\rho _c(q,p)`$ given by Eq. (10) or (17), we need the solution of the classical equation of motion

$$m\ddot{x}=D\frac{d}{dx}(1e^{ax})^2,$$
(59)
which has already been obtained exactly . Introducing the coordinate $`q=ax`$ and the variable $`\theta `$ given by

$$\theta =\omega _0t,$$
(60)
where $`\omega _0`$ is given by Eq. (41), Eq. (59) becomes

$$2\frac{d^2q}{d\theta ^2}=\frac{d}{dq}(1e^q)^2.$$
(61)
The solution of Eq. (61) for $`ϵ=E/D<1`$, where $`E`$ is the energy associated with the trajectory, is

$$q(E,t)=\mathrm{ln}\left\{\frac{1}{1ϵ}\left[1+\sqrt{ϵ}\mathrm{sin}\sqrt{1ϵ}(\omega _0t\theta _0)\right]\right\}.$$
(62)
The canonical momentum associated with $`q`$ is

$$p=\frac{m\dot{q}}{a^2},$$
(63)
which using Eqs. (62) and (41) gives

$$p(E,t)=\mathrm{}\lambda \frac{\sqrt{ϵ}\sqrt{1ϵ}\mathrm{cos}[\sqrt{1ϵ}(\omega _0t\theta _0)]}{1+\sqrt{ϵ}\mathrm{sin}[\sqrt{1ϵ}(\omega _0t\theta _0)]}.$$
(64)
The period associated with the orbit is

$$T(E)=\frac{2\pi }{\omega _0\sqrt{1ϵ}}.$$
(65)
Analising Eq. (62) one obtains that for $`ED`$ the trajectories are close to those of a harmonic oscillator of frequency $`\omega _0`$.

## IV Results

In this section we present the results of our calculations. We calculated the WDF $`\rho `$ and corresponding SDF $`\rho _c`$ for the ground state of the Morse oscillator choosing for the parameter $`\lambda `$ the values 1,2,4 and 10 corresponding to oscillators with 1,2,4 and 10 levels respectively. In the figures we used throughout the dimensionless coordinates $`Q`$ and the conjugate momenta $`P`$ (in units of $`\mathrm{}`$) defined by Eqs. (43) and (44). In Figs. 1,2,3 and 4 the potentials $`V(Q)`$ defined by Eq. (45) are displayed for $`\lambda =1,2,4`$ and 10 and the energy levels marked. We observe here that the potential given by Eq. (45) is independent of $`\lambda `$ for $`Q^2\lambda `$.

Figs. 5,6,7 and 8 reproduce our calculations of the WDF for $`\lambda =1,2,4`$ and 10 respectively through curves of constant density $`\rho (Q,P)`$. The value of $`\rho `$ varies from $`\rho 0.3`$ to $`\rho =0`$ in the region of phase-space corresponding approximately to the region of the bound classical particles. Then there occur adjacent regions on which $`\rho `$ alternates from negative to positive values and its magnitude decreases as the region gets farther away from the origin.

In the case $`\lambda =1`$ the minimum value of $`\rho `$ is (in units of $`\mathrm{}^1`$) of the order of $`10^2`$, as it can be seen from Fig. 5. This minimum approaches zero as $`\lambda `$ increases, which can be observed from Fig. 7 for $`\lambda =4`$ where the minimum of $`\rho `$ is about $`10^4`$.

The maximum of $`\rho `$, as it can be verified from Figs. 5—8, moves from the point $`(Q,P)=(1.2,0)`$ to $`(Q,P)=(0.3,0)`$ as $`\lambda `$ increases from 1 to 10 and its value increases slightly as $`\lambda `$ increases. Thus the curve for $`\rho =0.3`$, present for $`\lambda 2`$, does not occur for $`\lambda =1`$.

Another feature of the WDF is that the curves of constant density $`\rho `$ become more symmetric with respect to an axis parallel to the $`P`$axis as $`\lambda `$ increases, becoming close to the form of an ellipse. This tallies with the fact that for $`Q\sqrt{\lambda }`$, $`V(Q)`$ is the potential of a harmonic oscillator.

In Figs. 10,11,12 and 13 we present SDF curves for fixed $`\rho _c`$ superposed on WDF curves with $`\rho =\rho _c`$ for comparison. It will be seen that, as $`\lambda `$ increases the SDF approximation improves, which means also that the WDF curves of constant $`\rho `$ become closer to classical trajectories. In Fig. 9 we plot curves of constant $`\rho _c`$ for the case $`\lambda =1`$. For this oscillator our semiclassical approximation is anomalous, as the value of $`\rho _c`$ increases from $`\rho _c=0.145`$ to $`\rho _c=0.179`$ as the classical energy $`E`$ varies from $`E=0`$ to $`E=0.26\mathrm{}\omega _0`$ but decreases as $`E`$ increases further. For the other oscillators $`\rho _c`$ decreases as the energy $`E`$ increases until reaching the value $`\rho _c=0`$. As a consequence for $`\lambda =1`$ one has two closed curves with the same $`\rho _c`$ for $`0.145<\rho _c<0.179`$ whereas for $`\lambda 2`$, for each $`\rho _c`$ one has only one such curve. Also for the WDF there is only one closed curve for each $`\rho `$ from the maximum value of $`\rho `$ up to the curve $`\rho =0`$.

In Fig. 10 we compare the WDF curves of constant $`\rho `$ with the SDF curves for which $`\rho =\rho _c`$ in the case $`\lambda =1`$. One notices that the discrepancies between both curves are very large. In Figs. 11,12 and 13 we make the same comparison for $`\lambda =2,4`$ and 10 respectively.

One finds that for $`\rho <0.05`$ both curves are quite similar but displaced from each other. This displacement becomes less pronounced as $`\rho `$ gets closer to 0.05 so that the best agreement between $`\rho `$ and $`\rho _c`$ is reached for $`\rho 0.05`$. Also as $`\lambda `$ increases the displacement between both curves decreases, as it can be seen comparing the oscillators $`\lambda =4`$ and $`\lambda =10`$. For $`\rho \stackrel{>}{}0.05`$ the curves of fixed $`\rho _c`$ compared with the curves of fixed $`\rho `$ contain a certain amount of distortion which becomes less pronounced as $`\lambda `$ increases.

Finally we observe that curves with $`\rho _c=0.3`$ are absent, the largest values of $`\rho _c`$ being 0.227, 0.271 and 0.299 respectively for $`\lambda =2,4`$ and 10, values which are reached at the origin of phase-space.

## A Eigenfunctions and eigenvalues of the Morse potential

In this Appendix we solve briefly the eigenvalue equation for the Morse potential (Eq. (33) of Sec. III)

$$\left[\frac{d^2}{dz^2}\lambda ^2(1e^z)^2+\lambda ^2ϵ_n\right]\psi _n(z)=0.$$
(A1)
Making the substitution $`y=2\lambda e^z`$, Eq. (A1) is written

$$\left[y^2\frac{d^2}{dy^2}+y\frac{d}{dy}\lambda ^2(1\frac{y}{2\lambda })^2+\lambda ^2ϵ_n\right]\psi _n(y)=0,$$
(A2)
or,

$$\left[\frac{d^2}{dy^2}+\frac{1}{y}\frac{d}{dy}+\frac{(ϵ_n1)\lambda ^2}{y^2}+\frac{\lambda }{y}\frac{1}{4}\right]\psi _n(y)=0.$$
(A3)
Assuming a solution of the form

$$\psi _n(y)=e^{\frac{y}{2}}y^su_n(y),$$
(A4)
and writing $`s^2=(ϵ_n1)\lambda ^2`$, from Eq. (A3) we get

$$yu_n^{\prime \prime }+(2s+1y)u_n^{}+(\lambda s\frac{1}{2})u_n=0,$$
(A5)
where $`u_n^{}=du_n/dy`$ and $`u_n^{\prime \prime }=d^2u_n/dy^2`$. Making $`u_n=_{m=0}^na_my^m`$, Eq. (A5) may be written as

$$\underset{m=0}{\overset{n}{}}[(m+1)(m+2s+1)a_{m+1}+\left(\lambda sm\frac{1}{2}\right)a_m]y^m=0,$$
(A6)
which is the differential equation satisfied by the Laguerre generalized functions. Eq. (A6) is satisfied if

$$a_{m+1}=\frac{s+m\lambda +\frac{1}{2}}{(m+1)(m+2s+1)}a_m,$$
(A7)
which for $`m=0`$ is

$$a_1=\frac{s\lambda +\frac{1}{2}}{2s+1}a_0.$$
(A8)
Since for $`m=n`$, $`a_{m+1}=0`$ and $`a_m0`$, Eq. (A7) gives

$$s=\lambda n\frac{1}{2}.$$
(A9)
Taking into account Eq. (A9) we get

$$ϵ_n=\frac{2}{\lambda }\left(n+\frac{1}{2}\right)+\frac{1}{\lambda ^2}\left(n+\frac{1}{2}\right)^2,$$
(A10)
and the recurrence relation for the coefficients of the power series in Eq. (A7) becomes

$$a_{m+1}=\frac{mn}{(m+1)(2\lambda 2n+m)}a_m.$$
(A11)
Choosing

$$a_0=\frac{(2\lambda n1)!}{n!(2\lambda 2n1)!}N(\lambda ,n),$$
(A12)
we obtain

$$a_m=\frac{(1)^m(2\lambda n1)!N(\lambda ,n)}{m!(nm)!(2\lambda 2nm1)!},$$
(A13)
and we get for $`\psi _n`$ the expression given by Eq. (33) in Sec. III.

## B Numerical Evaluation of the Wigner function

Here we discuss the numerical calculation of the quantity $`K_\nu `$ defined in Eq. (53) of Sec. III,

$$K_\nu (\xi )=\frac{1}{2}_0^{\mathrm{}}e^{\frac{\xi }{2}\left(\tau +\frac{1}{\tau }\right)}\tau ^\nu \frac{d\tau }{\tau }.$$
(B1)
Here $`\nu `$ is complex,

$$\nu =N+2ik,$$
(B2)
where $`N`$ is an integer and $`k`$ is the dimensionless momentum $`p/\mathrm{}`$ and according to Eqs.(35) and (32)

$$\xi =2\lambda e^q,q=ax.$$
(B3)
Making the transformation

$$\tau =e^u,$$
(B4)
and considering that only the real part of $`K_\nu `$ enters into the expression (51) for the Wigner function we get from Eq. (B1), substituting also $`\nu `$ according to Eq. (B2),

$$Re[K_\nu (\xi )]=\frac{1}{2}_{\mathrm{}}^{\mathrm{}}e^{\xi \mathrm{cosh}u}e^{uN}\mathrm{cos}(2ku)𝑑u.$$
(B5)
Let us take initially $`N=0`$, which is the only value needed for the ground state of the Morse oscillator. Consider also $`k0`$ as the case $`k=0`$ is calculated separately. For convenience we introduce the new variable of integration

$$z=2ku,$$
(B6)
and use the fact that the integrand is an even function of $`u`$, obtaining from Eq. (B5)

$$Re[K_\nu (\xi )]=\frac{1}{2k}_0^{\mathrm{}}𝑑ze^{\xi \mathrm{cosh}(\frac{z}{2k})}\mathrm{cos}z.$$
(B7)
The integral in Eq. (B7) is of the form

$$I=_0^{\mathrm{}}f(z)\mathrm{cos}zdz,$$
(B8)
where $`f(z)`$ is a positive decreasing function of $`z`$. In order to avoid numerical cancellations arising from the change of sign of $`\mathrm{cos}(z)`$, we replace this integral by an integral in the interval $`[0,\frac{\pi }{2}]`$ of a series of functions.

We decompose $`I`$ as follows

$$I=I_1+I_2,$$
(B9)
where

$$I_1=_0^{\frac{\pi }{2}}f(z)\mathrm{cos}zdz,$$
(B10)
and

$$I_2=_{\frac{\pi }{2}}^{\mathrm{}}f(z)\mathrm{cos}zdz=_0^{\mathrm{}}f(y+\frac{\pi }{2})\mathrm{sin}ydy,$$
(B11)
where we made the substitution $`z=y+\frac{\pi }{2}`$ in $`I_2`$. The interval of integration of the integral $`I_2`$ is divided in the set of intervals $`[2\pi s,2\pi (s+1)]`$, $`s=0,1,\mathrm{}`$ so that $`I_2`$ is written

$$I_2=\underset{s=0}{\overset{\mathrm{}}{}}_{2\pi s}^{2\pi (s+1)}f(y+\frac{\pi }{2})\mathrm{sin}ydy.$$
(B12)
Making now the substitution

$$z=y2\pi s,$$
(B13)
for the integral in the interval $`[2\pi s,2\pi (s+1)]`$, one obtains for Eq. (B12)

$$I_2=\underset{s=0}{\overset{\mathrm{}}{}}_0^{2\pi }f(z+(2s+\frac{1}{2})\pi )\mathrm{sin}zdz.$$
(B14)
However $`\mathrm{sin}z`$ changes sign in the interval $`[0,2\pi ]`$ which leads to cancellation errors if $`f`$ is slowly varying in the interval. In order to eliminate the oscillation of sign of the integrand we decompose again the integral in the intervals $`[0,\pi ]`$ and $`[\pi ,2\pi ]`$ obtaining

$$I_2=\underset{s=0}{\overset{\mathrm{}}{}}\left\{_0^\pi \left[f\left(z+\left(2s+\frac{1}{2}\right)\pi \right)f\left(z+\left(2s+\frac{3}{2}\right)\pi \right)\right]\mathrm{sin}zdz\right\}.$$
(B15)
As $`f`$ is assumed to be monotonically decreasing each integrand in Eq. (B15) is now always positive. As the domain of integration of the integral $`I_1`$ is the interval $`[0,\frac{\pi }{2}]`$ we make the translation $`v=z\frac{\pi }{2}`$, obtaining

$$I_2=\underset{s=0}{\overset{\mathrm{}}{}}_{\frac{\pi }{2}}^{\frac{\pi }{2}}[f(v+(2s+1)\pi )f(v+(2s+2)\pi )]\mathrm{cos}vdv.$$
(B16)
Making $`vv`$ for the contribution from the interval $`[\frac{\pi }{2},0]`$ in Eq. (B16) and adding the contribution from $`I_1`$ one gets finally

$`{\displaystyle _0^{\mathrm{}}}f(x)\mathrm{cos}xdx`$ $`=`$ $`{\displaystyle _0^{\frac{\pi }{2}}}dx\mathrm{cos}x\{f(x){\displaystyle \underset{s=0}{\overset{\mathrm{}}{}}}[f((2s+1)\pi x)+f((2s+1)\pi +x)`$ (B17)
$``$ $`f((2s+2)\pi x)f((2s+2)\pi +x)]\}.`$ (B18)
For $`f(x)`$ monotonic each term of the series contributes with the same sign, however errors may arise from the subtraction of the sum of the series from $`f(x)`$.

In the general case in which $`N0`$ in Eq. (B5), the function $`f(z)`$ in Eq. (B8) becomes

$$f_N(z)=\frac{1}{2k}e^{\xi \mathrm{cosh}\frac{z}{2k}}\mathrm{cosh}\left(\frac{Nz}{2k}\right),z0,N=0,1,\mathrm{}$$
(B19)
The function $`f_N(z)`$ has $`M(MN)`$ extremes at the real positive points $`z_1^{(N)},z_2^{(N)},\mathrm{},z_M^{(N)}`$, labeled in the order of increasing magnitude. This function decreases monotonically for $`z>z_M^{(N)}`$ and it may still be useful to apply the decomposition (B18) in order to eliminate errors due to cancellations of contributions of opposite sign from the integrand $`f_N(z)\mathrm{cos}z`$. In the general case only the sign of the first $`s_N`$ terms of the series in Eq. (B18) may oscillate, where $`s_N`$ is given roughly by the smallest integer satisfying $`(2s_N+1)\pi >z_M^{(N)}`$. For $`s>s_N`$ the sign of the terms of the series in Eq. (B18) is always positive.

In fact, one may determine the extremes of $`f_N(z)`$ by expressing $`f_N(z)`$ in terms of $`\mathrm{cosh}(z/2k)`$ and applying the condition $`f_N(z)/z=0`$. By making this substitution one gets for Eq. (B19) the expression

$$f_N(z)=\frac{1}{2k}e^{\xi y}\underset{j=0}{\overset{[N/2]}{}}A_j(1)^jy^{N2j},y=\mathrm{cosh}(\frac{z}{2k}),$$
(B20)
where $`[N/2]`$ denotes the largest integer contained in $`N/2`$ and

$$A_j=\underset{i=j}{\overset{[N/2]}{}}\left(\begin{array}{c}N\\ 2i\end{array}\right)\left(\begin{array}{c}i\\ j\end{array}\right).$$
(B21)
Thus we get the extrema of $`f_N(z)`$ as the roots of a polynomial of degree $`N`$.

For $`N=1`$ the maximum occurs at

$$z_1^{(1)}=2k\mathrm{cosh}^1\left(\frac{1}{\xi }\right).$$
(B22)
For $`N=2`$ the maximum will be at

$$z_1^{(2)}=2k\mathrm{cosh}^1\left(\frac{2}{\xi }+\sqrt{\frac{4}{\xi ^2}+2}\right).$$
(B23)